libquadmath/math/csqrtq.c

   1 /* Complex square root of __float128 value.
   2    Copyright (C) 1997-2012 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 Lesser General Public
   9    License as published by the Free Software Foundation; either
  10    version 2.1 of the 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    Lesser General Public License for more details.
  16
  17    You should have received a copy of the GNU Lesser General Public
  18    License along with the GNU C Library; if not, see
  19    <http://www.gnu.org/licenses/>.  */
  20
  21 #include "quadmath-imp.h"
  22
  23 #ifdef HAVE_FENV_H
  24 # include <fenv.h>
  25 #endif
  26
  27
  28 __complex128
  29 csqrtq (__complex128 x)
  30 {
  31   __complex128 res;
  32   int rcls = fpclassifyq (__real__ x);
  33   int icls = fpclassifyq (__imag__ x);
  34
  35   if (__builtin_expect (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE, 0))
  36     {
  37       if (icls == QUADFP_INFINITE)
  38         {
  39           __real__ res = HUGE_VALQ;
  40           __imag__ res = __imag__ x;
  41         }
  42       else if (rcls == QUADFP_INFINITE)
  43         {
  44           if (__real__ x < 0.0Q)
  45             {
  46               __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
  47               __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
  48             }
  49           else
  50             {
  51               __real__ res = __real__ x;
  52               __imag__ res = (icls == QUADFP_NAN
  53                               ? nanq ("") : copysignq (0.0Q, __imag__ x));
  54             }
  55         }
  56       else
  57         {
  58           __real__ res = nanq ("");
  59           __imag__ res = nanq ("");
  60         }
  61     }
  62   else
  63     {
  64       if (__builtin_expect (icls == QUADFP_ZERO, 0))
  65         {
  66           if (__real__ x < 0.0Q)
  67             {
  68               __real__ res = 0.0Q;
  69               __imag__ res = copysignq (sqrtq (-__real__ x),
  70                                         __imag__ x);
  71             }
  72           else
  73             {
  74               __real__ res = fabsq (sqrtq (__real__ x));
  75               __imag__ res = copysignq (0.0Q, __imag__ x);
  76             }
  77         }
  78       else if (__builtin_expect (rcls == QUADFP_ZERO, 0))
  79         {
  80           __float128 r;
  81           if (fabsq (__imag__ x) >= 2.0Q * FLT128_MIN)
  82             r = sqrtq (0.5Q * fabsq (__imag__ x));
  83           else
  84             r = 0.5Q * sqrtq (2.0Q * fabsq (__imag__ x));
  85
  86           __real__ res = r;
  87           __imag__ res = copysignq (r, __imag__ x);
  88         }
  89       else
  90         {
  91           __float128 d, r, s;
  92           int scale = 0;
  93
  94           if (fabsq (__real__ x) > FLT128_MAX / 4.0Q)
  95             {
  96               scale = 1;
  97               __real__ x = scalbnq (__real__ x, -2 * scale);
  98               __imag__ x = scalbnq (__imag__ x, -2 * scale);
  99             }
 100           else if (fabsq (__imag__ x) > FLT128_MAX / 4.0Q)
 101             {
 102               scale = 1;
 103               if (fabsq (__real__ x) >= 4.0Q * FLT128_MIN)
 104                 __real__ x = scalbnq (__real__ x, -2 * scale);
 105               else
 106                 __real__ x = 0.0Q;
 107               __imag__ x = scalbnq (__imag__ x, -2 * scale);
 108             }
 109           else if (fabsq (__real__ x) < FLT128_MIN
 110                    && fabsq (__imag__ x) < FLT128_MIN)
 111             {
 112               scale = -(FLT128_MANT_DIG / 2);
 113               __real__ x = scalbnq (__real__ x, -2 * scale);
 114               __imag__ x = scalbnq (__imag__ x, -2 * scale);
 115             }
 116
 117           d = hypotq (__real__ x, __imag__ x);
 118           /* Use the identity   2  Re res  Im res = Im x
 119              to avoid cancellation error in  d +/- Re x.  */
 120           if (__real__ x > 0)
 121             {
 122               r = sqrtq (0.5Q * (d + __real__ x));
 123               s = 0.5Q * (__imag__ x / r);
 124             }
 125           else
 126             {
 127               s = sqrtq (0.5Q * (d - __real__ x));
 128               r = fabsq (0.5Q * (__imag__ x / s));
 129             }
 130
 131           if (scale)
 132             {
 133               r = scalbnq (r, scale);
 134               s = scalbnq (s, scale);
 135             }
 136
 137           __real__ res = r;
 138           __imag__ res = copysignq (s, __imag__ x);
 139         }
 140     }
 141
 142   return res;
 143 }