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