| 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 |