Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/cpow.c
+++ new/usr/src/lib/libm/common/complex/cpow.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 __cpow = cpow
31 32
32 -/* INDENT OFF */
33 +
33 34 /*
34 35 * dcomplex cpow(dcomplex z);
35 36 *
36 37 * z**w analytically equivalent to
37 38 *
38 39 * cpow(z,w) = cexp(w clog(z))
39 40 *
40 41 * Let z = x+iy, w = u+iv.
41 42 * Since
42 43 * _________
43 44 * / 2 2 -1 y
44 45 * log(x+iy) = log(\/ x + y ) + i tan (---)
45 46 * x
46 47 *
47 48 * 1 2 2 -1 y
48 49 * = --- log(x + y ) + i tan (---)
49 50 * 2 x
50 51 * u 2 2 -1 y
51 52 * (u+iv)* log(x+iy) = --- log(x + y ) - v tan (---) + (1)
52 53 * 2 x
53 54 *
54 55 * v 2 2 -1 y
55 56 * i * [ --- log(x + y ) + u tan (---) ] (2)
56 57 * 2 x
57 58 *
58 59 * = r + i q
59 60 *
60 61 * Therefore,
61 62 * w r+iq r
62 63 * z = e = e (cos(q)+i*sin(q))
63 64 * _______
64 65 * / 2 2
65 66 * r \/ x + y -v*atan2(y,x)
66 67 * Here e can be expressed as: u * e
67 68 *
68 69 * Special cases (in the order of appearance):
69 70 * 1. (anything) ** 0 is 1
70 71 * 2. (anything) ** 1 is itself
71 72 * 3. When v = 0, y = 0:
72 73 * If x is finite and negative, and u is finite, then
73 74 * x ** u = exp(u*pi i) * pow(|x|, u);
74 75 * otherwise,
75 76 * x ** u = pow(x, u);
76 77 * 4. When v = 0, x = 0 or |x| = |y| or x is inf or y is inf:
77 78 * (x + y i) ** u = r * exp(q i)
78 79 * where
|
↓ open down ↓ |
36 lines elided |
↑ open up ↑ |
79 80 * r = hypot(x,y) ** u
80 81 * q = u * atan2pi(y, x)
81 82 *
82 83 * 5. otherwise, z**w is NAN if any x, y, u, v is a Nan or inf
83 84 *
84 85 * Note: many results of special cases are obtained in terms of
85 86 * polar coordinate. In the conversion from polar to rectangle:
86 87 * r exp(q i) = r * cos(q) + r * sin(q) i,
87 88 * we regard r * 0 is 0 except when r is a NaN.
88 89 */
89 -/* INDENT ON */
90 90
91 -#include "libm.h" /* atan2/exp/fabs/hypot/log/pow/scalbn */
92 - /* atan2pi/exp2/sincos/sincospi/__k_clog_r/__k_atan2 */
91 +/*
92 + * atan2/exp/fabs/hypot/log/pow/scalbn
93 + * atan2pi/exp2/sincos/sincospi/__k_clog_r/__k_atan2
94 + */
95 +#include "libm.h"
93 96 #include "complex_wrapper.h"
94 97
95 98 extern void sincospi(double, double *, double *);
96 -
97 -static const double
98 - huge = 1e300,
99 +static const double huge = 1e300,
99 100 tiny = 1e-300,
100 101 invln2 = 1.44269504088896338700e+00,
101 - ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
102 - ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
102 + ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
103 + ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
103 104 one = 1.0,
104 105 zero = 0.0;
105 -
106 106 static const int hiinf = 0x7ff00000;
107 107 extern double atan2pi(double, double);
108 108
109 109 /*
110 110 * Assuming |t[0]| > |t[1]| and |t[2]| > |t[3]|, sum4fp subroutine
111 111 * compute t[0] + t[1] + t[2] + t[3] into two double fp numbers.
112 112 */
113 113 static double
114 -sum4fp(double ta[], double *w) {
114 +sum4fp(double ta[], double *w)
115 +{
115 116 double t1, t2, t3, t4, w1, w2, t;
116 - t1 = ta[0]; t2 = ta[1]; t3 = ta[2]; t4 = ta[3];
117 +
118 + t1 = ta[0];
119 + t2 = ta[1];
120 + t3 = ta[2];
121 + t4 = ta[3];
122 +
117 123 /*
118 124 * Rearrange ti so that |t1| >= |t2| >= |t3| >= |t4|
119 125 */
120 126 if (fabs(t4) > fabs(t1)) {
121 - t = t1; t1 = t3; t3 = t;
122 - t = t2; t2 = t4; t4 = t;
127 + t = t1;
128 + t1 = t3;
129 + t3 = t;
130 + t = t2;
131 + t2 = t4;
132 + t4 = t;
123 133 } else if (fabs(t3) > fabs(t1)) {
124 - t = t1; t1 = t3;
134 + t = t1;
135 + t1 = t3;
136 +
125 137 if (fabs(t4) > fabs(t2)) {
126 - t3 = t4; t4 = t2; t2 = t;
138 + t3 = t4;
139 + t4 = t2;
140 + t2 = t;
127 141 } else {
128 - t3 = t2; t2 = t;
142 + t3 = t2;
143 + t2 = t;
129 144 }
130 145 } else if (fabs(t3) > fabs(t2)) {
131 - t = t2; t2 = t3;
146 + t = t2;
147 + t2 = t3;
148 +
132 149 if (fabs(t4) > fabs(t2)) {
133 - t3 = t4; t4 = t;
134 - } else
150 + t3 = t4;
151 + t4 = t;
152 + } else {
135 153 t3 = t;
154 + }
136 155 }
156 +
137 157 /* summing r = t1 + t2 + t3 + t4 to w1 + w2 */
138 158 w1 = t3 + t4;
139 159 w2 = t4 - (w1 - t3);
140 - t = t2 + w1;
160 + t = t2 + w1;
141 161 w2 += w1 - (t - t2);
142 162 w1 = t + w2;
143 163 w2 += t - w1;
144 - t = t1 + w1;
164 + t = t1 + w1;
145 165 w2 += w1 - (t - t1);
146 166 w1 = t + w2;
147 167 *w = w2 - (w1 - t);
148 168 return (w1);
149 169 }
150 170
151 171 dcomplex
152 -cpow(dcomplex z, dcomplex w) {
172 +cpow(dcomplex z, dcomplex w)
173 +{
153 174 dcomplex ans;
154 175 double x, y, u, v, t, c, s, r, x2, y2;
155 176 double b[4], t1, t2, t3, t4, w1, w2, u1, v1, x1, y1;
156 177 int ix, iy, hx, lx, hy, ly, hv, hu, iu, iv, lu, lv;
157 178 int i, j, k;
158 179
159 180 x = D_RE(z);
160 181 y = D_IM(z);
161 182 u = D_RE(w);
162 183 v = D_IM(w);
163 - hx = ((int *) &x)[HIWORD];
164 - lx = ((int *) &x)[LOWORD];
165 - hy = ((int *) &y)[HIWORD];
166 - ly = ((int *) &y)[LOWORD];
167 - hu = ((int *) &u)[HIWORD];
168 - lu = ((int *) &u)[LOWORD];
169 - hv = ((int *) &v)[HIWORD];
170 - lv = ((int *) &v)[LOWORD];
184 + hx = ((int *)&x)[HIWORD];
185 + lx = ((int *)&x)[LOWORD];
186 + hy = ((int *)&y)[HIWORD];
187 + ly = ((int *)&y)[LOWORD];
188 + hu = ((int *)&u)[HIWORD];
189 + lu = ((int *)&u)[LOWORD];
190 + hv = ((int *)&v)[HIWORD];
191 + lv = ((int *)&v)[LOWORD];
171 192 ix = hx & 0x7fffffff;
172 193 iy = hy & 0x7fffffff;
173 194 iu = hu & 0x7fffffff;
174 195 iv = hv & 0x7fffffff;
175 196
176 197 j = 0;
177 - if ((iv | lv) == 0) { /* z**(real) */
198 +
199 + if ((iv | lv) == 0) { /* z**(real) */
178 200 if (((hu - 0x3ff00000) | lu) == 0) { /* z ** 1 = z */
179 201 D_RE(ans) = x;
180 202 D_IM(ans) = y;
181 - } else if ((iu | lu) == 0) { /* z ** 0 = 1 */
203 + } else if ((iu | lu) == 0) { /* z ** 0 = 1 */
182 204 D_RE(ans) = one;
183 205 D_IM(ans) = zero;
184 - } else if ((iy | ly) == 0) { /* (real)**(real) */
206 + } else if ((iy | ly) == 0) { /* (real)**(real) */
185 207 D_IM(ans) = zero;
208 +
186 209 if (hx < 0 && ix < hiinf && iu < hiinf) {
187 210 /* -x ** u is exp(i*pi*u)*pow(x,u) */
188 211 r = pow(-x, u);
189 212 sincospi(u, &s, &c);
190 - D_RE(ans) = (c == zero)? c: c * r;
191 - D_IM(ans) = (s == zero)? s: s * r;
192 - } else
213 + D_RE(ans) = (c == zero) ? c : c *r;
214 + D_IM(ans) = (s == zero) ? s : s *r;
215 + } else {
193 216 D_RE(ans) = pow(x, u);
217 + }
194 218 } else if (((ix | lx) == 0) || ix >= hiinf || iy >= hiinf) {
195 - if (isnan(x) || isnan(y) || isnan(u))
219 + if (isnan(x) || isnan(y) || isnan(u)) {
196 220 D_RE(ans) = D_IM(ans) = x + y + u;
197 - else {
221 + } else {
198 222 if ((ix | lx) == 0)
199 223 r = fabs(y);
200 224 else
201 225 r = fabs(x) + fabs(y);
226 +
202 227 t = atan2pi(y, x);
203 228 sincospi(t * u, &s, &c);
204 - D_RE(ans) = (c == zero)? c: c * r;
205 - D_IM(ans) = (s == zero)? s: s * r;
229 + D_RE(ans) = (c == zero) ? c : c *r;
230 + D_IM(ans) = (s == zero) ? s : s *r;
206 231 }
207 - } else if (((ix - iy) | (lx - ly)) == 0) { /* |x| = |y| */
232 + } else if (((ix - iy) | (lx - ly)) == 0) { /* |x| = |y| */
208 233 if (hx >= 0) {
209 - t = (hy >= 0)? 0.25 : -0.25;
234 + t = (hy >= 0) ? 0.25 : -0.25;
210 235 sincospi(t * u, &s, &c);
211 236 } else if ((lu & 3) == 0) {
212 - t = (hy >= 0)? 0.75 : -0.75;
237 + t = (hy >= 0) ? 0.75 : -0.75;
213 238 sincospi(t * u, &s, &c);
214 239 } else {
215 - r = (hy >= 0)? u : -u;
240 + r = (hy >= 0) ? u : -u;
216 241 t = -0.25 * r;
217 242 w1 = r + t;
218 243 w2 = t - (w1 - r);
219 244 sincospi(w1, &t1, &t2);
220 245 sincospi(w2, &t3, &t4);
221 246 s = t1 * t4 + t3 * t2;
222 247 c = t2 * t4 - t1 * t3;
223 248 }
249 +
224 250 if (ix < 0x3fe00000) /* |x| < 1/2 */
225 251 r = pow(fabs(x + x), u) * exp2(-0.5 * u);
226 252 else if (ix >= 0x3ff00000 || iu < 0x408ff800)
227 253 /* |x| >= 1 or |u| < 1023 */
228 254 r = pow(fabs(x), u) * exp2(0.5 * u);
229 - else /* special treatment */
255 + else /* special treatment */
230 256 j = 2;
257 +
231 258 if (j == 0) {
232 - D_RE(ans) = (c == zero)? c: c * r;
233 - D_IM(ans) = (s == zero)? s: s * r;
259 + D_RE(ans) = (c == zero) ? c : c *r;
260 + D_IM(ans) = (s == zero) ? s : s *r;
234 261 }
235 - } else
262 + } else {
236 263 j = 1;
264 + }
265 +
237 266 if (j == 0)
238 267 return (ans);
239 268 }
269 +
240 270 if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
241 271 /*
242 272 * non-zero imag part(s) with inf component(s) yields NaN
243 273 */
244 274 t = fabs(x) + fabs(y) + fabs(u) + fabs(v);
245 275 D_RE(ans) = D_IM(ans) = t - t;
246 276 } else {
247 - k = 0; /* no scaling */
277 + k = 0; /* no scaling */
278 +
248 279 if (iu > 0x7f000000 || iv > 0x7f000000) {
249 280 u *= .0009765625; /* scale 2**-10 to avoid overflow */
250 281 v *= .0009765625;
251 - k = 1; /* scale by 2**-10 */
282 + k = 1; /* scale by 2**-10 */
252 283 }
284 +
253 285 /*
254 286 * Use similated higher precision arithmetic to compute:
255 287 * r = u * log(hypot(x, y)) - v * atan2(y, x)
256 288 * q = u * atan2(y, x) + v * log(hypot(x, y))
257 289 */
258 290 t1 = __k_clog_r(x, y, &t2);
259 291 t3 = __k_atan2(y, x, &t4);
260 292 x1 = t1;
261 293 y1 = t3;
262 294 u1 = u;
263 295 v1 = v;
264 - ((int *) &u1)[LOWORD] &= 0xf8000000;
265 - ((int *) &v1)[LOWORD] &= 0xf8000000;
266 - ((int *) &x1)[LOWORD] &= 0xf8000000;
267 - ((int *) &y1)[LOWORD] &= 0xf8000000;
296 + ((int *)&u1)[LOWORD] &= 0xf8000000;
297 + ((int *)&v1)[LOWORD] &= 0xf8000000;
298 + ((int *)&x1)[LOWORD] &= 0xf8000000;
299 + ((int *)&y1)[LOWORD] &= 0xf8000000;
268 300 x2 = t2 - (x1 - t1); /* log(hypot(x,y)) = x1 + x2 */
269 301 y2 = t4 - (y1 - t3); /* atan2(y,x) = y1 + y2 */
302 +
270 303 /* compute q = u * atan2(y, x) + v * log(hypot(x, y)) */
271 304 if (j != 2) {
272 305 b[0] = u1 * y1;
273 306 b[1] = (u - u1) * y1 + u * y2;
307 +
274 308 if (j == 1) { /* v = 0 */
275 309 w1 = b[0] + b[1];
276 310 w2 = b[1] - (w1 - b[0]);
277 311 } else {
278 312 b[2] = v1 * x1;
279 313 b[3] = (v - v1) * x1 + v * x2;
280 314 w1 = sum4fp(b, &w2);
281 315 }
316 +
282 317 sincos(w1, &t1, &t2);
283 318 sincos(w2, &t3, &t4);
284 319 s = t1 * t4 + t3 * t2;
285 320 c = t2 * t4 - t1 * t3;
286 - if (k == 1)
287 - /*
288 - * square (cos(q) + i sin(q)) k times to get
289 - * (cos(2^k * q + i sin(2^k * q)
290 - */
321 +
322 + if (k == 1) {
323 + /*
324 + * square (cos(q) + i sin(q)) k times to get
325 + * (cos(2^k * q + i sin(2^k * q)
326 + */
291 327 for (i = 0; i < 10; i++) {
292 328 t1 = s * c;
293 329 c = (c + s) * (c - s);
294 330 s = t1 + t1;
295 331 }
332 + }
296 333 }
334 +
297 335 /* compute r = u * (t1, t2) - v * (t3, t4) */
298 336 b[0] = u1 * x1;
299 337 b[1] = (u - u1) * x1 + u * x2;
300 - if (j == 1) { /* v = 0 */
338 +
339 + if (j == 1) { /* v = 0 */
301 340 w1 = b[0] + b[1];
302 341 w2 = b[1] - (w1 - b[0]);
303 342 } else {
304 343 b[2] = -v1 * y1;
305 344 b[3] = (v1 - v) * y1 - v * y2;
306 345 w1 = sum4fp(b, &w2);
307 346 }
347 +
308 348 /* check over/underflow for exp(w1 + w2) */
309 349 if (k && fabs(w1) < 1000.0) {
310 - w1 *= 1024; w2 *= 1024; k = 0;
350 + w1 *= 1024;
351 + w2 *= 1024;
352 + k = 0;
311 353 }
312 - hx = ((int *) &w1)[HIWORD];
313 - lx = ((int *) &w1)[LOWORD];
354 +
355 + hx = ((int *)&w1)[HIWORD];
356 + lx = ((int *)&w1)[LOWORD];
314 357 ix = hx & 0x7fffffff;
358 +
315 359 /* compute exp(w1 + w2) */
316 - if (ix < 0x3c900000) /* exp(tiny < 2**-54) = 1 */
360 + if (ix < 0x3c900000) { /* exp(tiny < 2**-54) = 1 */
317 361 r = one;
318 - else if (ix >= 0x40880000) /* overflow/underflow */
319 - r = (hx < 0)? tiny * tiny : huge * huge;
320 - else { /* compute exp(w1 + w2) */
321 - k = (int) (invln2 * w1 + ((hx >= 0)? 0.5 : -0.5));
322 - t1 = (double) k;
362 + } else if (ix >= 0x40880000) { /* overflow/underflow */
363 + r = (hx < 0) ? tiny * tiny : huge * huge;
364 + } else { /* compute exp(w1 + w2) */
365 + k = (int)(invln2 * w1 + ((hx >= 0) ? 0.5 : -0.5));
366 + t1 = (double)k;
323 367 t2 = w1 - t1 * ln2hi;
324 368 t3 = w2 - t1 * ln2lo;
325 369 r = exp(t2 + t3);
326 370 }
327 - if (c != zero) c *= r;
328 - if (s != zero) s *= r;
371 +
372 + if (c != zero)
373 + c *= r;
374 +
375 + if (s != zero)
376 + s *= r;
377 +
329 378 if (k != 0) {
330 379 c = scalbn(c, k);
331 380 s = scalbn(s, k);
332 381 }
382 +
333 383 D_RE(ans) = c;
334 384 D_IM(ans) = s;
335 385 }
386 +
336 387 return (ans);
337 388 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX