xref: /illumos-gate/usr/src/lib/libm/common/complex/csqrt.c (revision ddc0e0b53c661f6e439e3b7072b3ef353eadb4af)
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
csqrt(dcomplex z)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 }
211