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 |
395 /* log(x) | 541 /* log(x) |
396 * Return the logrithm of x | 542 * Return the logrithm of x |
397 * | 543 * |
398 * Method : | 544 * Method : |
399 * 1. Argument Reduction: find k and f such that | 545 * 1. Argument Reduction: find k and f such that |
400 * x = 2^k * (1+f), | 546 * x = 2^k * (1+f), |
401 * where sqrt(2)/2 < 1+f < sqrt(2) . | 547 * where sqrt(2)/2 < 1+f < sqrt(2) . |
402 * | 548 * |
403 * 2. Approximation of log(1+f). | 549 * 2. Approximation of log(1+f). |
404 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) | 550 * 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... |
982 SET_HIGH_WORD(x, hx); | 1128 SET_HIGH_WORD(x, hx); |
983 SET_LOW_WORD(x, lx); | 1129 SET_LOW_WORD(x, lx); |
984 | 1130 |
985 double z = y * log10_2lo + ivln10 * log(x); | 1131 double z = y * log10_2lo + ivln10 * log(x); |
986 return z + y * log10_2hi; | 1132 return z + y * log10_2hi; |
987 } | 1133 } |
988 | 1134 |
989 } // namespace ieee754 | 1135 } // namespace ieee754 |
990 } // namespace base | 1136 } // namespace base |
991 } // namespace v8 | 1137 } // namespace v8 |
OLD | NEW |