OLD | NEW |
1 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm). | 1 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm). |
2 // | 2 // |
3 // ==================================================== | 3 // ==================================================== |
4 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 // | 5 // |
6 // Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 // Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 // Permission to use, copy, modify, and distribute this | 7 // Permission to use, copy, modify, and distribute this |
8 // software is freely granted, provided that this notice | 8 // software is freely granted, provided that this notice |
9 // is preserved. | 9 // is preserved. |
10 // ==================================================== | 10 // ==================================================== |
(...skipping 374 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
385 return z; /* atan(+,+) */ | 385 return z; /* atan(+,+) */ |
386 case 1: | 386 case 1: |
387 return -z; /* atan(-,+) */ | 387 return -z; /* atan(-,+) */ |
388 case 2: | 388 case 2: |
389 return pi - (z - pi_lo); /* atan(+,-) */ | 389 return pi - (z - pi_lo); /* atan(+,-) */ |
390 default: /* case 3 */ | 390 default: /* case 3 */ |
391 return (z - pi_lo) - pi; /* atan(-,-) */ | 391 return (z - pi_lo) - pi; /* atan(-,-) */ |
392 } | 392 } |
393 } | 393 } |
394 | 394 |
395 /* exp(x) | |
396 * Returns the exponential of x. | |
397 * | |
398 * Method | |
399 * 1. Argument reduction: | |
400 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. | |
401 * Given x, find r and integer k such that | |
402 * | |
403 * x = k*ln2 + r, |r| <= 0.5*ln2. | |
404 * | |
405 * Here r will be represented as r = hi-lo for better | |
406 * accuracy. | |
407 * | |
408 * 2. Approximation of exp(r) by a special rational function on | |
409 * the interval [0,0.34658]: | |
410 * Write | |
411 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... | |
412 * We use a special Remes algorithm on [0,0.34658] to generate | |
413 * a polynomial of degree 5 to approximate R. The maximum error | |
414 * of this polynomial approximation is bounded by 2**-59. In | |
415 * other words, | |
416 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 | |
417 * (where z=r*r, and the values of P1 to P5 are listed below) | |
418 * and | |
419 * | 5 | -59 | |
420 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 | |
421 * | | | |
422 * The computation of exp(r) thus becomes | |
423 * 2*r | |
424 * exp(r) = 1 + ------- | |
425 * R - r | |
426 * r*R1(r) | |
427 * = 1 + r + ----------- (for better accuracy) | |
428 * 2 - R1(r) | |
429 * where | |
430 * 2 4 10 | |
431 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). | |
432 * | |
433 * 3. Scale back to obtain exp(x): | |
434 * From step 1, we have | |
435 * exp(x) = 2^k * exp(r) | |
436 * | |
437 * Special cases: | |
438 * exp(INF) is INF, exp(NaN) is NaN; | |
439 * exp(-INF) is 0, and | |
440 * for finite argument, only exp(0)=1 is exact. | |
441 * | |
442 * Accuracy: | |
443 * according to an error analysis, the error is always less than | |
444 * 1 ulp (unit in the last place). | |
445 * | |
446 * Misc. info. | |
447 * For IEEE double | |
448 * if x > 7.09782712893383973096e+02 then exp(x) overflow | |
449 * if x < -7.45133219101941108420e+02 then exp(x) underflow | |
450 * | |
451 * Constants: | |
452 * The hexadecimal values are the intended ones for the following | |
453 * constants. The decimal values may be used, provided that the | |
454 * compiler will convert from decimal to binary accurately enough | |
455 * to produce the hexadecimal values shown. | |
456 */ | |
457 double exp(double x) { | |
458 static const double | |
459 one = 1.0, | |
460 halF[2] = {0.5, -0.5}, | |
461 o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ | |
462 u_threshold = -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ | |
463 ln2HI[2] = {6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ | |
464 -6.93147180369123816490e-01}, /* 0xbfe62e42, 0xfee00000 */ | |
465 ln2LO[2] = {1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ | |
466 -1.90821492927058770002e-10}, /* 0xbdea39ef, 0x35793c76 */ | |
467 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ | |
468 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ | |
469 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ | |
470 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ | |
471 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ | |
472 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ | |
473 | |
474 static volatile double | |
475 huge = 1.0e+300, | |
476 twom1000 = 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ | |
477 two1023 = 8.988465674311579539e307; /* 0x1p1023 */ | |
478 | |
479 double y, hi = 0.0, lo = 0.0, c, t, twopk; | |
480 int32_t k = 0, xsb; | |
481 u_int32_t hx; | |
482 | |
483 GET_HIGH_WORD(hx, x); | |
484 xsb = (hx >> 31) & 1; /* sign bit of x */ | |
485 hx &= 0x7fffffff; /* high word of |x| */ | |
486 | |
487 /* filter out non-finite argument */ | |
488 if (hx >= 0x40862E42) { /* if |x|>=709.78... */ | |
489 if (hx >= 0x7ff00000) { | |
490 u_int32_t lx; | |
491 GET_LOW_WORD(lx, x); | |
492 if (((hx & 0xfffff) | lx) != 0) | |
493 return x + x; /* NaN */ | |
494 else | |
495 return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */ | |
496 } | |
497 if (x > o_threshold) return huge * huge; /* overflow */ | |
498 if (x < u_threshold) return twom1000 * twom1000; /* underflow */ | |
499 } | |
500 | |
501 /* argument reduction */ | |
502 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ | |
503 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ | |
504 hi = x - ln2HI[xsb]; | |
505 lo = ln2LO[xsb]; | |
506 k = 1 - xsb - xsb; | |
507 } else { | |
508 k = static_cast<int>(invln2 * x + halF[xsb]); | |
509 t = k; | |
510 hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */ | |
511 lo = t * ln2LO[0]; | |
512 } | |
513 STRICT_ASSIGN(double, x, hi - lo); | |
514 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */ | |
515 if (huge + x > one) return one + x; /* trigger inexact */ | |
516 } else { | |
517 k = 0; | |
518 } | |
519 | |
520 /* x is now in primary range */ | |
521 t = x * x; | |
522 if (k >= -1021) { | |
523 INSERT_WORDS(twopk, 0x3ff00000 + (k << 20), 0); | |
524 } else { | |
525 INSERT_WORDS(twopk, 0x3ff00000 + ((k + 1000) << 20), 0); | |
526 } | |
527 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); | |
528 if (k == 0) { | |
529 return one - ((x * c) / (c - 2.0) - x); | |
530 } else { | |
531 y = one - ((lo - (x * c) / (2.0 - c)) - hi); | |
532 } | |
533 if (k >= -1021) { | |
534 if (k == 1024) return y * 2.0 * two1023; | |
535 return y * twopk; | |
536 } else { | |
537 return y * twopk * twom1000; | |
538 } | |
539 } | |
540 | |
541 /* log(x) | 395 /* log(x) |
542 * Return the logrithm of x | 396 * Return the logrithm of x |
543 * | 397 * |
544 * Method : | 398 * Method : |
545 * 1. Argument Reduction: find k and f such that | 399 * 1. Argument Reduction: find k and f such that |
546 * x = 2^k * (1+f), | 400 * x = 2^k * (1+f), |
547 * where sqrt(2)/2 < 1+f < sqrt(2) . | 401 * where sqrt(2)/2 < 1+f < sqrt(2) . |
548 * | 402 * |
549 * 2. Approximation of log(1+f). | 403 * 2. Approximation of log(1+f). |
550 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) | 404 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) |
(...skipping 577 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1128 SET_HIGH_WORD(x, hx); | 982 SET_HIGH_WORD(x, hx); |
1129 SET_LOW_WORD(x, lx); | 983 SET_LOW_WORD(x, lx); |
1130 | 984 |
1131 double z = y * log10_2lo + ivln10 * log(x); | 985 double z = y * log10_2lo + ivln10 * log(x); |
1132 return z + y * log10_2hi; | 986 return z + y * log10_2hi; |
1133 } | 987 } |
1134 | 988 |
1135 } // namespace ieee754 | 989 } // namespace ieee754 |
1136 } // namespace base | 990 } // namespace base |
1137 } // namespace v8 | 991 } // namespace v8 |
OLD | NEW |