math/s_csqrt_template.c

   1 /* Complex square root of a float type.
   2    Copyright (C) 1997-2023 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 #include <complex.h>
  20 #include <math.h>
  21 #include <math_private.h>
  22 #include <math-underflow.h>
  23 #include <float.h>
  24
  25 CFLOAT
  26 M_DECL_FUNC (__csqrt) (CFLOAT x)
  27 {
  28   CFLOAT res;
  29   int rcls = fpclassify (__real__ x);
  30   int icls = fpclassify (__imag__ x);
  31
  32   if (__glibc_unlikely (rcls <= FP_INFINITE || icls <= FP_INFINITE))
  33     {
  34       if (icls == FP_INFINITE)
  35         {
  36           __real__ res = M_HUGE_VAL;
  37           __imag__ res = __imag__ x;
  38         }
  39       else if (rcls == FP_INFINITE)
  40         {
  41           if (__real__ x < 0)
  42             {
  43               __real__ res = icls == FP_NAN ? M_NAN : 0;
  44               __imag__ res = M_COPYSIGN (M_HUGE_VAL, __imag__ x);
  45             }
  46           else
  47             {
  48               __real__ res = __real__ x;
  49               __imag__ res = (icls == FP_NAN
  50                               ? M_NAN : M_COPYSIGN (0, __imag__ x));
  51             }
  52         }
  53       else
  54         {
  55           __real__ res = M_NAN;
  56           __imag__ res = M_NAN;
  57         }
  58     }
  59   else
  60     {
  61       if (__glibc_unlikely (icls == FP_ZERO))
  62         {
  63           if (__real__ x < 0)
  64             {
  65               __real__ res = 0;
  66               __imag__ res = M_COPYSIGN (M_SQRT (-__real__ x), __imag__ x);
  67             }
  68           else
  69             {
  70               __real__ res = M_FABS (M_SQRT (__real__ x));
  71               __imag__ res = M_COPYSIGN (0, __imag__ x);
  72             }
  73         }
  74       else if (__glibc_unlikely (rcls == FP_ZERO))
  75         {
  76           FLOAT r;
  77           if (M_FABS (__imag__ x) >= 2 * M_MIN)
  78             r = M_SQRT (M_LIT (0.5) * M_FABS (__imag__ x));
  79           else
  80             r = M_LIT (0.5) * M_SQRT (2 * M_FABS (__imag__ x));
  81
  82           __real__ res = r;
  83           __imag__ res = M_COPYSIGN (r, __imag__ x);
  84         }
  85       else
  86         {
  87           FLOAT d, r, s;
  88           int scale = 0;
  89
  90           if (M_FABS (__real__ x) > M_MAX / 4)
  91             {
  92               scale = 1;
  93               __real__ x = M_SCALBN (__real__ x, -2 * scale);
  94               __imag__ x = M_SCALBN (__imag__ x, -2 * scale);
  95             }
  96           else if (M_FABS (__imag__ x) > M_MAX / 4)
  97             {
  98               scale = 1;
  99               if (M_FABS (__real__ x) >= 4 * M_MIN)
 100                 __real__ x = M_SCALBN (__real__ x, -2 * scale);
 101               else
 102                 __real__ x = 0;
 103               __imag__ x = M_SCALBN (__imag__ x, -2 * scale);
 104             }
 105           else if (M_FABS (__real__ x) < 2 * M_MIN
 106                    && M_FABS (__imag__ x) < 2 * M_MIN)
 107             {
 108               scale = -((M_MANT_DIG + 1) / 2);
 109               __real__ x = M_SCALBN (__real__ x, -2 * scale);
 110               __imag__ x = M_SCALBN (__imag__ x, -2 * scale);
 111             }
 112
 113           d = M_HYPOT (__real__ x, __imag__ x);
 114           /* Use the identity   2  Re res  Im res = Im x
 115              to avoid cancellation error in  d +/- Re x.  */
 116           if (__real__ x > 0)
 117             {
 118               r = M_SQRT (M_LIT (0.5) * (d + __real__ x));
 119               if (scale == 1 && M_FABS (__imag__ x) < 1)
 120                 {
 121                   /* Avoid possible intermediate underflow.  */
 122                   s = __imag__ x / r;
 123                   r = M_SCALBN (r, scale);
 124                   scale = 0;
 125                 }
 126               else
 127                 s = M_LIT (0.5) * (__imag__ x / r);
 128             }
 129           else
 130             {
 131               s = M_SQRT (M_LIT (0.5) * (d - __real__ x));
 132               if (scale == 1 && M_FABS (__imag__ x) < 1)
 133                 {
 134                   /* Avoid possible intermediate underflow.  */
 135                   r = M_FABS (__imag__ x / s);
 136                   s = M_SCALBN (s, scale);
 137                   scale = 0;
 138                 }
 139               else
 140                 r = M_FABS (M_LIT (0.5) * (__imag__ x / s));
 141             }
 142
 143           if (scale)
 144             {
 145               r = M_SCALBN (r, scale);
 146               s = M_SCALBN (s, scale);
 147             }
 148
 149           math_check_force_underflow (r);
 150           math_check_force_underflow (s);
 151
 152           __real__ res = r;
 153           __imag__ res = M_COPYSIGN (s, __imag__ x);
 154         }
 155     }
 156
 157   return res;
 158 }
 159 declare_mgen_alias (__csqrt, csqrt)