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
csqrtq(__complex128 x)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 }
157