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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak csqrt = __csqrt
31
32 /* INDENT OFF */
33 /*
34 * dcomplex csqrt(dcomplex z);
35 *
36 * 2 2 2
37 * Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y.
38 *
39 * Hence x = r*r-s*s, y = 2sr.
40 *
41 * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have
42 * ________
43 * 2 2 / 2 2
44 * (1) r + s = \/ x + y ,
45 *
46 * 2 2
47 * (2) r - s = x
48 *
49 * (3) 2sr = y.
50 *
51 * Perform (1)-(2) and (1)+(2), we obtain
52 *
53 * 2
54 * (4) 2 r = hypot(x,y)+x,
55 *
56 * 2
57 * (5) 2*s = hypot(x,y)-x
58 * ________
59 * / 2 2
60 * where hypot(x,y) = \/ x + y .
61 *
62 * In order to avoid numerical cancellation, we use formula (4) for
63 * positive x, and (5) for negative x. The other component is then
64 * computed by formula (3).
65 *
66 *
67 * ALGORITHM
68 * ------------------
69 *
70 * (assume x and y are of medium size, i.e., no over/underflow in squaring)
71 *
72 * If x >=0 then
73 * ________
74 * / 2 2
75 * 2 \/ x + y + x y
76 * r = ---------------------, s = -------; (6)
77 * 2 2 r
78 *
79 * (note that we choose sign(s) = sign(y) to force r >=0).
80 * Otherwise,
81 * ________
82 * / 2 2
83 * 2 \/ x + y - x y
84 * s = ---------------------, r = -------; (7)
85 * 2 2 s
86 *
87 * EXCEPTION:
88 *
89 * One may use the polar coordinate of a complex number to justify the
90 * following exception cases:
91 *
92 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
93 * csqrt(+-0+ i 0 ) = 0 + i 0
94 * csqrt( x + i inf ) = inf + i inf for all x (including NaN)
95 * csqrt( x + i NaN ) = NaN + i NaN with invalid for finite x
96 * csqrt(-inf+ iy ) = 0 + i inf for finite positive-signed y
97 * csqrt(+inf+ iy ) = inf + i 0 for finite positive-signed y
98 * csqrt(-inf+ i NaN) = NaN +-i inf
99 * csqrt(+inf+ i NaN) = inf + i NaN
100 * csqrt(NaN + i y ) = NaN + i NaN for finite y
101 * csqrt(NaN + i NaN) = NaN + i NaN
102 */
103 /* INDENT ON */
104
105 #include "libm.h" /* fabs/sqrt */
106 #include "complex_wrapper.h"
107
108 /* INDENT OFF */
109 static const double
110 two300 = 2.03703597633448608627e+90,
111 twom300 = 4.90909346529772655310e-91,
112 two599 = 2.07475778444049647926e+180,
113 twom601 = 1.20495993255144205887e-181,
114 two = 2.0,
115 zero = 0.0,
116 half = 0.5;
117 /* INDENT ON */
118
119 dcomplex
120 csqrt(dcomplex z) {
121 dcomplex ans;
122 double x, y, t, ax, ay;
123 int n, ix, iy, hx, hy, lx, ly;
124
125 x = D_RE(z);
126 y = D_IM(z);
127 hx = HI_WORD(x);
128 lx = LO_WORD(x);
129 hy = HI_WORD(y);
130 ly = LO_WORD(y);
131 ix = hx & 0x7fffffff;
132 iy = hy & 0x7fffffff;
133 ay = fabs(y);
134 ax = fabs(x);
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 } else if ((iy | ly) == 0) { /* y = 0 */
150 if (hx >= 0) {
151 D_RE(ans) = sqrt(ax);
152 D_IM(ans) = zero;
153 } else {
154 D_IM(ans) = sqrt(ax);
155 D_RE(ans) = zero;
156 }
157 } else if (ix >= iy) {
158 n = (ix - iy) >> 20;
159 if (n >= 30) { /* x >> y or y=0 */
160 t = sqrt(ax);
161 } else if (ix >= 0x5f300000) { /* x > 2**500 */
162 ax *= twom601;
163 y *= twom601;
164 t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
165 } else if (iy < 0x20b00000) { /* y < 2**-500 */
166 ax *= two599;
167 y *= two599;
168 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
169 } else
170 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
171 if (hx >= 0) {
172 D_RE(ans) = t;
173 D_IM(ans) = ay / (t + t);
174 } else {
175 D_IM(ans) = t;
176 D_RE(ans) = ay / (t + t);
177 }
178 } else {
179 n = (iy - ix) >> 20;
180 if (n >= 30) { /* y >> x */
181 if (n >= 60)
182 t = sqrt(half * ay);
183 else if (iy >= 0x7fe00000)
184 t = sqrt(half * ay + half * ax);
185 else if (ix <= 0x00100000)
186 t = half * sqrt(two * (ay + ax));
187 else
188 t = sqrt(half * (ay + ax));
189 } else if (iy >= 0x5f300000) { /* y > 2**500 */
190 ax *= twom601;
191 y *= twom601;
192 t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
193 } else if (ix < 0x20b00000) { /* x < 2**-500 */
194 ax *= two599;
195 y *= two599;
196 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
197 } else
198 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
199 if (hx >= 0) {
200 D_RE(ans) = t;
201 D_IM(ans) = ay / (t + t);
202 } else {
203 D_IM(ans) = t;
204 D_RE(ans) = ay / (t + t);
205 }
206 }
207 if (hy < 0)
208 D_IM(ans) = -D_IM(ans);
209 return (ans);
210 }