Back to index

glibc  2.9
s_csqrtf.c
Go to the documentation of this file.
00001 /* Complex square root of float value.
00002    Copyright (C) 1997, 1998, 2005 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
00005    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
00006 
00007    The GNU C Library is free software; you can redistribute it and/or
00008    modify it under the terms of the GNU Lesser General Public
00009    License as published by the Free Software Foundation; either
00010    version 2.1 of the License, or (at your option) any later version.
00011 
00012    The GNU C Library is distributed in the hope that it will be useful,
00013    but WITHOUT ANY WARRANTY; without even the implied warranty of
00014    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015    Lesser General Public License for more details.
00016 
00017    You should have received a copy of the GNU Lesser General Public
00018    License along with the GNU C Library; if not, write to the Free
00019    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00020    02111-1307 USA.  */
00021 
00022 #include <complex.h>
00023 #include <math.h>
00024 
00025 #include "math_private.h"
00026 
00027 
00028 __complex__ float
00029 __csqrtf (__complex__ float x)
00030 {
00031   __complex__ float res;
00032   int rcls = fpclassify (__real__ x);
00033   int icls = fpclassify (__imag__ x);
00034 
00035   if (rcls <= FP_INFINITE || icls <= FP_INFINITE)
00036     {
00037       if (icls == FP_INFINITE)
00038        {
00039          __real__ res = HUGE_VALF;
00040          __imag__ res = __imag__ x;
00041        }
00042       else if (rcls == FP_INFINITE)
00043        {
00044          if (__real__ x < 0.0)
00045            {
00046              __real__ res = icls == FP_NAN ? __nanf ("") : 0;
00047              __imag__ res = __copysignf (HUGE_VALF, __imag__ x);
00048            }
00049          else
00050            {
00051              __real__ res = __real__ x;
00052              __imag__ res = (icls == FP_NAN
00053                            ? __nanf ("") : __copysignf (0.0, __imag__ x));
00054            }
00055        }
00056       else
00057        {
00058          __real__ res = __nanf ("");
00059          __imag__ res = __nanf ("");
00060        }
00061     }
00062   else
00063     {
00064       if (icls == FP_ZERO)
00065        {
00066          if (__real__ x < 0.0)
00067            {
00068              __real__ res = 0.0;
00069              __imag__ res = __copysignf (__ieee754_sqrtf (-__real__ x),
00070                                      __imag__ x);
00071            }
00072          else
00073            {
00074              __real__ res = fabsf (__ieee754_sqrtf (__real__ x));
00075              __imag__ res = __copysignf (0.0, __imag__ x);
00076            }
00077        }
00078       else if (rcls == FP_ZERO)
00079        {
00080          float r = __ieee754_sqrtf (0.5 * fabsf (__imag__ x));
00081 
00082          __real__ res = r;
00083          __imag__ res = __copysignf (r, __imag__ x);
00084        }
00085       else
00086        {
00087          float d, r, s;
00088 
00089          d = __ieee754_hypotf (__real__ x, __imag__ x);
00090          /* Use the identity   2  Re res  Im res = Im x
00091             to avoid cancellation error in  d +/- Re x.  */
00092          if (__real__ x > 0)
00093            {
00094              r = __ieee754_sqrtf (0.5f * d + 0.5f * __real__ x);
00095              s = (0.5f * __imag__ x) / r;
00096            }
00097          else
00098            {
00099              s = __ieee754_sqrtf (0.5f * d - 0.5f * __real__ x);
00100              r = fabsf ((0.5f * __imag__ x) / s);
00101            }
00102 
00103          __real__ res = r;
00104          __imag__ res = __copysignf (s, __imag__ x);
00105        }
00106     }
00107 
00108   return res;
00109 }
00110 #ifndef __csqrtf
00111 weak_alias (__csqrtf, csqrtf)
00112 #endif