sysdeps/libm-ieee754/s_csqrt.c

   1 /* Complex square root of double value.
   2    Copyright (C) 1997 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
   5    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
   6
   7    The GNU C Library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Library General Public License as
   9    published by the Free Software Foundation; either version 2 of the
  10    License, or (at your option) any later version.
  11
  12    The GNU C Library is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    Library General Public License for more details.
  16
  17    You should have received a copy of the GNU Library General Public
  18    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  19    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20    Boston, MA 02111-1307, USA.  */
  21
  22 #include <complex.h>
  23 #include <math.h>
  24
  25 #include "math_private.h"
  26
  27
  28 __complex__ double
  29 __csqrt (__complex__ double x)
  30 {
  31   __complex__ double res;
  32   int rcls = fpclassify (__real__ x);
  33   int icls = fpclassify (__imag__ x);
  34
  35   if (rcls <= FP_INFINITE || icls <= FP_INFINITE)
  36     {
  37       if (icls == FP_INFINITE)
  38         {
  39           __real__ res = HUGE_VAL;
  40           __imag__ res = __imag__ x;
  41         }
  42       else if (rcls == FP_INFINITE)
  43         {
  44           if (__real__ x < 0.0)
  45             {
  46               __real__ res = icls == FP_NAN ? __nan ("") : 0;
  47               __imag__ res = __copysign (HUGE_VAL, __imag__ x);
  48             }
  49           else
  50             {
  51               __real__ res = __real__ x;
  52               __imag__ res = (icls == FP_NAN
  53                               ? __nan ("") : __copysign (0.0, __imag__ x));
  54             }
  55         }
  56       else
  57         {
  58           __real__ res = __nan ("");
  59           __imag__ res = __nan ("");
  60         }
  61     }
  62   else
  63     {
  64       if (icls == FP_ZERO)
  65         {
  66           if (__real__ x < 0.0)
  67             {
  68               __real__ res = 0.0;
  69               __imag__ res = __copysign (__ieee754_sqrt (-__real__ x),
  70                                          __imag__ x);
  71             }
  72           else
  73             {
  74               __real__ res = fabs (__ieee754_sqrt (__real__ x));
  75               __imag__ res = __copysign (0.0, __imag__ x);
  76             }
  77         }
  78       else if (rcls == FP_ZERO)
  79         {
  80           double r = __ieee754_sqrt (0.5 * fabs (__imag__ x));
  81
  82           __real__ res = __copysign (r, __imag__ x);
  83           __imag__ res = r;
  84         }
  85       else
  86         {
  87           __complex__ double q;
  88           double t, r;
  89
  90           if (fabs (__imag__ x) < 2.0e-4 * fabs (__real__ x))
  91             t = 0.25 * __imag__ x * (__imag__ x / __real__ x);
  92           else
  93             t = 0.5 * (__ieee754_hypot (__real__ x, __imag__ x) - __real__ x);
  94
  95           r = __ieee754_sqrt (t);
  96
  97           __real__ q = __imag__ x / (2.0 * r);
  98           __imag__ q = r;
  99
 100           /* Heron iteration in complex arithmetic.  */
 101           res = 0.5 * (q + q / x);
 102         }
 103     }
 104
 105   return res;
 106 }
 107 weak_alias (__csqrt, csqrt)
 108 #ifdef NO_LONG_DOUBLE
 109 strong_alias (__csqrt, __csqrtl)
 110 weak_alias (__csqrt, csqrtl)
 111 #endif