1 /* Copyright JS Foundation and other contributors, http://js.foundation
2 *
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 * This file is based on work under the following copyright and permission
16 * notice:
17 *
18 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
19 *
20 * Permission to use, copy, modify, and distribute this
21 * software is freely granted, provided that this notice
22 * is preserved.
23 *
24 * @(#)e_pow.c 1.5 04/04/22
25 */
26
27 #include "jerry-libm-internal.h"
28
29 /* pow(x,y) return x**y
30 *
31 * n
32 * Method: Let x = 2 * (1+f)
33 * 1. Compute and return log2(x) in two pieces:
34 * log2(x) = w1 + w2,
35 * where w1 has 53-24 = 29 bit trailing zeros.
36 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
37 * arithmetic, where |y'|<=0.5.
38 * 3. Return x**y = 2**n*exp(y'*log2)
39 *
40 * Special cases:
41 * 0. +1 ** (anything) is 1
42 * 1. (anything) ** 0 is 1
43 * 2. (anything) ** 1 is itself
44 * 3. (anything) ** NAN is NAN
45 * 4. NAN ** (anything except 0) is NAN
46 * 5. +-(|x| > 1) ** +INF is +INF
47 * 6. +-(|x| > 1) ** -INF is +0
48 * 7. +-(|x| < 1) ** +INF is +0
49 * 8. +-(|x| < 1) ** -INF is +INF
50 * 9. -1 ** +-INF is 1
51 * 10. +0 ** (+anything except 0, NAN) is +0
52 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
53 * 12. +0 ** (-anything except 0, NAN) is +INF
54 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
55 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
56 * 15. +INF ** (+anything except 0,NAN) is +INF
57 * 16. +INF ** (-anything except 0,NAN) is +0
58 * 17. -INF ** (anything) = -0 ** (-anything)
59 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
60 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
61 *
62 * Accuracy:
63 * pow(x,y) returns x**y nearly rounded. In particular
64 * pow(integer,integer)
65 * always returns the correct integer provided it is
66 * representable.
67 *
68 * Constants:
69 * The hexadecimal values are the intended ones for the following
70 * constants. The decimal values may be used, provided that the
71 * compiler will convert from decimal to binary accurately enough
72 * to produce the hexadecimal values shown.
73 */
74
75 static const double bp[] =
76 {
77 1.0,
78 1.5,
79 };
80 static const double dp_h[] =
81 {
82 0.0,
83 5.84962487220764160156e-01, /* 0x3FE2B803, 0x40000000 */
84 };
85 static const double dp_l[] =
86 {
87 0.0,
88 1.35003920212974897128e-08, /* 0x3E4CFDEB, 0x43CFD006 */
89 };
90
91 #define zero 0.0
92 #define one 1.0
93 #define two 2.0
94 #define two53 9007199254740992.0 /* 0x43400000, 0x00000000 */
95 #define huge 1.0e300
96 #define tiny 1.0e-300
97 /* poly coefs for (3/2) * (log(x) - 2s - 2/3 * s**3 */
98 #define L1 5.99999999999994648725e-01 /* 0x3FE33333, 0x33333303 */
99 #define L2 4.28571428578550184252e-01 /* 0x3FDB6DB6, 0xDB6FABFF */
100 #define L3 3.33333329818377432918e-01 /* 0x3FD55555, 0x518F264D */
101 #define L4 2.72728123808534006489e-01 /* 0x3FD17460, 0xA91D4101 */
102 #define L5 2.30660745775561754067e-01 /* 0x3FCD864A, 0x93C9DB65 */
103 #define L6 2.06975017800338417784e-01 /* 0x3FCA7E28, 0x4A454EEF */
104 #define P1 1.66666666666666019037e-01 /* 0x3FC55555, 0x5555553E */
105 #define P2 -2.77777777770155933842e-03 /* 0xBF66C16C, 0x16BEBD93 */
106 #define P3 6.61375632143793436117e-05 /* 0x3F11566A, 0xAF25DE2C */
107 #define P4 -1.65339022054652515390e-06 /* 0xBEBBBD41, 0xC5D26BF1 */
108 #define P5 4.13813679705723846039e-08 /* 0x3E663769, 0x72BEA4D0 */
109 #define lg2 6.93147180559945286227e-01 /* 0x3FE62E42, 0xFEFA39EF */
110 #define lg2_h 6.93147182464599609375e-01 /* 0x3FE62E43, 0x00000000 */
111 #define lg2_l -1.90465429995776804525e-09 /* 0xBE205C61, 0x0CA86C39 */
112 #define ovt 8.0085662595372944372e-0017 /* -(1024-log2(ovfl+.5ulp)) */
113 #define cp 9.61796693925975554329e-01 /* 0x3FEEC709, 0xDC3A03FD = 2 / (3 ln2) */
114 #define cp_h 9.61796700954437255859e-01 /* 0x3FEEC709, 0xE0000000 = (float) cp */
115 #define cp_l -7.02846165095275826516e-09 /* 0xBE3E2FE0, 0x145B01F5 = tail of cp_h */
116 #define ivln2 1.44269504088896338700e+00 /* 0x3FF71547, 0x652B82FE = 1 / ln2 */
117 #define ivln2_h 1.44269502162933349609e+00 /* 0x3FF71547, 0x60000000 = 24b 1 / ln2 */
118 #define ivln2_l 1.92596299112661746887e-08 /* 0x3E54AE0B, 0xF85DDF44 = 1 / ln2 tail */
119
120 double
pow(double x,double y)121 pow (double x, double y)
122 {
123 double_accessor t1, ax, p_h, y1, t, z;
124 double z_h, z_l, p_l;
125 double t2, r, s, u, v, w;
126 int i, j, k, yisint, n;
127 int hx, hy, ix, iy;
128 unsigned lx, ly;
129
130 hx = __HI (x);
131 lx = __LO (x);
132 hy = __HI (y);
133 ly = __LO (y);
134 ix = hx & 0x7fffffff;
135 iy = hy & 0x7fffffff;
136
137 /* x == one: 1**y = 1 */
138 if (((hx - 0x3ff00000) | lx) == 0)
139 {
140 return one;
141 }
142
143 /* y == zero: x**0 = 1 */
144 if ((iy | ly) == 0)
145 {
146 return one;
147 }
148
149 /* +-NaN return x + y */
150 if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) || iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
151 {
152 return x + y;
153 }
154
155 /* determine if y is an odd int when x < 0
156 * yisint = 0 ... y is not an integer
157 * yisint = 1 ... y is an odd int
158 * yisint = 2 ... y is an even int
159 */
160 yisint = 0;
161 if (hx < 0)
162 {
163 if (iy >= 0x43400000) /* even integer y */
164 {
165 yisint = 2;
166 }
167 else if (iy >= 0x3ff00000)
168 {
169 k = (iy >> 20) - 0x3ff; /* exponent */
170 if (k > 20)
171 {
172 j = ly >> (52 - k);
173 if ((j << (52 - k)) == ly)
174 {
175 yisint = 2 - (j & 1);
176 }
177 }
178 else if (ly == 0)
179 {
180 j = iy >> (20 - k);
181 if ((j << (20 - k)) == iy)
182 {
183 yisint = 2 - (j & 1);
184 }
185 }
186 }
187 }
188
189 /* special value of y */
190 if (ly == 0)
191 {
192 if (iy == 0x7ff00000) /* y is +-inf */
193 {
194 if (((ix - 0x3ff00000) | lx) == 0) /* +-1**+-inf is 1 */
195 {
196 return one;
197 }
198 else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
199 {
200 return (hy >= 0) ? y : zero;
201 }
202 else /* (|x|<1)**-,+inf = inf,0 */
203 {
204 return (hy < 0) ? -y : zero;
205 }
206 }
207 if (iy == 0x3ff00000) /* y is +-1 */
208 {
209 if (hy < 0)
210 {
211 return one / x;
212 }
213 else
214 {
215 return x;
216 }
217 }
218 if (hy == 0x40000000) /* y is 2 */
219 {
220 return x * x;
221 }
222 if (hy == 0x3fe00000) /* y is 0.5 */
223 {
224 if (hx >= 0) /* x >= +0 */
225 {
226 return sqrt (x);
227 }
228 }
229 }
230
231 ax.dbl = fabs (x);
232 /* special value of x */
233 if (lx == 0)
234 {
235 if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000)
236 {
237 z.dbl = ax.dbl; /* x is +-0,+-inf,+-1 */
238 if (hy < 0)
239 {
240 z.dbl = one / z.dbl; /* z = (1 / |x|) */
241 }
242 if (hx < 0)
243 {
244 if (((ix - 0x3ff00000) | yisint) == 0)
245 {
246 z.dbl = NAN; /* (-1)**non-int is NaN */
247 }
248 else if (yisint == 1)
249 {
250 z.dbl = -z.dbl; /* (x<0)**odd = -(|x|**odd) */
251 }
252 }
253 return z.dbl;
254 }
255 }
256
257 n = (hx < 0) ? 0 : 1;
258
259 /* (x<0)**(non-int) is NaN */
260 if ((n | yisint) == 0)
261 {
262 return NAN;
263 }
264
265 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
266 if ((n | (yisint - 1)) == 0)
267 {
268 s = -one; /* (-ve)**(odd int) */
269 }
270
271 /* |y| is huge */
272 if (iy > 0x41e00000) /* if |y| > 2**31 */
273 {
274 if (iy > 0x43f00000) /* if |y| > 2**64, must o/uflow */
275 {
276 if (ix <= 0x3fefffff)
277 {
278 return (hy < 0) ? huge * huge : tiny * tiny;
279 }
280 if (ix >= 0x3ff00000)
281 {
282 return (hy > 0) ? huge * huge : tiny * tiny;
283 }
284 }
285 /* over/underflow if x is not close to one */
286 if (ix < 0x3fefffff)
287 {
288 return (hy < 0) ? s * huge * huge : s * tiny * tiny;
289 }
290 if (ix > 0x3ff00000)
291 {
292 return (hy > 0) ? s * huge * huge : s * tiny * tiny;
293 }
294 /* now |1 - x| is tiny <= 2**-20, suffice to compute
295 log(x) by x - x^2 / 2 + x^3 / 3 - x^4 / 4 */
296 t.dbl = ax.dbl - one; /* t has 20 trailing zeros */
297 w = (t.dbl * t.dbl) * (0.5 - t.dbl * (0.3333333333333333333333 - t.dbl * 0.25));
298 u = ivln2_h * t.dbl; /* ivln2_h has 21 sig. bits */
299 v = t.dbl * ivln2_l - w * ivln2;
300 t1.dbl = u + v;
301 t1.as_int.lo = 0;
302 t2 = v - (t1.dbl - u);
303 }
304 else
305 {
306 double_accessor s_h, t_h;
307 double ss, s2, s_l, t_l;
308
309 n = 0;
310 /* take care subnormal number */
311 if (ix < 0x00100000)
312 {
313 ax.dbl *= two53;
314 n -= 53;
315 ix = ax.as_int.hi;
316 }
317 n += ((ix) >> 20) - 0x3ff;
318 j = ix & 0x000fffff;
319 /* determine interval */
320 ix = j | 0x3ff00000; /* normalize ix */
321 if (j <= 0x3988E) /* |x| < sqrt(3/2) */
322 {
323 k = 0;
324 }
325 else if (j < 0xBB67A) /* |x| < sqrt(3) */
326 {
327 k = 1;
328 }
329 else
330 {
331 k = 0;
332 n += 1;
333 ix -= 0x00100000;
334 }
335 ax.as_int.hi = ix;
336
337 /* compute ss = s_h + s_l = (x - 1) / (x + 1) or (x - 1.5) / (x + 1.5) */
338 u = ax.dbl - bp[k]; /* bp[0] = 1.0, bp[1] = 1.5 */
339 v = one / (ax.dbl + bp[k]);
340 ss = u * v;
341 s_h.dbl = ss;
342 s_h.as_int.lo = 0;
343 /* t_h = ax + bp[k] High */
344 t_h.dbl = zero;
345 t_h.as_int.hi = ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18);
346 t_l = ax.dbl - (t_h.dbl - bp[k]);
347 s_l = v * ((u - s_h.dbl * t_h.dbl) - s_h.dbl * t_l);
348 /* compute log(ax) */
349 s2 = ss * ss;
350 r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
351 r += s_l * (s_h.dbl + ss);
352 s2 = s_h.dbl * s_h.dbl;
353 t_h.dbl = 3.0 + s2 + r;
354 t_h.as_int.lo = 0;
355 t_l = r - ((t_h.dbl - 3.0) - s2);
356 /* u + v = ss * (1 + ...) */
357 u = s_h.dbl * t_h.dbl;
358 v = s_l * t_h.dbl + t_l * ss;
359 /* 2 / (3 * log2) * (ss + ...) */
360 p_h.dbl = u + v;
361 p_h.as_int.lo = 0;
362 p_l = v - (p_h.dbl - u);
363 z_h = cp_h * p_h.dbl; /* cp_h + cp_l = 2 / (3 * log2) */
364 z_l = cp_l * p_h.dbl + p_l * cp + dp_l[k];
365 /* log2(ax) = (ss + ...) * 2 / (3 * log2) = n + dp_h + z_h + z_l */
366 t.dbl = (double) n;
367 t1.dbl = (((z_h + z_l) + dp_h[k]) + t.dbl);
368 t1.as_int.lo = 0;
369 t2 = z_l - (((t1.dbl - t.dbl) - dp_h[k]) - z_h);
370 }
371
372 /* split up y into y1 + y2 and compute (y1 + y2) * (t1 + t2) */
373 y1.dbl = y;
374 y1.as_int.lo = 0;
375 p_l = (y - y1.dbl) * t1.dbl + y * t2;
376 p_h.dbl = y1.dbl * t1.dbl;
377 z.dbl = p_l + p_h.dbl;
378 j = z.as_int.hi;
379 i = z.as_int.lo;
380 if (j >= 0x40900000) /* z >= 1024 */
381 {
382 if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
383 {
384 return s * huge * huge; /* overflow */
385 }
386 else
387 {
388 if (p_l + ovt > z.dbl - p_h.dbl)
389 {
390 return s * huge * huge; /* overflow */
391 }
392 }
393 }
394 else if ((j & 0x7fffffff) >= 0x4090cc00) /* z <= -1075 */
395 {
396 if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */
397 {
398 return s * tiny * tiny; /* underflow */
399 }
400 else
401 {
402 if (p_l <= z.dbl - p_h.dbl)
403 {
404 return s * tiny * tiny; /* underflow */
405 }
406 }
407 }
408 /*
409 * compute 2**(p_h + p_l)
410 */
411 i = j & 0x7fffffff;
412 k = (i >> 20) - 0x3ff;
413 n = 0;
414 if (i > 0x3fe00000) /* if |z| > 0.5, set n = [z + 0.5] */
415 {
416 n = j + (0x00100000 >> (k + 1));
417 k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
418 t.dbl = zero;
419 t.as_int.hi = (n & ~(0x000fffff >> k));
420 n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
421 if (j < 0)
422 {
423 n = -n;
424 }
425 p_h.dbl -= t.dbl;
426 }
427 t.dbl = p_l + p_h.dbl;
428 t.as_int.lo = 0;
429 u = t.dbl * lg2_h;
430 v = (p_l - (t.dbl - p_h.dbl)) * lg2 + t.dbl * lg2_l;
431 z.dbl = u + v;
432 w = v - (z.dbl - u);
433 t.dbl = z.dbl * z.dbl;
434 t1.dbl = z.dbl - t.dbl * (P1 + t.dbl * (P2 + t.dbl * (P3 + t.dbl * (P4 + t.dbl * P5))));
435 r = (z.dbl * t1.dbl) / (t1.dbl - two) - (w + z.dbl * w);
436 z.dbl = one - (r - z.dbl);
437 j = z.as_int.hi;
438 j += (n << 20);
439 if ((j >> 20) <= 0) /* subnormal output */
440 {
441 z.dbl = scalbn (z.dbl, n);
442 }
443 else
444 {
445 z.as_int.hi += (n << 20);
446 }
447 return s * z.dbl;
448 } /* pow */
449
450 #undef zero
451 #undef one
452 #undef two
453 #undef two53
454 #undef huge
455 #undef tiny
456 #undef L1
457 #undef L2
458 #undef L3
459 #undef L4
460 #undef L5
461 #undef L6
462 #undef P1
463 #undef P2
464 #undef P3
465 #undef P4
466 #undef P5
467 #undef lg2
468 #undef lg2_h
469 #undef lg2_l
470 #undef ovt
471 #undef cp
472 #undef cp_h
473 #undef cp_l
474 #undef ivln2
475 #undef ivln2_h
476 #undef ivln2_l
477