Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/cpowf.c
+++ new/usr/src/lib/libm/common/complex/cpowf.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
|
↓ open down ↓ |
10 lines elided |
↑ open up ↑ |
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
20 20 */
21 +
21 22 /*
22 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 24 */
25 +
24 26 /*
25 27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 28 * Use is subject to license terms.
27 29 */
28 30
29 31 #pragma weak __cpowf = cpowf
30 32
31 33 #include "libm.h"
32 34 #include "complex_wrapper.h"
33 35
34 36 extern void sincospi(double, double *, double *);
35 37 extern void sincospif(float, float *, float *);
36 38 extern double atan2pi(double, double);
37 39 extern float atan2pif(float, float);
38 40
39 41 #if defined(__i386) && !defined(__amd64)
40 42 extern int __swapRP(int);
41 43 #endif
42 44
43 -static const double
44 - dpi = 3.1415926535897931160E0, /* Hex 2^ 1 * 1.921FB54442D18 */
45 - dhalf = 0.5,
45 +/* Hex 2^ 1 * 1.921FB54442D18 */
46 +static const double dpi = 3.1415926535897931160E0;
47 +
48 +static const double dhalf = 0.5,
46 49 dsqrt2 = 1.41421356237309514547, /* 3FF6A09E 667F3BCD */
47 50 dinvpi = 0.3183098861837906715377675;
48 51
49 52 static const float one = 1.0F, zero = 0.0F;
50 53
51 54 #define hiinf 0x7f800000
52 55
53 56 fcomplex
54 -cpowf(fcomplex z, fcomplex w) {
55 - fcomplex ans;
56 - float x, y, u, v, t, c, s;
57 - double dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
58 - int ix, iy, hx, hy, hv, hu, iu, iv, j;
57 +cpowf(fcomplex z, fcomplex w)
58 +{
59 + fcomplex ans;
60 + float x, y, u, v, t, c, s;
61 + double dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
62 + int ix, iy, hx, hy, hv, hu, iu, iv, j;
59 63
60 64 x = F_RE(z);
61 65 y = F_IM(z);
62 66 u = F_RE(w);
63 67 v = F_IM(w);
64 68 hx = THE_WORD(x);
65 69 hy = THE_WORD(y);
66 70 hu = THE_WORD(u);
67 71 hv = THE_WORD(v);
68 72 ix = hx & 0x7fffffff;
69 73 iy = hy & 0x7fffffff;
70 74 iu = hu & 0x7fffffff;
71 75 iv = hv & 0x7fffffff;
72 76
73 77 j = 0;
74 - if (iv == 0) { /* z**(real) */
78 +
79 + if (iv == 0) { /* z**(real) */
75 80 if (hu == 0x3f800000) { /* (anything) ** 1 is itself */
76 81 F_RE(ans) = x;
77 82 F_IM(ans) = y;
78 83 } else if (iu == 0) { /* (anything) ** 0 is 1 */
79 84 F_RE(ans) = one;
80 85 F_IM(ans) = zero;
81 86 } else if (iy == 0) { /* (real)**(real) */
82 87 F_IM(ans) = zero;
88 +
83 89 if (hx < 0 && ix < hiinf && iu < hiinf) {
84 90 /* -x ** u is exp(i*pi*u)*pow(x,u) */
85 91 t = powf(-x, u);
86 92 sincospif(u, &s, &c);
87 - F_RE(ans) = (c == zero)? c: c * t;
88 - F_IM(ans) = (s == zero)? s: s * t;
93 + F_RE(ans) = (c == zero) ? c : c *t;
94 + F_IM(ans) = (s == zero) ? s : s *t;
89 95 } else {
90 96 F_RE(ans) = powf(x, u);
91 97 }
92 98 } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
93 99 if (ix > hiinf || iy > hiinf || iu > hiinf) {
94 100 F_RE(ans) = F_IM(ans) = x + y + u;
95 101 } else {
96 102 v = fabsf(y);
103 +
97 104 if (ix != 0)
98 105 v += fabsf(x);
106 +
99 107 t = atan2pif(y, x);
100 108 sincospif(t * u, &s, &c);
101 - F_RE(ans) = (c == zero)? c: c * v;
102 - F_IM(ans) = (s == zero)? s: s * v;
109 + F_RE(ans) = (c == zero) ? c : c *v;
110 + F_IM(ans) = (s == zero) ? s : s *v;
103 111 }
104 112 } else if (ix == iy) { /* if |x| == |y| */
105 113 #if defined(__i386) && !defined(__amd64)
106 - int rp = __swapRP(fp_extended);
114 + int rp = __swapRP(fp_extended);
107 115 #endif
108 116 dx = (double)x;
109 117 du = (double)u;
110 - dt = (hx >= 0)? 0.25 : 0.75;
118 + dt = (hx >= 0) ? 0.25 : 0.75;
119 +
111 120 if (hy < 0)
112 121 dt = -dt;
122 +
113 123 dr = pow(dsqrt2 * dx, du);
114 124 sincospi(dt * du, &ds, &dc);
115 125 F_RE(ans) = (float)(dr * dc);
116 126 F_IM(ans) = (float)(dr * ds);
117 127 #if defined(__i386) && !defined(__amd64)
118 128 if (rp != fp_extended)
119 129 (void) __swapRP(rp);
120 130 #endif
121 131 } else {
122 132 j = 1;
123 133 }
134 +
124 135 if (j == 0)
125 136 return (ans);
126 137 }
138 +
127 139 if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
128 140 /*
129 141 * non-zero imaginery part(s) with inf component(s) yields NaN
130 142 */
131 143 t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
132 144 F_RE(ans) = F_IM(ans) = t - t;
133 145 } else {
134 146 #if defined(__i386) && !defined(__amd64)
135 - int rp = __swapRP(fp_extended);
147 + int rp = __swapRP(fp_extended);
136 148 #endif
137 - /* INDENT OFF */
149 +
138 150 /*
139 151 * r = u*log(hypot(x,y))-v*atan2(y,x),
140 152 * q = u*atan2(y,x)+v*log(hypot(x,y))
141 153 * or
142 154 * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
143 155 * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
144 156 * ans = exp(r)*(cospi(q/pi) + i sinpi(q/pi))
145 157 */
146 - /* INDENT ON */
147 158 dx = (double)x;
148 159 dy = (double)y;
149 160 du = (double)u;
150 161 dv = (double)v;
162 +
151 163 if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */
152 164 dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
153 165 else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
154 166 dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);
155 167 else
156 168 dt = dhalf * log(dx * dx + dy * dy);
169 +
157 170 dp = atan2pi(dy, dx);
158 - if (iv == 0) { /* dv = 0 */
171 +
172 + if (iv == 0) { /* dv = 0 */
159 173 dr = exp(du * dt);
160 174 dq = du * dp;
161 175 } else {
162 176 dr = exp(du * dt - dv * dp * dpi);
163 177 dq = du * dp + dv * dt * dinvpi;
164 178 }
179 +
165 180 sincospi(dq, &ds, &dc);
166 181 F_RE(ans) = (float)(dr * dc);
167 182 F_IM(ans) = (float)(dr * ds);
168 183 #if defined(__i386) && !defined(__amd64)
169 184 if (rp != fp_extended)
170 185 (void) __swapRP(rp);
171 186 #endif
172 187 }
188 +
173 189 return (ans);
174 190 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX