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