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)