1 /* Complex square root of double value.
2 Copyright (C) 1997-2015 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 <http://www.gnu.org/licenses/>. */
23 #include <math_private.h>
27 __csqrt (__complex__
double x
)
29 __complex__
double res
;
30 int rcls
= fpclassify (__real__ x
);
31 int icls
= fpclassify (__imag__ x
);
33 if (__glibc_unlikely (rcls
<= FP_INFINITE
|| icls
<= FP_INFINITE
))
35 if (icls
== FP_INFINITE
)
37 __real__ res
= HUGE_VAL
;
38 __imag__ res
= __imag__ x
;
40 else if (rcls
== FP_INFINITE
)
44 __real__ res
= icls
== FP_NAN
? __nan ("") : 0;
45 __imag__ res
= __copysign (HUGE_VAL
, __imag__ x
);
49 __real__ res
= __real__ x
;
50 __imag__ res
= (icls
== FP_NAN
51 ? __nan ("") : __copysign (0.0, __imag__ x
));
56 __real__ res
= __nan ("");
57 __imag__ res
= __nan ("");
62 if (__glibc_unlikely (icls
== FP_ZERO
))
67 __imag__ res
= __copysign (__ieee754_sqrt (-__real__ x
),
72 __real__ res
= fabs (__ieee754_sqrt (__real__ x
));
73 __imag__ res
= __copysign (0.0, __imag__ x
);
76 else if (__glibc_unlikely (rcls
== FP_ZERO
))
79 if (fabs (__imag__ x
) >= 2.0 * DBL_MIN
)
80 r
= __ieee754_sqrt (0.5 * fabs (__imag__ x
));
82 r
= 0.5 * __ieee754_sqrt (2.0 * fabs (__imag__ x
));
85 __imag__ res
= __copysign (r
, __imag__ x
);
92 if (fabs (__real__ x
) > DBL_MAX
/ 4.0)
95 __real__ x
= __scalbn (__real__ x
, -2 * scale
);
96 __imag__ x
= __scalbn (__imag__ x
, -2 * scale
);
98 else if (fabs (__imag__ x
) > DBL_MAX
/ 4.0)
101 if (fabs (__real__ x
) >= 4.0 * DBL_MIN
)
102 __real__ x
= __scalbn (__real__ x
, -2 * scale
);
105 __imag__ x
= __scalbn (__imag__ x
, -2 * scale
);
107 else if (fabs (__real__ x
) < 2.0 * DBL_MIN
108 && fabs (__imag__ x
) < 2.0 * DBL_MIN
)
110 scale
= -((DBL_MANT_DIG
+ 1) / 2);
111 __real__ x
= __scalbn (__real__ x
, -2 * scale
);
112 __imag__ x
= __scalbn (__imag__ x
, -2 * scale
);
115 d
= __ieee754_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
= __ieee754_sqrt (0.5 * (d
+ __real__ x
));
121 if (scale
== 1 && fabs (__imag__ x
) < 1.0)
123 /* Avoid possible intermediate underflow. */
125 r
= __scalbn (r
, scale
);
129 s
= 0.5 * (__imag__ x
/ r
);
133 s
= __ieee754_sqrt (0.5 * (d
- __real__ x
));
134 if (scale
== 1 && fabs (__imag__ x
) < 1.0)
136 /* Avoid possible intermediate underflow. */
137 r
= fabs (__imag__ x
/ s
);
138 s
= __scalbn (s
, scale
);
142 r
= fabs (0.5 * (__imag__ x
/ s
));
147 r
= __scalbn (r
, scale
);
148 s
= __scalbn (s
, scale
);
151 math_check_force_underflow (r
);
152 math_check_force_underflow (s
);
155 __imag__ res
= __copysign (s
, __imag__ x
);
161 weak_alias (__csqrt
, csqrt
)
162 #ifdef NO_LONG_DOUBLE
163 strong_alias (__csqrt
, __csqrtl
)
164 weak_alias (__csqrt
, csqrtl
)