libquadmath/math/csqrtq.c

   1 /* Complex square root of a float type.
   2    Copyright (C) 1997-2018 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 __complex128
  24 csqrtq (__complex128 x)
  25 {
  26   __complex128 res;
  27   int rcls = fpclassifyq (__real__ x);
  28   int icls = fpclassifyq (__imag__ x);
  29
  30   if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
  31     {
  32       if (icls == QUADFP_INFINITE)
  33         {
  34           __real__ res = HUGE_VALQ;
  35           __imag__ res = __imag__ x;
  36         }
  37       else if (rcls == QUADFP_INFINITE)
  38         {
  39           if (__real__ x < 0)
  40             {
  41               __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
  42               __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
  43             }
  44           else
  45             {
  46               __real__ res = __real__ x;
  47               __imag__ res = (icls == QUADFP_NAN
  48                               ? nanq ("") : copysignq (0, __imag__ x));
  49             }
  50         }
  51       else
  52         {
  53           __real__ res = nanq ("");
  54           __imag__ res = nanq ("");
  55         }
  56     }
  57   else
  58     {
  59       if (__glibc_unlikely (icls == QUADFP_ZERO))
  60         {
  61           if (__real__ x < 0)
  62             {
  63               __real__ res = 0;
  64               __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x);
  65             }
  66           else
  67             {
  68               __real__ res = fabsq (sqrtq (__real__ x));
  69               __imag__ res = copysignq (0, __imag__ x);
  70             }
  71         }
  72       else if (__glibc_unlikely (rcls == QUADFP_ZERO))
  73         {
  74           __float128 r;
  75           if (fabsq (__imag__ x) >= 2 * FLT128_MIN)
  76             r = sqrtq (0.5Q * fabsq (__imag__ x));
  77           else
  78             r = 0.5Q * sqrtq (2 * fabsq (__imag__ x));
  79
  80           __real__ res = r;
  81           __imag__ res = copysignq (r, __imag__ x);
  82         }
  83       else
  84         {
  85           __float128 d, r, s;
  86           int scale = 0;
  87
  88           if (fabsq (__real__ x) > FLT128_MAX / 4)
  89             {
  90               scale = 1;
  91               __real__ x = scalbnq (__real__ x, -2 * scale);
  92               __imag__ x = scalbnq (__imag__ x, -2 * scale);
  93             }
  94           else if (fabsq (__imag__ x) > FLT128_MAX / 4)
  95             {
  96               scale = 1;
  97               if (fabsq (__real__ x) >= 4 * FLT128_MIN)
  98                 __real__ x = scalbnq (__real__ x, -2 * scale);
  99               else
 100                 __real__ x = 0;
 101               __imag__ x = scalbnq (__imag__ x, -2 * scale);
 102             }
 103           else if (fabsq (__real__ x) < 2 * FLT128_MIN
 104                    && fabsq (__imag__ x) < 2 * FLT128_MIN)
 105             {
 106               scale = -((FLT128_MANT_DIG + 1) / 2);
 107               __real__ x = scalbnq (__real__ x, -2 * scale);
 108               __imag__ x = scalbnq (__imag__ x, -2 * scale);
 109             }
 110
 111           d = hypotq (__real__ x, __imag__ x);
 112           /* Use the identity   2  Re res  Im res = Im x
 113              to avoid cancellation error in  d +/- Re x.  */
 114           if (__real__ x > 0)
 115             {
 116               r = sqrtq (0.5Q * (d + __real__ x));
 117               if (scale == 1 && fabsq (__imag__ x) < 1)
 118                 {
 119                   /* Avoid possible intermediate underflow.  */
 120                   s = __imag__ x / r;
 121                   r = scalbnq (r, scale);
 122                   scale = 0;
 123                 }
 124               else
 125                 s = 0.5Q * (__imag__ x / r);
 126             }
 127           else
 128             {
 129               s = sqrtq (0.5Q * (d - __real__ x));
 130               if (scale == 1 && fabsq (__imag__ x) < 1)
 131                 {
 132                   /* Avoid possible intermediate underflow.  */
 133                   r = fabsq (__imag__ x / s);
 134                   s = scalbnq (s, scale);
 135                   scale = 0;
 136                 }
 137               else
 138                 r = fabsq (0.5Q * (__imag__ x / s));
 139             }
 140
 141           if (scale)
 142             {
 143               r = scalbnq (r, scale);
 144               s = scalbnq (s, scale);
 145             }
 146
 147           math_check_force_underflow (r);
 148           math_check_force_underflow (s);
 149
 150           __real__ res = r;
 151           __imag__ res = copysignq (s, __imag__ x);
 152         }
 153     }
 154
 155   return res;
 156 }