Back to index

glibc  2.9
s_csqrtl.c
Go to the documentation of this file.
00001 /* Complex square root of long double 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__ long double
00029 __csqrtl (__complex__ long double x)
00030 {
00031   __complex__ long double 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_VALL;
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 ? __nanl ("") : 0;
00047              __imag__ res = __copysignl (HUGE_VALL, __imag__ x);
00048            }
00049          else
00050            {
00051              __real__ res = __real__ x;
00052              __imag__ res = (icls == FP_NAN
00053                            ? __nanl ("") : __copysignl (0.0, __imag__ x));
00054            }
00055        }
00056       else
00057        {
00058          __real__ res = __nanl ("");
00059          __imag__ res = __nanl ("");
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 = __copysignl (__ieee754_sqrtl (-__real__ x),
00070                                      __imag__ x);
00071            }
00072          else
00073            {
00074              __real__ res = fabsl (__ieee754_sqrtl (__real__ x));
00075              __imag__ res = __copysignl (0.0, __imag__ x);
00076            }
00077        }
00078       else if (rcls == FP_ZERO)
00079        {
00080          long double r = __ieee754_sqrtl (0.5 * fabsl (__imag__ x));
00081 
00082          __real__ res = r;
00083          __imag__ res = __copysignl (r, __imag__ x);
00084        }
00085       else
00086        {
00087          long double d, r, s;
00088 
00089          d = __ieee754_hypotl (__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_sqrtl (0.5L * d + 0.5L * __real__ x);
00095              s = (0.5L * __imag__ x) / r;
00096            }
00097          else
00098            {
00099              s = __ieee754_sqrtl (0.5L * d - 0.5L * __real__ x);
00100              r = fabsl ((0.5L * __imag__ x) / s);
00101            }
00102 
00103          __real__ res = r;
00104          __imag__ res = __copysignl (s, __imag__ x);
00105        }
00106     }
00107 
00108   return res;
00109 }
00110 weak_alias (__csqrtl, csqrtl)