1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __csqrt = csqrt
32
33
34 /*
35 * dcomplex csqrt(dcomplex z);
36 *
37 * 2 2 2
38 * Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y.
39 *
40 * Hence x = r*r-s*s, y = 2sr.
41 *
42 * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have
43 * ________
44 * 2 2 / 2 2
45 * (1) r + s = \/ x + y ,
46 *
47 * 2 2
48 * (2) r - s = x
49 *
50 * (3) 2sr = y.
51 *
52 * Perform (1)-(2) and (1)+(2), we obtain
53 *
54 * 2
55 * (4) 2 r = hypot(x,y)+x,
56 *
57 * 2
58 * (5) 2*s = hypot(x,y)-x
59 * ________
60 * / 2 2
61 * where hypot(x,y) = \/ x + y .
62 *
63 * In order to avoid numerical cancellation, we use formula (4) for
64 * positive x, and (5) for negative x. The other component is then
65 * computed by formula (3).
66 *
67 *
68 * ALGORITHM
69 * ------------------
70 *
71 * (assume x and y are of medium size, i.e., no over/underflow in squaring)
72 *
73 * If x >=0 then
74 * ________
75 * / 2 2
76 * 2 \/ x + y + x y
77 * r = ---------------------, s = -------; (6)
78 * 2 2 r
79 *
80 * (note that we choose sign(s) = sign(y) to force r >=0).
81 * Otherwise,
82 * ________
83 * / 2 2
84 * 2 \/ x + y - x y
85 * s = ---------------------, r = -------; (7)
86 * 2 2 s
87 *
88 * EXCEPTION:
89 *
90 * One may use the polar coordinate of a complex number to justify the
91 * following exception cases:
92 *
93 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
94 * csqrt(+-0+ i 0 ) = 0 + i 0
95 * csqrt( x + i inf ) = inf + i inf for all x (including NaN)
96 * csqrt( x + i NaN ) = NaN + i NaN with invalid for finite x
97 * csqrt(-inf+ iy ) = 0 + i inf for finite positive-signed y
98 * csqrt(+inf+ iy ) = inf + i 0 for finite positive-signed y
99 * csqrt(-inf+ i NaN) = NaN +-i inf
100 * csqrt(+inf+ i NaN) = inf + i NaN
101 * csqrt(NaN + i y ) = NaN + i NaN for finite y
102 * csqrt(NaN + i NaN) = NaN + i NaN
103 */
104
105 #include "libm.h" /* fabs/sqrt */
106 #include "complex_wrapper.h"
107
108 static const double two300 = 2.03703597633448608627e+90,
109 twom300 = 4.90909346529772655310e-91,
110 two599 = 2.07475778444049647926e+180,
111 twom601 = 1.20495993255144205887e-181,
112 two = 2.0,
113 zero = 0.0,
114 half = 0.5;
115
116
117 dcomplex
118 csqrt(dcomplex z)
119 {
120 dcomplex ans;
121 double x, y, t, ax, ay;
122 int n, ix, iy, hx, hy, lx, ly;
123
124 x = D_RE(z);
125 y = D_IM(z);
126 hx = HI_WORD(x);
127 lx = LO_WORD(x);
128 hy = HI_WORD(y);
129 ly = LO_WORD(y);
130 ix = hx & 0x7fffffff;
131 iy = hy & 0x7fffffff;
132 ay = fabs(y);
133 ax = fabs(x);
134
135 if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
136 /* x or y is Inf or NaN */
137 if (ISINF(iy, ly)) {
138 D_IM(ans) = D_RE(ans) = ay;
139 } else if (ISINF(ix, lx)) {
140 if (hx > 0) {
141 D_RE(ans) = ax;
142 D_IM(ans) = ay * zero;
143 } else {
144 D_RE(ans) = ay * zero;
145 D_IM(ans) = ax;
146 }
147 } else {
148 D_IM(ans) = D_RE(ans) = ax + ay;
149 }
150 } else if ((iy | ly) == 0) { /* y = 0 */
151 if (hx >= 0) {
152 D_RE(ans) = sqrt(ax);
153 D_IM(ans) = zero;
154 } else {
155 D_IM(ans) = sqrt(ax);
156 D_RE(ans) = zero;
157 }
158 } else if (ix >= iy) {
159 n = (ix - iy) >> 20;
160
161 if (n >= 30) { /* x >> y or y=0 */
162 t = sqrt(ax);
163 } else if (ix >= 0x5f300000) { /* x > 2**500 */
164 ax *= twom601;
165 y *= twom601;
166 t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
167 } else if (iy < 0x20b00000) { /* y < 2**-500 */
168 ax *= two599;
169 y *= two599;
170 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
171 } else {
172 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
173 }
174
175 if (hx >= 0) {
176 D_RE(ans) = t;
177 D_IM(ans) = ay / (t + t);
178 } else {
179 D_IM(ans) = t;
180 D_RE(ans) = ay / (t + t);
181 }
182 } else {
183 n = (iy - ix) >> 20;
184
185 if (n >= 30) { /* y >> x */
186 if (n >= 60)
187 t = sqrt(half * ay);
188 else if (iy >= 0x7fe00000)
189 t = sqrt(half * ay + half * ax);
190 else if (ix <= 0x00100000)
191 t = half * sqrt(two * (ay + ax));
192 else
193 t = sqrt(half * (ay + ax));
194 } else if (iy >= 0x5f300000) { /* y > 2**500 */
195 ax *= twom601;
196 y *= twom601;
197 t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
198 } else if (ix < 0x20b00000) { /* x < 2**-500 */
199 ax *= two599;
200 y *= two599;
201 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
202 } else {
203 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
204 }
205
206 if (hx >= 0) {
207 D_RE(ans) = t;
208 D_IM(ans) = ay / (t + t);
209 } else {
210 D_IM(ans) = t;
211 D_RE(ans) = ay / (t + t);
212 }
213 }
214
215 if (hy < 0)
216 D_IM(ans) = -D_IM(ans);
217
218 return (ans);
219 }