sysdeps/libm-ieee754/s_csqrtf.c

   1 /* Complex square root of float 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__ float
  29 __csqrtf (__complex__ float x)
  30 {
  31   __complex__ float 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_VALF;
  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 ? __nanf ("") : 0;
  47               __imag__ res = __copysignf (HUGE_VALF, __imag__ x);
  48             }
  49           else
  50             {
  51               __real__ res = __real__ x;
  52               __imag__ res = (icls == FP_NAN
  53                               ? __nanf ("") : __copysignf (0.0, __imag__ x));
  54             }
  55         }
  56       else
  57         {
  58           __real__ res = __nanf ("");
  59           __imag__ res = __nanf ("");
  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 = __copysignf (__ieee754_sqrtf (-__real__ x),
  70                                           __imag__ x);
  71             }
  72           else
  73             {
  74               __real__ res = fabsf (__ieee754_sqrtf (__real__ x));
  75               __imag__ res = __copysignf (0.0, __imag__ x);
  76             }
  77         }
  78       else if (rcls == FP_ZERO)
  79         {
  80           float r = __ieee754_sqrtf (0.5 * fabsf (__imag__ x));
  81
  82           __real__ res = __copysignf (r, __imag__ x);
  83           __imag__ res = r;
  84         }
  85       else
  86         {
  87 #if 0 /* FIXME: this is broken. */
  88           __complex__ float q;
  89           float t, r;
  90
  91           if (fabsf (__imag__ x) < 2.0e-4 * fabsf (__real__ x))
  92             t = 0.25 * __imag__ x * (__imag__ x / __real__ x);
  93           else
  94             t = 0.5 * (__ieee754_hypotf (__real__ x, __imag__ x) - __real__ x);
  95
  96           r = __ieee754_sqrtf (t);
  97
  98           __real__ q = __imag__ x / (2.0 * r);
  99           __imag__ q = r;
 100
 101           /* Heron iteration in complex arithmetic.  */
 102           res = 0.5 * (q + q / x);
 103 #else
 104           float d, imag;
 105
 106           d = __ieee754_hypotf (__real__ x, __imag__ x);
 107           imag = __ieee754_sqrtf (0.5 * (d - __real__ x));
 108
 109           __real__ res = __ieee754_sqrtf (0.5 * (d + __real__ x));
 110           __imag__ res = __copysignf (imag, __imag__ x);
 111 #endif
 112         }
 113     }
 114
 115   return res;
 116 }
 117 weak_alias (__csqrtf, csqrtf)