Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/pow.c
+++ new/usr/src/lib/libm/common/C/pow.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 __pow = pow
31 32
32 33 /*
33 34 * pow(x,y) return x**y
34 35 * n
35 36 * Method: Let x = 2 * (1+f)
36 37 * 1. Compute and return log2(x) in two pieces:
37 38 * log2(x) = w1 + w2,
38 39 * where w1 has 24 bits trailing zero.
39 40 * 2. Perform y*log2(x) by simulating muti-precision arithmetic
40 41 * 3. Return x**y = exp2(y*log(x))
41 42 *
42 43 * Special cases:
43 44 * 1. (anything) ** +-0 is 1
44 45 * 1'. 1 ** (anything) is 1 (C99; 1 ** +-INF/NAN used to be NAN)
45 46 * 2. (anything) ** 1 is itself
46 47 * 3. (anything except 1) ** NAN is NAN ("except 1" is C99)
47 48 * 4. NAN ** (anything except 0) is NAN
48 49 * 5. +-(|x| > 1) ** +INF is +INF
49 50 * 6. +-(|x| > 1) ** -INF is +0
50 51 * 7. +-(|x| < 1) ** +INF is +0
51 52 * 8. +-(|x| < 1) ** -INF is +INF
52 53 * 9. -1 ** +-INF is 1 (C99; -1 ** +-INF used to be NAN)
53 54 * 10. +0 ** (+anything except 0, NAN) is +0
54 55 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
55 56 * 12. +0 ** (-anything except 0, NAN) is +INF
56 57 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
57 58 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
58 59 * 15. +INF ** (+anything except 0,NAN) is +INF
59 60 * 16. +INF ** (-anything except 0,NAN) is +0
60 61 * 17. -INF ** (anything) = -0 ** (-anything)
|
↓ open down ↓ |
26 lines elided |
↑ open up ↑ |
61 62 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
62 63 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
63 64 *
64 65 * Accuracy:
65 66 * pow(x,y) returns x**y nearly rounded. In particular
66 67 * pow(integer,integer)
67 68 * always returns the correct integer provided it is representable.
68 69 */
69 70
70 71 #include "libm.h"
71 -#include "xpg6.h" /* __xpg6 */
72 +#include "xpg6.h" /* __xpg6 */
72 73 #define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int
73 74
74 -static const double zero = 0.0, one = 1.0, two = 2.0;
75 +static const double zero = 0.0,
76 + one = 1.0,
77 + two = 2.0;
75 78
76 79 extern const double _TBL_log2_hi[], _TBL_log2_lo[];
77 -static const double
78 - two53 = 9007199254740992.0,
80 +
81 +static const double two53 = 9007199254740992.0,
79 82 A1_hi = 2.8853900432586669921875,
80 83 A1_lo = 3.8519259825035041963606002e-8,
81 84 A1 = 2.885390081777926817222541963606002026086e+0000,
82 85 A2 = 9.617966939207270828380543979852286255862e-0001,
83 86 A3 = 5.770807680887875964868853124873696201995e-0001,
84 87 B0_hi = 2.8853900432586669921875,
85 88 B0_lo = 3.8519259822532793056374320585e-8,
86 89 B0 = 2.885390081777926814720293056374320585689e+0000,
87 90 B1 = 9.617966939259755138949202350396200257632e-0001,
88 91 B2 = 5.770780163585687000782112776448797953382e-0001,
89 92 B3 = 4.121985488948771523290174512461778354953e-0001,
90 93 B4 = 3.207590534812432970433641789022666850193e-0001;
91 94
92 95 static double
93 -log2_x(double x, double *w) {
96 +log2_x(double x, double *w)
97 +{
94 98 double f, s, z, qn, h, t;
95 - int *px = (int *) &x;
96 - int *pz = (int *) &z;
99 + int *px = (int *)&x;
100 + int *pz = (int *)&z;
97 101 int i, j, ix, n;
98 102
99 103 n = 0;
100 104 ix = px[HIWORD];
105 +
101 106 if (ix >= 0x3fef03f1 && ix < 0x3ff08208) { /* 65/63 > x > 63/65 */
102 107 double f1, v;
108 +
103 109 f = x - one;
110 +
104 111 if (((ix - 0x3ff00000) | px[LOWORD]) == 0) {
105 112 *w = zero;
106 - return (zero); /* log2(1)= +0 */
113 + return (zero); /* log2(1)= +0 */
107 114 }
115 +
108 116 qn = one / (two + f);
109 - s = f * qn; /* |s|<2**-6 */
117 + s = f * qn; /* |s|<2**-6 */
110 118 v = s * s;
111 - h = (double) ((float) s);
112 - f1 = (double) ((float) f);
119 + h = (double)((float)s);
120 + f1 = (double)((float)f);
113 121 t = qn * (((f - two * h) - h * f1) - h * (f - f1));
114 - /* s = h+t */
122 + /* s = h+t */
115 123 f1 = h * B0_lo + s * (v * (B1 + v * (B2 + v * (B3 + v * B4))));
116 124 t = f1 + t * B0;
117 125 h *= B0_hi;
118 - s = (double) ((float) (h + t));
126 + s = (double)((float)(h + t));
119 127 *w = t - (s - h);
120 128 return (s);
121 129 }
122 - if (ix < 0x00100000) { /* subnormal x */
130 +
131 + if (ix < 0x00100000) { /* subnormal x */
123 132 x *= two53;
124 133 n = -53;
125 134 ix = px[HIWORD];
126 135 }
136 +
127 137 /* LARGE N */
128 138 n += ((ix + 0x1000) >> 20) - 0x3ff;
129 - ix = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
139 + ix = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
130 140 px[HIWORD] = ix;
131 141 i = ix + 0x1000;
132 142 pz[HIWORD] = i & 0xffffe000;
133 143 pz[LOWORD] = 0;
134 144 qn = one / (x + z);
135 145 f = x - z;
136 146 s = f * qn;
137 - h = (double) ((float) s);
147 + h = (double)((float)s);
138 148 t = qn * ((f - (h + h) * z) - h * f);
139 149 j = (i >> 13) & 0x7f;
140 150 f = s * s;
141 151 t = t * A1 + h * A1_lo;
142 152 t += (s * f) * (A2 + f * A3);
143 153 qn = h * A1_hi;
144 154 s = n + _TBL_log2_hi[j];
145 155 h = qn + s;
146 156 t += _TBL_log2_lo[j] - ((h - s) - qn);
147 - f = (double) ((float) (h + t));
157 + f = (double)((float)(h + t));
148 158 *w = t - (f - h);
149 159 return (f);
150 160 }
151 161
152 162 extern const double _TBL_exp2_hi[], _TBL_exp2_lo[];
153 163 static const double /* poly app of 2^x-1 on [-1e-10,2^-7+1e-10] */
154 164 E1 = 6.931471805599453100674958533810346197328e-0001,
155 165 E2 = 2.402265069587779347846769151717493815979e-0001,
156 166 E3 = 5.550410866475410512631124892773937864699e-0002,
157 167 E4 = 9.618143209991026824853712740162451423355e-0003,
158 168 E5 = 1.333357676549940345096774122231849082991e-0003;
159 169
160 170 double
161 -pow(double x, double y) {
171 +pow(double x, double y)
172 +{
162 173 double z, ax;
163 174 double y1, y2, w1, w2;
164 175 int sbx, sby, j, k, yisint;
165 176 int hx, hy, ahx, ahy;
166 177 unsigned lx, ly;
167 - int *pz = (int *) &z;
178 + int *pz = (int *)&z;
168 179
169 - hx = ((int *) &x)[HIWORD];
170 - lx = ((unsigned *) &x)[LOWORD];
171 - hy = ((int *) &y)[HIWORD];
172 - ly = ((unsigned *) &y)[LOWORD];
180 + hx = ((int *)&x)[HIWORD];
181 + lx = ((unsigned *)&x)[LOWORD];
182 + hy = ((int *)&y)[HIWORD];
183 + ly = ((unsigned *)&y)[LOWORD];
173 184 ahx = hx & ~0x80000000;
174 185 ahy = hy & ~0x80000000;
175 - if ((ahy | ly) == 0) { /* y==zero */
186 +
187 + if ((ahy | ly) == 0) { /* y==zero */
176 188 if ((ahx | lx) == 0)
177 189 z = _SVID_libm_err(x, y, 20); /* +-0**+-0 */
178 190 else if ((ahx | (((lx | -lx) >> 31) & 1)) > 0x7ff00000)
179 191 z = _SVID_libm_err(x, y, 42); /* NaN**+-0 */
180 192 else
181 193 z = one; /* x**+-0 = 1 */
194 +
182 195 return (z);
183 - } else if (hx == 0x3ff00000 && lx == 0 &&
184 - (__xpg6 & _C99SUSv3_pow) != 0)
185 - return (one); /* C99: 1**anything = 1 */
186 - else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) ||
187 - ahy > 0x7ff00000 || (ahy == 0x7ff00000 && ly != 0))
196 + } else if (hx == 0x3ff00000 && lx == 0 && (__xpg6 & _C99SUSv3_pow) !=
197 + 0) {
198 + return (one); /* C99: 1**anything = 1 */
199 + } else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) || ahy >
200 + 0x7ff00000 || (ahy == 0x7ff00000 && ly != 0)) {
188 201 return (x * y); /* +-NaN return x*y; + -> * for Cheetah */
189 - /* includes Sun: 1**NaN = NaN */
190 - sbx = (unsigned) hx >> 31;
191 - sby = (unsigned) hy >> 31;
202 + }
203 +
204 + /* includes Sun: 1**NaN = NaN */
205 + sbx = (unsigned)hx >> 31;
206 + sby = (unsigned)hy >> 31;
192 207 ax = fabs(x);
193 208
194 209 /*
195 210 * determine if y is an odd int when x < 0
196 211 * yisint = 0 ... y is not an integer
197 212 * yisint = 1 ... y is an odd int
198 213 * yisint = 2 ... y is an even int
199 214 */
200 215 yisint = 0;
216 +
201 217 if (sbx) {
202 - if (ahy >= 0x43400000)
203 - yisint = 2; /* even integer y */
204 - else if (ahy >= 0x3ff00000) {
218 + if (ahy >= 0x43400000) {
219 + yisint = 2; /* even integer y */
220 + } else if (ahy >= 0x3ff00000) {
205 221 k = (ahy >> 20) - 0x3ff; /* exponent */
222 +
206 223 if (k > 20) {
207 224 j = ly >> (52 - k);
225 +
208 226 if ((j << (52 - k)) == ly)
209 227 yisint = 2 - (j & 1);
210 228 } else if (ly == 0) {
211 229 j = ahy >> (20 - k);
230 +
212 231 if ((j << (20 - k)) == ahy)
213 232 yisint = 2 - (j & 1);
214 233 }
215 234 }
216 235 }
236 +
217 237 /* special value of y */
218 238 if (ly == 0) {
219 239 if (ahy == 0x7ff00000) { /* y is +-inf */
220 240 if (((ahx - 0x3ff00000) | lx) == 0) {
221 241 if ((__xpg6 & _C99SUSv3_pow) != 0)
222 242 return (one);
223 - /* C99: (-1)**+-inf = 1 */
243 + /* C99: (-1)**+-inf = 1 */
224 244 else
225 245 return (y - y);
226 - /* Sun: (+-1)**+-inf = NaN */
227 - } else if (ahx >= 0x3ff00000)
228 - /* (|x|>1)**+,-inf = inf,0 */
246 +
247 + /* Sun: (+-1)**+-inf = NaN */
248 + } else if (ahx >= 0x3ff00000) {
249 + /* (|x|>1)**+,-inf = inf,0 */
229 250 return (sby == 0 ? y : zero);
230 - else /* (|x|<1)**-,+inf = inf,0 */
251 + } else { /* (|x|<1)**-,+inf = inf,0 */
231 252 return (sby != 0 ? -y : zero);
253 + }
232 254 }
255 +
233 256 if (ahy == 0x3ff00000) { /* y is +-1 */
234 - if (sby != 0) { /* y is -1 */
257 + if (sby != 0) { /* y is -1 */
235 258 if (x == zero) /* divided by zero */
236 259 return (_SVID_libm_err(x, y, 23));
237 260 else if (ahx < 0x40000 || ((ahx - 0x40000) |
238 - lx) == 0) /* overflow */
261 + lx) == 0) /* overflow */
239 262 return (_SVID_libm_err(x, y, 21));
240 263 else
241 264 return (one / x);
242 - } else
265 + } else {
243 266 return (x);
267 + }
244 268 }
245 - if (hy == 0x40000000) { /* y is 2 */
269 +
270 + if (hy == 0x40000000) { /* y is 2 */
246 271 if (ahx >= 0x5ff00000 && ahx < 0x7ff00000)
247 272 return (_SVID_libm_err(x, y, 21));
248 - /* x*x overflow */
273 + /* x*x overflow */
249 274 else if ((ahx < 0x1e56a09e && (ahx | lx) != 0) ||
250 - (ahx == 0x1e56a09e && lx < 0x667f3bcd))
275 + (ahx == 0x1e56a09e && lx < 0x667f3bcd))
251 276 return (_SVID_libm_err(x, y, 22));
252 - /* x*x underflow */
277 + /* x*x underflow */
253 278 else
254 279 return (x * x);
255 280 }
281 +
256 282 if (hy == 0x3fe00000) {
257 283 if (!((ahx | lx) == 0 || ((ahx - 0x7ff00000) | lx) ==
258 - 0 || sbx == 1))
284 + 0 || sbx == 1))
259 285 return (sqrt(x)); /* y is 0.5 and x > 0 */
260 286 }
261 287 }
288 +
262 289 /* special value of x */
263 290 if (lx == 0) {
264 291 if (ahx == 0x7ff00000 || ahx == 0 || ahx == 0x3ff00000) {
265 292 /* x is +-0,+-inf,-1 */
266 293 z = ax;
294 +
267 295 if (sby == 1) {
268 296 z = one / z; /* z = |x|**y */
297 +
269 298 if (ahx == 0)
270 299 return (_SVID_libm_err(x, y, 23));
271 300 }
301 +
272 302 if (sbx == 1) {
273 303 if (ahx == 0x3ff00000 && yisint == 0)
274 304 z = _SVID_libm_err(x, y, 24);
275 - /* neg**non-integral is NaN + invalid */
305 + /* neg**non-integral is NaN + invalid */
276 306 else if (yisint == 1)
277 307 z = -z; /* (x<0)**odd = -(|x|**odd) */
278 308 }
309 +
279 310 return (z);
280 311 }
281 312 }
313 +
282 314 /* (x<0)**(non-int) is NaN */
283 315 if (sbx == 1 && yisint == 0)
284 316 return (_SVID_libm_err(x, y, 24));
285 - /* Now ax is finite, y is finite */
286 - /* first compute log2(ax) = w1+w2, with 24 bits w1 */
317 +
318 + /*
319 + * Now ax is finite, y is finite
320 + * first compute log2(ax) = w1+w2, with 24 bits w1
321 + */
287 322 w1 = log2_x(ax, &w2);
288 323
289 324 /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
290 - if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 ||
291 - ahy <= 0x38100000) {
325 + if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 || ahy <=
326 + 0x38100000) {
292 327 /* no need to split if y is short or too large or too small */
293 328 y1 = y * w1;
294 329 y2 = y * w2;
295 330 } else {
296 - y1 = (double) ((float) y);
331 + y1 = (double)((float)y);
297 332 y2 = (y - y1) * w1 + y * w2;
298 333 y1 *= w1;
299 334 }
335 +
300 336 z = y1 + y2;
301 337 j = pz[HIWORD];
302 - if (j >= 0x40900000) { /* z >= 1024 */
303 - if (!(j == 0x40900000 && pz[LOWORD] == 0)) /* z > 1024 */
338 +
339 + if (j >= 0x40900000) { /* z >= 1024 */
340 + if (!(j == 0x40900000 && pz[LOWORD] == 0)) { /* z > 1024 */
304 341 return (_SVID_libm_err(x, y, 21)); /* overflow */
305 - else {
342 + } else {
306 343 w2 = y1 - z;
307 344 w2 += y2;
308 - /* rounded to inf */
345 +
346 + /* rounded to inf */
309 347 if (w2 >= -8.008566259537296567160e-17)
310 348 return (_SVID_libm_err(x, y, 21));
311 - /* overflow */
349 +
350 + /* overflow */
312 351 }
313 - } else if ((j & ~0x80000000) >= 0x4090cc00) { /* z <= -1075 */
314 - if (!(j == 0xc090cc00 && pz[LOWORD] == 0)) /* z < -1075 */
352 + } else if ((j & ~0x80000000) >= 0x4090cc00) { /* z <= -1075 */
353 + if (!(j == 0xc090cc00 && pz[LOWORD] == 0)) { /* z < -1075 */
315 354 return (_SVID_libm_err(x, y, 22)); /* underflow */
316 - else {
355 + } else {
317 356 w2 = y1 - z;
318 357 w2 += y2;
319 - if (w2 <= zero) /* underflow */
358 +
359 + if (w2 <= zero) /* underflow */
320 360 return (_SVID_libm_err(x, y, 22));
321 361 }
322 362 }
363 +
323 364 /*
324 365 * compute 2**(k+f[j]+g)
325 366 */
326 - k = (int) (z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
367 + k = (int)(z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
327 368 j = k & 63;
328 - w1 = y2 - ((double) k * 0.015625 - y1);
369 + w1 = y2 - ((double)k * 0.015625 - y1);
329 370 w2 = _TBL_exp2_hi[j];
330 - z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 *
331 - (E4 + w1 * E5))));
371 + z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 * (E4 +
372 + w1 * E5))));
332 373 z += w2;
333 374 k >>= 6;
375 +
334 376 if (k < -1021)
335 377 z = scalbn(z, k);
336 - else /* subnormal output */
378 + else /* subnormal output */
337 379 pz[HIWORD] += k << 20;
380 +
338 381 if (sbx == 1 && yisint == 1)
339 - z = -z; /* (-ve)**(odd int) */
382 + z = -z; /* (-ve)**(odd int) */
383 +
340 384 return (z);
341 385 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX