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 *
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 }
|
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 *
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 }
|