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