sysdeps/libm-ieee754/s_csqrt.c

   1 /* Complex square root of double value.
   2    Copyright (C) 1997, 1998 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           double d, r, s;
  88
  89           d = __ieee754_hypot (__real__ x, __imag__ x);
  90           /* Use the identity   2  Re res  Im res = Im x
  91              to avoid cancellation error in  d +/- Re x.  */
  92           if (__real__ x > 0)
  93             {
  94               r = __ieee754_sqrt (0.5 * d + 0.5 * __real__ x);
  95               s = (0.5 * __imag__ x) / r;
  96             }
  97           else
  98             {
  99               s = __ieee754_sqrt (0.5 * d - 0.5 * __real__ x);
 100               r = fabs ((0.5 * __imag__ x) / s);
 101             }
 102
 103           __real__ res = r;
 104           __imag__ res = __copysign (s, __imag__ x);
 105         }
 106     }
 107
 108   return res;
 109 }
 110 weak_alias (__csqrt, csqrt)
 111 #ifdef NO_LONG_DOUBLE
 112 strong_alias (__csqrt, __csqrtl)
 113 weak_alias (__csqrt, csqrtl)
 114 #endif