Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/R/powf.c
+++ new/usr/src/lib/libm/common/R/powf.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 __powf = powf
30 32
31 33 #include "libm.h"
32 -#include "xpg6.h" /* __xpg6 */
34 +#include "xpg6.h" /* __xpg6 */
33 35 #define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int
34 36
35 37 #if defined(__i386) && !defined(__amd64)
36 38 extern int __swapRP(int);
37 39 #endif
38 40
39 -/* INDENT OFF */
40 41 static const double
41 - ln2 = 6.93147180559945286227e-01, /* 0x3fe62e42, 0xfefa39ef */
42 - invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
42 + ln2 = 6.93147180559945286227e-01, /* 0x3fe62e42, 0xfefa39ef */
43 + invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
43 44 dtwo = 2.0,
44 45 done = 1.0,
45 46 dhalf = 0.5,
46 47 d32 = 32.0,
47 48 d1_32 = 0.03125,
48 49 A0 = 1.999999999813723303647511146995966439250e+0000,
49 50 A1 = 6.666910817935858533770138657139665608610e-0001,
50 51 t0 = 2.000000000004777489262405315073203746943e+0000,
51 52 t1 = 1.666663408349926379873111932994250726307e-0001;
52 53
53 54 static const double S[] = {
54 55 1.00000000000000000000e+00, /* 3FF0000000000000 */
55 56 1.02189714865411662714e+00, /* 3FF059B0D3158574 */
56 57 1.04427378242741375480e+00, /* 3FF0B5586CF9890F */
57 58 1.06714040067682369717e+00, /* 3FF11301D0125B51 */
58 59 1.09050773266525768967e+00, /* 3FF172B83C7D517B */
59 60 1.11438674259589243221e+00, /* 3FF1D4873168B9AA */
60 61 1.13878863475669156458e+00, /* 3FF2387A6E756238 */
61 62 1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */
62 63 1.18920711500272102690e+00, /* 3FF306FE0A31B715 */
63 64 1.21524735998046895524e+00, /* 3FF371A7373AA9CB */
64 65 1.24185781207348400201e+00, /* 3FF3DEA64C123422 */
65 66 1.26905095719173321989e+00, /* 3FF44E086061892D */
66 67 1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */
67 68 1.32523664315974132322e+00, /* 3FF5342B569D4F82 */
68 69 1.35425554693689265129e+00, /* 3FF5AB07DD485429 */
69 70 1.38390988196383202258e+00, /* 3FF6247EB03A5585 */
70 71 1.41421356237309514547e+00, /* 3FF6A09E667F3BCD */
71 72 1.44518080697704665027e+00, /* 3FF71F75E8EC5F74 */
72 73 1.47682614593949934623e+00, /* 3FF7A11473EB0187 */
73 74 1.50916442759342284141e+00, /* 3FF82589994CCE13 */
74 75 1.54221082540794074411e+00, /* 3FF8ACE5422AA0DB */
75 76 1.57598084510788649659e+00, /* 3FF93737B0CDC5E5 */
76 77 1.61049033194925428347e+00, /* 3FF9C49182A3F090 */
77 78 1.64575547815396494578e+00, /* 3FFA5503B23E255D */
78 79 1.68179283050742900407e+00, /* 3FFAE89F995AD3AD */
|
↓ open down ↓ |
26 lines elided |
↑ open up ↑ |
79 80 1.71861929812247793414e+00, /* 3FFB7F76F2FB5E47 */
80 81 1.75625216037329945351e+00, /* 3FFC199BDD85529C */
81 82 1.79470907500310716820e+00, /* 3FFCB720DCEF9069 */
82 83 1.83400808640934243066e+00, /* 3FFD5818DCFBA487 */
83 84 1.87416763411029996256e+00, /* 3FFDFC97337B9B5F */
84 85 1.91520656139714740007e+00, /* 3FFEA4AFA2A490DA */
85 86 1.95714412417540017941e+00, /* 3FFF50765B6E4540 */
86 87 };
87 88
88 89 static const double TBL[] = {
89 - 0.00000000000000000e+00,
90 - 3.07716586667536873e-02,
91 - 6.06246218164348399e-02,
92 - 8.96121586896871380e-02,
93 - 1.17783035656383456e-01,
94 - 1.45182009844497889e-01,
95 - 1.71850256926659228e-01,
96 - 1.97825743329919868e-01,
97 - 2.23143551314209765e-01,
98 - 2.47836163904581269e-01,
99 - 2.71933715483641758e-01,
100 - 2.95464212893835898e-01,
101 - 3.18453731118534589e-01,
102 - 3.40926586970593193e-01,
103 - 3.62905493689368475e-01,
104 - 3.84411698910332056e-01,
105 - 4.05465108108164385e-01,
106 - 4.26084395310900088e-01,
107 - 4.46287102628419530e-01,
108 - 4.66089729924599239e-01,
109 - 4.85507815781700824e-01,
110 - 5.04556010752395312e-01,
111 - 5.23248143764547868e-01,
112 - 5.41597282432744409e-01,
113 - 5.59615787935422659e-01,
114 - 5.77315365034823613e-01,
115 - 5.94707107746692776e-01,
116 - 6.11801541105992941e-01,
117 - 6.28608659422374094e-01,
118 - 6.45137961373584701e-01,
119 - 6.61398482245365016e-01,
120 - 6.77398823591806143e-01,
90 + 0.00000000000000000e+00, 3.07716586667536873e-02,
91 + 6.06246218164348399e-02, 8.96121586896871380e-02,
92 + 1.17783035656383456e-01, 1.45182009844497889e-01,
93 + 1.71850256926659228e-01, 1.97825743329919868e-01,
94 + 2.23143551314209765e-01, 2.47836163904581269e-01,
95 + 2.71933715483641758e-01, 2.95464212893835898e-01,
96 + 3.18453731118534589e-01, 3.40926586970593193e-01,
97 + 3.62905493689368475e-01, 3.84411698910332056e-01,
98 + 4.05465108108164385e-01, 4.26084395310900088e-01,
99 + 4.46287102628419530e-01, 4.66089729924599239e-01,
100 + 4.85507815781700824e-01, 5.04556010752395312e-01,
101 + 5.23248143764547868e-01, 5.41597282432744409e-01,
102 + 5.59615787935422659e-01, 5.77315365034823613e-01,
103 + 5.94707107746692776e-01, 6.11801541105992941e-01,
104 + 6.28608659422374094e-01, 6.45137961373584701e-01,
105 + 6.61398482245365016e-01, 6.77398823591806143e-01,
121 106 };
122 107
123 108 static const float zero = 0.0F, one = 1.0F, huge = 1.0e25f, tiny = 1.0e-25f;
124 -/* INDENT ON */
109 +
125 110
126 111 float
127 -powf(float x, float y) {
128 - float fx = x, fy = y;
129 - float fz;
130 - int ix, iy, jx, jy, k, iw, yisint;
112 +powf(float x, float y)
113 +{
114 + float fx = x, fy = y;
115 + float fz;
116 + int ix, iy, jx, jy, k, iw, yisint;
131 117
132 118 ix = *(int *)&x;
133 119 iy = *(int *)&y;
134 120 jx = ix & ~0x80000000;
135 121 jy = iy & ~0x80000000;
136 122
137 123 if (jy == 0)
138 - return (one); /* x**+-0 = 1 */
124 + return (one); /* x**+-0 = 1 */
139 125 else if (ix == 0x3f800000 && (__xpg6 & _C99SUSv3_pow) != 0)
140 - return (one); /* C99: 1**anything = 1 */
126 + return (one); /* C99: 1**anything = 1 */
141 127 else if (((0x7f800000 - jx) | (0x7f800000 - jy)) < 0)
142 128 return (fx * fy); /* at least one of x or y is NaN */
143 - /* includes Sun: 1**NaN = NaN */
144 - /* INDENT OFF */
129 +
130 + /*
131 + * includes Sun: 1**NaN = NaN
132 + */
133 +
145 134 /*
146 135 * determine if y is an odd int
147 136 * yisint = 0 ... y is not an integer
148 137 * yisint = 1 ... y is an odd int
149 138 * yisint = 2 ... y is an even int
150 139 */
151 - /* INDENT ON */
152 140 yisint = 0;
141 +
153 142 if (ix < 0) {
154 143 if (jy >= 0x4b800000) {
155 - yisint = 2; /* |y|>=2**24: y must be even */
144 + yisint = 2; /* |y|>=2**24: y must be even */
156 145 } else if (jy >= 0x3f800000) {
157 146 k = (jy >> 23) - 0x7f; /* exponent */
158 147 iw = jy >> (23 - k);
148 +
159 149 if ((iw << (23 - k)) == jy)
160 150 yisint = 2 - (iw & 1);
161 151 }
162 152 }
163 153
164 154 /* special value of y */
165 155 if ((jy & ~0x7f800000) == 0) {
166 - if (jy == 0x7f800000) { /* y is +-inf */
156 + if (jy == 0x7f800000) { /* y is +-inf */
167 157 if (jx == 0x3f800000) {
168 158 if ((__xpg6 & _C99SUSv3_pow) != 0)
169 159 fz = one;
170 - /* C99: (-1)**+-inf is 1 */
160 + /* C99: (-1)**+-inf is 1 */
171 161 else
172 162 fz = fy - fy;
173 - /* Sun: (+-1)**+-inf = NaN */
163 +
164 + /* Sun: (+-1)**+-inf = NaN */
174 165 } else if (jx > 0x3f800000) {
175 - /* (|x|>1)**+,-inf = inf,0 */
166 + /* (|x|>1)**+,-inf = inf,0 */
176 167 if (iy > 0)
177 168 fz = fy;
178 169 else
179 170 fz = zero;
180 - } else { /* (|x|<1)**-,+inf = inf,0 */
171 + } else { /* (|x|<1)**-,+inf = inf,0 */
181 172 if (iy < 0)
182 173 fz = -fy;
183 174 else
184 175 fz = zero;
185 176 }
177 +
186 178 return (fz);
187 179 } else if (jy == 0x3f800000) { /* y is +-1 */
188 180 if (iy < 0)
189 181 fx = one / fx; /* y is -1 */
182 +
190 183 return (fx);
191 184 } else if (iy == 0x40000000) { /* y is 2 */
192 185 return (fx * fx);
193 186 } else if (iy == 0x3f000000) { /* y is 0.5 */
194 187 if (jx != 0 && jx != 0x7f800000)
195 188 return (sqrtf(x));
196 189 }
197 190 }
198 191
199 192 /* special value of x */
200 193 if ((jx & ~0x7f800000) == 0) {
201 194 if (jx == 0x7f800000 || jx == 0 || jx == 0x3f800000) {
202 195 /* x is +-0,+-inf,-1; set fz = |x|**y */
203 196 *(int *)&fz = jx;
197 +
204 198 if (iy < 0)
205 199 fz = one / fz;
200 +
206 201 if (ix < 0) {
207 202 if (jx == 0x3f800000 && yisint == 0) {
208 203 /* (-1)**non-int is NaN */
209 204 fz = zero;
210 205 fz /= fz;
211 206 } else if (yisint == 1) {
212 207 /* (x<0)**odd = -(|x|**odd) */
213 208 fz = -fz;
214 209 }
215 210 }
211 +
216 212 return (fz);
217 213 }
218 214 }
219 215
220 216 /* (x<0)**(non-int) is NaN */
221 217 if (ix < 0 && yisint == 0) {
222 218 fz = zero;
223 219 return (fz / fz);
224 220 }
225 221
226 222 /*
227 223 * compute exp(y*log(|x|))
228 224 * fx = *(float *) &jx;
229 225 * fz = (float) exp(((double) fy) * log((double) fx));
230 226 */
231 227 {
232 - double dx, dy, dz, ds;
233 - int *px = (int *)&dx, *pz = (int *)&dz, i, n, m;
228 + double dx, dy, dz, ds;
229 + int *px = (int *)&dx, *pz = (int *)&dz, i, n, m;
230 +
234 231 #if defined(__i386) && !defined(__amd64)
235 - int rp = __swapRP(fp_extended);
232 + int rp = __swapRP(fp_extended);
236 233 #endif
237 234
238 235 fx = *(float *)&jx;
239 236 dx = (double)fx;
240 237
241 238 /* compute log(x)/ln2 */
242 239 i = px[HIWORD] + 0x4000;
243 240 n = (i >> 20) - 0x3ff;
244 241 pz[HIWORD] = i & 0xffff8000;
245 242 pz[LOWORD] = 0;
246 243 ds = (dx - dz) / (dx + dz);
247 244 i = (i >> 15) & 0x1f;
248 245 dz = ds * ds;
249 246 dy = invln2 * (TBL[i] + ds * (A0 + dz * A1));
247 +
250 248 if (n == 0)
251 249 dz = (double)fy * dy;
252 250 else
253 251 dz = (double)fy * (dy + (double)n);
254 252
255 253 /* compute exp2(dz=y*ln(x)) */
256 254 i = pz[HIWORD];
255 +
257 256 if ((i & ~0x80000000) >= 0x40640000) { /* |z| >= 160.0 */
258 - fz = (i > 0)? huge : tiny;
257 + fz = (i > 0) ? huge : tiny;
258 +
259 259 if (ix < 0 && yisint == 1)
260 260 fz *= -fz; /* (-ve)**(odd int) */
261 261 else
262 262 fz *= fz;
263 +
263 264 #if defined(__i386) && !defined(__amd64)
264 265 if (rp != fp_extended)
265 266 (void) __swapRP(rp);
266 267 #endif
267 268 return (fz);
268 269 }
269 270
270 271 n = (int)(d32 * dz + (i > 0 ? dhalf : -dhalf));
271 272 i = n & 0x1f;
272 273 m = n >> 5;
273 274 dy = ln2 * (dz - d1_32 * (double)n);
274 275 dx = S[i] * (done - (dtwo * dy) / (dy * (done - dy * t1) - t0));
276 +
275 277 if (m != 0)
276 278 px[HIWORD] += m << 20;
279 +
277 280 fz = (float)dx;
278 281 #if defined(__i386) && !defined(__amd64)
279 282 if (rp != fp_extended)
280 283 (void) __swapRP(rp);
281 284 #endif
282 285 }
283 286
284 287 /* end of computing exp(y*log(x)) */
285 288 if (ix < 0 && yisint == 1)
286 - fz = -fz; /* (-ve)**(odd int) */
289 + fz = -fz; /* (-ve)**(odd int) */
290 +
287 291 return (fz);
288 292 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX