1 /* Complex square root of a float type.
2 Copyright (C) 1997-2020 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.
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.
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.
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 <https://www.gnu.org/licenses/>. */
23 #include <math_private.h>
24 #include <math-underflow.h>
28 M_DECL_FUNC (__csqrt
) (CFLOAT x
)
31 int rcls
= fpclassify (__real__ x
);
32 int icls
= fpclassify (__imag__ x
);
34 if (__glibc_unlikely (rcls
<= FP_INFINITE
|| icls
<= FP_INFINITE
))
36 if (icls
== FP_INFINITE
)
38 __real__ res
= M_HUGE_VAL
;
39 __imag__ res
= __imag__ x
;
41 else if (rcls
== FP_INFINITE
)
45 __real__ res
= icls
== FP_NAN
? M_NAN
: 0;
46 __imag__ res
= M_COPYSIGN (M_HUGE_VAL
, __imag__ x
);
50 __real__ res
= __real__ x
;
51 __imag__ res
= (icls
== FP_NAN
52 ? M_NAN
: M_COPYSIGN (0, __imag__ x
));
63 if (__glibc_unlikely (icls
== FP_ZERO
))
68 __imag__ res
= M_COPYSIGN (M_SQRT (-__real__ x
), __imag__ x
);
72 __real__ res
= M_FABS (M_SQRT (__real__ x
));
73 __imag__ res
= M_COPYSIGN (0, __imag__ x
);
76 else if (__glibc_unlikely (rcls
== FP_ZERO
))
79 if (M_FABS (__imag__ x
) >= 2 * M_MIN
)
80 r
= M_SQRT (M_LIT (0.5) * M_FABS (__imag__ x
));
82 r
= M_LIT (0.5) * M_SQRT (2 * M_FABS (__imag__ x
));
85 __imag__ res
= M_COPYSIGN (r
, __imag__ x
);
92 if (M_FABS (__real__ x
) > M_MAX
/ 4)
95 __real__ x
= M_SCALBN (__real__ x
, -2 * scale
);
96 __imag__ x
= M_SCALBN (__imag__ x
, -2 * scale
);
98 else if (M_FABS (__imag__ x
) > M_MAX
/ 4)
101 if (M_FABS (__real__ x
) >= 4 * M_MIN
)
102 __real__ x
= M_SCALBN (__real__ x
, -2 * scale
);
105 __imag__ x
= M_SCALBN (__imag__ x
, -2 * scale
);
107 else if (M_FABS (__real__ x
) < 2 * M_MIN
108 && M_FABS (__imag__ x
) < 2 * M_MIN
)
110 scale
= -((M_MANT_DIG
+ 1) / 2);
111 __real__ x
= M_SCALBN (__real__ x
, -2 * scale
);
112 __imag__ x
= M_SCALBN (__imag__ x
, -2 * scale
);
115 d
= M_HYPOT (__real__ x
, __imag__ x
);
116 /* Use the identity 2 Re res Im res = Im x
117 to avoid cancellation error in d +/- Re x. */
120 r
= M_SQRT (M_LIT (0.5) * (d
+ __real__ x
));
121 if (scale
== 1 && M_FABS (__imag__ x
) < 1)
123 /* Avoid possible intermediate underflow. */
125 r
= M_SCALBN (r
, scale
);
129 s
= M_LIT (0.5) * (__imag__ x
/ r
);
133 s
= M_SQRT (M_LIT (0.5) * (d
- __real__ x
));
134 if (scale
== 1 && M_FABS (__imag__ x
) < 1)
136 /* Avoid possible intermediate underflow. */
137 r
= M_FABS (__imag__ x
/ s
);
138 s
= M_SCALBN (s
, scale
);
142 r
= M_FABS (M_LIT (0.5) * (__imag__ x
/ s
));
147 r
= M_SCALBN (r
, scale
);
148 s
= M_SCALBN (s
, scale
);
151 math_check_force_underflow (r
);
152 math_check_force_underflow (s
);
155 __imag__ res
= M_COPYSIGN (s
, __imag__ x
);
161 declare_mgen_alias (__csqrt
, csqrt
)