Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/cpowl.c
+++ new/usr/src/lib/libm/common/complex/cpowl.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 __cpowl = cpowl
31 32
32 -#include "libm.h" /* __k_clog_rl/__k_atan2l */
33 +#include "libm.h" /* __k_clog_rl/__k_atan2l */
33 34 /* atan2l/atan2pil/exp2l/expl/fabsl/hypotl/isinfl/logl/powl/sincosl/sincospil */
34 35 #include "complex_wrapper.h"
35 36 #include "longdouble.h"
36 37
37 38 #if defined(__sparc)
38 -#define HALF(x) ((int *) &x)[3] = 0; ((int *) &x)[2] &= 0xfe000000
39 -#define LAST(x) ((int *) &x)[3]
39 +#define HALF(x) ((int *)&x)[3] = 0; ((int *)&x)[2] &= 0xfe000000
40 +#define LAST(x) ((int *)&x)[3]
40 41 #elif defined(__x86)
41 -#define HALF(x) ((int *) &x)[0] = 0
42 -#define LAST(x) ((int *) &x)[0]
42 +#define HALF(x) ((int *)&x)[0] = 0
43 +#define LAST(x) ((int *)&x)[0]
43 44 #endif
44 45
45 -/* INDENT OFF */
46 46 static const int hiinf = 0x7fff0000;
47 -static const long double
48 - tiny = 1.0e-4000L,
47 +static const long double tiny = 1.0e-4000L,
49 48 huge = 1.0e4000L,
50 49 #if defined(__x86)
51 - /* 43 significant bits, 21 trailing zeros */
50 +/* 43 significant bits, 21 trailing zeros */
52 51 ln2hil = 0.693147180559890330187045037746429443359375L,
53 52 ln2lol = 5.497923018708371174712471612513436025525412068e-14L,
54 -#else /* sparc */
55 - /* 0x3FF962E4 2FEFA39E F35793C7 00000000 */
53 +#else /* sparc */
54 +/* 0x3FF962E4 2FEFA39E F35793C7 00000000 */
56 55 ln2hil = 0.693147180559945309417231592858066493070671489074L,
57 56 ln2lol = 5.28600110075004828645286235820646730106802446566153e-25L,
58 57 #endif
59 - invln2 = 1.442695040888963407359924681001892137427e+0000L,
58 + invln2 = 1.442695040888963407359924681001892137427e+0000L,
60 59 one = 1.0L,
61 60 zero = 0.0L;
62 -/* INDENT ON */
63 61
64 62 /*
65 63 * Assuming |t[0]| > |t[1]| and |t[2]| > |t[3]|, sum4fpl subroutine
66 64 * compute t[0] + t[1] + t[2] + t[3] into two long double fp numbers.
67 65 */
68 -static long double sum4fpl(long double ta[], long double *w)
66 +static long double
67 +sum4fpl(long double ta[], long double *w)
69 68 {
70 69 long double t1, t2, t3, t4, w1, w2, t;
71 - t1 = ta[0]; t2 = ta[1]; t3 = ta[2]; t4 = ta[3];
70 +
71 + t1 = ta[0];
72 + t2 = ta[1];
73 + t3 = ta[2];
74 + t4 = ta[3];
75 +
72 76 /*
73 77 * Rearrange ti so that |t1| >= |t2| >= |t3| >= |t4|
74 78 */
75 79 if (fabsl(t4) > fabsl(t1)) {
76 - t = t1; t1 = t3; t3 = t;
77 - t = t2; t2 = t4; t4 = t;
80 + t = t1;
81 + t1 = t3;
82 + t3 = t;
83 + t = t2;
84 + t2 = t4;
85 + t4 = t;
78 86 } else if (fabsl(t3) > fabsl(t1)) {
79 - t = t1; t1 = t3;
87 + t = t1;
88 + t1 = t3;
89 +
80 90 if (fabsl(t4) > fabsl(t2)) {
81 - t3 = t4; t4 = t2; t2 = t;
91 + t3 = t4;
92 + t4 = t2;
93 + t2 = t;
82 94 } else {
83 - t3 = t2; t2 = t;
95 + t3 = t2;
96 + t2 = t;
84 97 }
85 98 } else if (fabsl(t3) > fabsl(t2)) {
86 - t = t2; t2 = t3;
99 + t = t2;
100 + t2 = t3;
101 +
87 102 if (fabsl(t4) > fabsl(t2)) {
88 - t3 = t4; t4 = t;
89 - } else
103 + t3 = t4;
104 + t4 = t;
105 + } else {
90 106 t3 = t;
107 + }
91 108 }
109 +
92 110 /* summing r = t1 + t2 + t3 + t4 to w1 + w2 */
93 111 w1 = t3 + t4;
94 112 w2 = t4 - (w1 - t3);
95 - t = t2 + w1;
113 + t = t2 + w1;
96 114 w2 += w1 - (t - t2);
97 115 w1 = t + w2;
98 116 w2 += t - w1;
99 - t = t1 + w1;
117 + t = t1 + w1;
100 118 w2 += w1 - (t - t1);
101 119 w1 = t + w2;
102 120 *w = w2 - (w1 - t);
103 121 return (w1);
104 122 }
105 123
106 124 ldcomplex
107 -cpowl(ldcomplex z, ldcomplex w) {
125 +cpowl(ldcomplex z, ldcomplex w)
126 +{
108 127 ldcomplex ans;
109 128 long double x, y, u, v, t, c, s, r;
110 129 long double t1, t2, t3, t4, x1, x2, y1, y2, u1, v1, b[4], w1, w2;
111 130 int ix, iy, hx, hy, hv, hu, iu, iv, i, j, k;
112 131
113 132 x = LD_RE(z);
114 133 y = LD_IM(z);
115 134 u = LD_RE(w);
116 135 v = LD_IM(w);
117 136 hx = HI_XWORD(x);
118 137 hy = HI_XWORD(y);
119 138 hu = HI_XWORD(u);
120 139 hv = HI_XWORD(v);
121 140 ix = hx & 0x7fffffff;
122 141 iy = hy & 0x7fffffff;
123 142 iu = hu & 0x7fffffff;
124 143 iv = hv & 0x7fffffff;
125 144
126 145 j = 0;
127 - if (v == zero) { /* z**(real) */
128 - if (u == one) { /* (anything) ** 1 is itself */
146 +
147 + if (v == zero) { /* z**(real) */
148 + if (u == one) { /* (anything) ** 1 is itself */
129 149 LD_RE(ans) = x;
130 150 LD_IM(ans) = y;
131 151 } else if (u == zero) { /* (anything) ** 0 is 1 */
132 152 LD_RE(ans) = one;
133 153 LD_IM(ans) = zero;
134 154 } else if (y == zero) { /* real ** real */
135 155 LD_IM(ans) = zero;
156 +
136 157 if (hx < 0 && ix < hiinf && iu < hiinf) {
137 - /* -x ** u is exp(i*pi*u)*pow(x,u) */
158 + /* -x ** u is exp(i*pi*u)*pow(x,u) */
138 159 r = powl(-x, u);
139 160 sincospil(u, &s, &c);
140 - LD_RE(ans) = (c == zero)? c: c * r;
141 - LD_IM(ans) = (s == zero)? s: s * r;
142 - } else
161 + LD_RE(ans) = (c == zero) ? c : c *r;
162 + LD_IM(ans) = (s == zero) ? s : s *r;
163 + } else {
143 164 LD_RE(ans) = powl(x, u);
165 + }
144 166 } else if (x == zero || ix >= hiinf || iy >= hiinf) {
145 - if (isnanl(x) || isnanl(y) || isnanl(u))
167 + if (isnanl(x) || isnanl(y) || isnanl(u)) {
146 168 LD_RE(ans) = LD_IM(ans) = x + y + u;
147 - else {
169 + } else {
148 170 if (x == zero)
149 171 r = fabsl(y);
150 172 else
151 173 r = fabsl(x) + fabsl(y);
174 +
152 175 t = atan2pil(y, x);
153 176 sincospil(t * u, &s, &c);
154 - LD_RE(ans) = (c == zero)? c: c * r;
155 - LD_IM(ans) = (s == zero)? s: s * r;
177 + LD_RE(ans) = (c == zero) ? c : c *r;
178 + LD_IM(ans) = (s == zero) ? s : s *r;
156 179 }
157 - } else if (fabsl(x) == fabsl(y)) { /* |x| = |y| */
180 + } else if (fabsl(x) == fabsl(y)) { /* |x| = |y| */
158 181 if (hx >= 0) {
159 - t = (hy >= 0)? 0.25L : -0.25L;
182 + t = (hy >= 0) ? 0.25L : -0.25L;
160 183 sincospil(t * u, &s, &c);
161 184 } else if ((LAST(u) & 3) == 0) {
162 - t = (hy >= 0)? 0.75L : -0.75L;
185 + t = (hy >= 0) ? 0.75L : -0.75L;
163 186 sincospil(t * u, &s, &c);
164 187 } else {
165 - r = (hy >= 0)? u : -u;
188 + r = (hy >= 0) ? u : -u;
166 189 t = -0.25L * r;
167 190 w1 = r + t;
168 191 w2 = t - (w1 - r);
169 192 sincospil(w1, &t1, &t2);
170 193 sincospil(w2, &t3, &t4);
171 194 s = t1 * t4 + t3 * t2;
172 195 c = t2 * t4 - t1 * t3;
173 196 }
197 +
174 198 if (ix < 0x3ffe0000) /* |x| < 1/2 */
175 199 r = powl(fabsl(x + x), u) * exp2l(-0.5L * u);
176 200 else if (ix >= 0x3fff0000 || iu < 0x400cfff8)
177 201 /* |x| >= 1 or |u| < 16383 */
178 202 r = powl(fabsl(x), u) * exp2l(0.5L * u);
179 - else /* special treatment */
203 + else /* special treatment */
180 204 j = 2;
205 +
181 206 if (j == 0) {
182 - LD_RE(ans) = (c == zero)? c: c * r;
183 - LD_IM(ans) = (s == zero)? s: s * r;
207 + LD_RE(ans) = (c == zero) ? c : c *r;
208 + LD_IM(ans) = (s == zero) ? s : s *r;
184 209 }
185 - } else
210 + } else {
186 211 j = 1;
212 + }
213 +
187 214 if (j == 0)
188 215 return (ans);
189 216 }
217 +
190 218 if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
191 219 /*
192 220 * non-zero imag part(s) with inf component(s) yields NaN
193 221 */
194 222 t = fabsl(x) + fabsl(y) + fabsl(u) + fabsl(v);
195 223 LD_RE(ans) = LD_IM(ans) = t - t;
196 224 } else {
197 - k = 0; /* no scaling */
225 + k = 0; /* no scaling */
226 +
198 227 if (iu > 0x7ffe0000 || iv > 0x7ffe0000) {
199 228 u *= 1.52587890625000000000e-05L;
200 229 v *= 1.52587890625000000000e-05L;
201 - k = 1; /* scale u and v by 2**-16 */
230 + k = 1; /* scale u and v by 2**-16 */
202 231 }
232 +
203 233 /*
204 234 * Use similated higher precision arithmetic to compute:
205 235 * r = u * log(hypot(x, y)) - v * atan2(y, x)
206 236 * q = u * atan2(y, x) + v * log(hypot(x, y))
207 237 */
208 238
209 239 t1 = __k_clog_rl(x, y, &t2);
210 240 t3 = __k_atan2l(y, x, &t4);
211 - x1 = t1; HALF(x1);
212 - y1 = t3; HALF(y1);
213 - u1 = u; HALF(u1);
214 - v1 = v; HALF(v1);
215 - x2 = t2 - (x1 - t1); /* log(hypot(x,y)) = x1 + x2 */
216 - y2 = t4 - (y1 - t3); /* atan2(y,x) = y1 + y2 */
241 + x1 = t1;
242 + HALF(x1);
243 + y1 = t3;
244 + HALF(y1);
245 + u1 = u;
246 + HALF(u1);
247 + v1 = v;
248 + HALF(v1);
249 + x2 = t2 - (x1 - t1); /* log(hypot(x,y)) = x1 + x2 */
250 + y2 = t4 - (y1 - t3); /* atan2(y,x) = y1 + y2 */
251 +
217 252 /* compute q = u * atan2(y, x) + v * log(hypot(x, y)) */
218 253 if (j != 2) {
219 254 b[0] = u1 * y1;
220 255 b[1] = (u - u1) * y1 + u * y2;
256 +
221 257 if (j == 1) { /* v = 0 */
222 258 w1 = b[0] + b[1];
223 259 w2 = b[1] - (w1 - b[0]);
224 260 } else {
225 261 b[2] = v1 * x1;
226 262 b[3] = (v - v1) * x1 + v * x2;
227 263 w1 = sum4fpl(b, &w2);
228 264 }
265 +
229 266 sincosl(w1, &t1, &t2);
230 267 sincosl(w2, &t3, &t4);
231 268 s = t1 * t4 + t3 * t2;
232 269 c = t2 * t4 - t1 * t3;
233 - if (k == 1) /* square j times */
270 +
271 + if (k == 1) { /* square j times */
234 272 for (i = 0; i < 10; i++) {
235 273 t1 = s * c;
236 274 c = (c + s) * (c - s);
237 275 s = t1 + t1;
238 276 }
277 + }
239 278 }
279 +
240 280 /* compute r = u * (t1, t2) - v * (t3, t4) */
241 281 b[0] = u1 * x1;
242 282 b[1] = (u - u1) * x1 + u * x2;
243 - if (j == 1) { /* v = 0 */
283 +
284 + if (j == 1) { /* v = 0 */
244 285 w1 = b[0] + b[1];
245 286 w2 = b[1] - (w1 - b[0]);
246 287 } else {
247 288 b[2] = -v1 * y1;
248 289 b[3] = (v1 - v) * y1 - v * y2;
249 290 w1 = sum4fpl(b, &w2);
250 291 }
292 +
251 293 /* scale back unless w1 is large enough to cause exception */
252 294 if (k != 0 && fabsl(w1) < 20000.0L) {
253 - w1 *= 65536.0L; w2 *= 65536.0L;
295 + w1 *= 65536.0L;
296 + w2 *= 65536.0L;
254 297 }
298 +
255 299 hx = HI_XWORD(w1);
256 300 ix = hx & 0x7fffffff;
257 301 /* compute exp(w1 + w2) */
258 302 k = 0;
259 - if (ix < 0x3f8c0000) /* exp(tiny < 2**-115) = 1 */
303 +
304 + if (ix < 0x3f8c0000) { /* exp(tiny < 2**-115) = 1 */
260 305 r = one;
261 - else if (ix >= 0x400c6760) /* overflow/underflow */
262 - r = (hx < 0)? tiny * tiny : huge * huge;
263 - else { /* compute exp(w1 + w2) */
264 - k = (int) (invln2 * w1 + ((hx >= 0)? 0.5L : -0.5L));
265 - t1 = (long double) k;
306 + } else if (ix >= 0x400c6760) { /* overflow/underflow */
307 + r = (hx < 0) ? tiny * tiny : huge * huge;
308 + } else { /* compute exp(w1 + w2) */
309 + k = (int)(invln2 * w1 + ((hx >= 0) ? 0.5L : -0.5L));
310 + t1 = (long double)k;
266 311 t2 = w1 - t1 * ln2hil;
267 312 t3 = w2 - t1 * ln2lol;
268 313 r = expl(t2 + t3);
269 314 }
270 - if (c != zero) c *= r;
271 - if (s != zero) s *= r;
315 +
316 + if (c != zero)
317 + c *= r;
318 +
319 + if (s != zero)
320 + s *= r;
321 +
272 322 if (k != 0) {
273 323 c = scalbnl(c, k);
274 324 s = scalbnl(s, k);
275 325 }
326 +
276 327 LD_RE(ans) = c;
277 328 LD_IM(ans) = s;
278 329 }
330 +
279 331 return (ans);
280 332 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX