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-2004 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 // ==================================================== |
11 // | 11 // |
12 // The original source code covered by the above license above has been | 12 // The original source code covered by the above license above has been |
13 // modified significantly by Google Inc. | 13 // modified significantly by Google Inc. |
14 // Copyright 2014 the V8 project authors. All rights reserved. | 14 // Copyright 2014 the V8 project authors. All rights reserved. |
15 // | 15 // |
16 // The following is a straightforward translation of fdlibm routines | 16 // The following is a straightforward translation of fdlibm routines |
17 // by Raymond Toy (rtoy@google.com). | 17 // by Raymond Toy (rtoy@google.com). |
18 | 18 |
19 | 19 // Double constants that do not have empty lower 32 bits are found in fdlibm.cc |
20 var kMath; // Initialized to a Float64Array during genesis and is not writable. | 20 // and exposed through kMath as typed array. We assume the compiler to convert |
21 // from decimal to binary accurately enough to produce the intended values. | |
22 // kMath is initialized to a Float64Array during genesis and not writable. | |
23 var kMath; | |
21 | 24 |
22 const INVPIO2 = kMath[0]; | 25 const INVPIO2 = kMath[0]; |
23 const PIO2_1 = kMath[1]; | 26 const PIO2_1 = kMath[1]; |
24 const PIO2_1T = kMath[2]; | 27 const PIO2_1T = kMath[2]; |
25 const PIO2_2 = kMath[3]; | 28 const PIO2_2 = kMath[3]; |
26 const PIO2_2T = kMath[4]; | 29 const PIO2_2T = kMath[4]; |
27 const PIO2_3 = kMath[5]; | 30 const PIO2_3 = kMath[5]; |
28 const PIO2_3T = kMath[6]; | 31 const PIO2_3T = kMath[6]; |
29 const PIO4 = kMath[32]; | 32 const PIO4 = kMath[32]; |
30 const PIO4LO = kMath[33]; | 33 const PIO4LO = kMath[33]; |
(...skipping 369 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
400 // Special cases: | 403 // Special cases: |
401 // log1p(x) is NaN with signal if x < -1 (including -INF) ; | 404 // log1p(x) is NaN with signal if x < -1 (including -INF) ; |
402 // log1p(+INF) is +INF; log1p(-1) is -INF with signal; | 405 // log1p(+INF) is +INF; log1p(-1) is -INF with signal; |
403 // log1p(NaN) is that NaN with no signal. | 406 // log1p(NaN) is that NaN with no signal. |
404 // | 407 // |
405 // Accuracy: | 408 // Accuracy: |
406 // according to an error analysis, the error is always less than | 409 // according to an error analysis, the error is always less than |
407 // 1 ulp (unit in the last place). | 410 // 1 ulp (unit in the last place). |
408 // | 411 // |
409 // Constants: | 412 // Constants: |
410 // The hexadecimal values are the intended ones for the following | 413 // Constants are found in fdlibm.cc. We assume the C++ compiler to convert |
411 // constants. The decimal values may be used, provided that the | 414 // from decimal to binary accurately enough to produce the intended values. |
412 // compiler will convert from decimal to binary accurately enough | |
413 // to produce the hexadecimal values shown. | |
414 // | 415 // |
415 // Note: Assuming log() return accurate answer, the following | 416 // Note: Assuming log() return accurate answer, the following |
416 // algorithm can be used to compute log1p(x) to within a few ULP: | 417 // algorithm can be used to compute log1p(x) to within a few ULP: |
417 // | 418 // |
418 // u = 1+x; | 419 // u = 1+x; |
419 // if (u==1.0) return x ; else | 420 // if (u==1.0) return x ; else |
420 // return log(u)*(x/(u-1.0)); | 421 // return log(u)*(x/(u-1.0)); |
421 // | 422 // |
422 // See HP-15C Advanced Functions Handbook, p.193. | 423 // See HP-15C Advanced Functions Handbook, p.193. |
423 // | 424 // |
424 const LN2_HI = kMath[34]; | 425 const LN2_HI = kMath[34]; |
425 const LN2_LO = kMath[35]; | 426 const LN2_LO = kMath[35]; |
426 const TWO54 = kMath[36]; | 427 const TWO54 = kMath[36]; |
427 const TWO_THIRD = kMath[37]; | 428 const TWO_THIRD = kMath[37]; |
428 macro KLOGP1(x) | 429 macro KLOG1P(x) |
Raymond Toy
2014/08/13 16:25:17
What changed here? And how is this related to expm
| |
429 (kMath[38+x]) | 430 (kMath[38+x]) |
430 endmacro | 431 endmacro |
431 | 432 |
432 function MathLog1p(x) { | 433 function MathLog1p(x) { |
433 x = x * 1; // Convert to number. | 434 x = x * 1; // Convert to number. |
434 var hx = %_DoubleHi(x); | 435 var hx = %_DoubleHi(x); |
435 var ax = hx & 0x7fffffff; | 436 var ax = hx & 0x7fffffff; |
436 var k = 1; | 437 var k = 1; |
437 var f = x; | 438 var f = x; |
438 var hu = 1; | 439 var hu = 1; |
(...skipping 61 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
500 var R = hfsq * (1 - TWO_THIRD * f); | 501 var R = hfsq * (1 - TWO_THIRD * f); |
501 if (k === 0) { | 502 if (k === 0) { |
502 return f - R; | 503 return f - R; |
503 } else { | 504 } else { |
504 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); | 505 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); |
505 } | 506 } |
506 } | 507 } |
507 | 508 |
508 var s = f / (2 + f); | 509 var s = f / (2 + f); |
509 var z = s * s; | 510 var z = s * s; |
510 var R = z * (KLOGP1(0) + z * (KLOGP1(1) + z * | 511 var R = z * (KLOG1P(0) + z * (KLOG1P(1) + z * |
511 (KLOGP1(2) + z * (KLOGP1(3) + z * | 512 (KLOG1P(2) + z * (KLOG1P(3) + z * |
512 (KLOGP1(4) + z * (KLOGP1(5) + z * KLOGP1(6))))))); | 513 (KLOG1P(4) + z * (KLOG1P(5) + z * KLOG1P(6))))))); |
Raymond Toy
2014/08/13 16:25:17
This seems unrelated to expm1. And I can't see wha
Yang
2014/08/20 14:18:24
it used to be KLOGP1. I realized that it should be
Raymond Toy
2014/08/20 16:18:20
Sounds good. Sorry I missed that in the review of
| |
513 if (k === 0) { | 514 if (k === 0) { |
514 return f - (hfsq - s * (hfsq + R)); | 515 return f - (hfsq - s * (hfsq + R)); |
515 } else { | 516 } else { |
516 return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); | 517 return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); |
517 } | 518 } |
518 } | 519 } |
520 | |
521 // ES6 draft 09-27-13, section 20.2.2.14. | |
522 // Math.expm1 | |
523 // Returns exp(x)-1, the exponential of x minus 1. | |
524 // | |
525 // Method | |
526 // 1. Argument reduction: | |
527 // Given x, find r and integer k such that | |
528 // | |
529 // x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 | |
530 // | |
531 // Here a correction term c will be computed to compensate | |
532 // the error in r when rounded to a floating-point number. | |
533 // | |
534 // 2. Approximating expm1(r) by a special rational function on | |
535 // the interval [0,0.34658]: | |
536 // Since | |
537 // r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ... | |
538 // we define R1(r*r) by | |
539 // r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r) | |
540 // That is, | |
541 // R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) | |
542 // = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) | |
543 // = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... | |
544 // We use a special Remes algorithm on [0,0.347] to generate | |
545 // a polynomial of degree 5 in r*r to approximate R1. The | |
546 // maximum error of this polynomial approximation is bounded | |
547 // by 2**-61. In other words, | |
548 // R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 | |
549 // where Q1 = -1.6666666666666567384E-2, | |
550 // Q2 = 3.9682539681370365873E-4, | |
551 // Q3 = -9.9206344733435987357E-6, | |
552 // Q4 = 2.5051361420808517002E-7, | |
553 // Q5 = -6.2843505682382617102E-9; | |
554 // (where z=r*r, and the values of Q1 to Q5 are listed below) | |
555 // with error bounded by | |
556 // | 5 | -61 | |
557 // | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 | |
558 // | | | |
559 // | |
560 // expm1(r) = exp(r)-1 is then computed by the following | |
561 // specific way which minimize the accumulation rounding error: | |
562 // 2 3 | |
563 // r r [ 3 - (R1 + R1*r/2) ] | |
564 // expm1(r) = r + --- + --- * [--------------------] | |
565 // 2 2 [ 6 - r*(3 - R1*r/2) ] | |
566 // | |
567 // To compensate the error in the argument reduction, we use | |
568 // expm1(r+c) = expm1(r) + c + expm1(r)*c | |
569 // ~ expm1(r) + c + r*c | |
570 // Thus c+r*c will be added in as the correction terms for | |
571 // expm1(r+c). Now rearrange the term to avoid optimization | |
572 // screw up: | |
573 // ( 2 2 ) | |
574 // ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) | |
575 // expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) | |
576 // ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) | |
577 // ( ) | |
578 // | |
579 // = r - E | |
580 // 3. Scale back to obtain expm1(x): | |
581 // From step 1, we have | |
582 // expm1(x) = either 2^k*[expm1(r)+1] - 1 | |
583 // = or 2^k*[expm1(r) + (1-2^-k)] | |
584 // 4. Implementation notes: | |
585 // (A). To save one multiplication, we scale the coefficient Qi | |
586 // to Qi*2^i, and replace z by (x^2)/2. | |
587 // (B). To achieve maximum accuracy, we compute expm1(x) by | |
588 // (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) | |
589 // (ii) if k=0, return r-E | |
590 // (iii) if k=-1, return 0.5*(r-E)-0.5 | |
591 // (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) | |
592 // else return 1.0+2.0*(r-E); | |
593 // (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) | |
594 // (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else | |
595 // (vii) return 2^k(1-((E+2^-k)-r)) | |
596 // | |
597 // Special cases: | |
598 // expm1(INF) is INF, expm1(NaN) is NaN; | |
599 // expm1(-INF) is -1, and | |
600 // for finite argument, only expm1(0)=0 is exact. | |
601 // | |
602 // Accuracy: | |
603 // according to an error analysis, the error is always less than | |
604 // 1 ulp (unit in the last place). | |
605 // | |
606 // Misc. info. | |
607 // For IEEE double | |
608 // if x > 7.09782712893383973096e+02 then expm1(x) overflow | |
609 // | |
610 const KEXPM1_OVERFLOW = kMath[45]; | |
611 const INVLN2 = kMath[46]; | |
612 macro KEXPM1(x) | |
613 (kMath[47+x]) | |
614 endmacro | |
615 | |
616 function MathExpm1(x) { | |
617 x = x * 1; // Convert to number. | |
618 var y; | |
619 var hi; | |
620 var lo; | |
621 var k; | |
622 var t; | |
623 var c; | |
624 | |
625 var hx = %_DoubleHi(x); | |
626 var xsb = hx & 0x80000000; // Sign bit of x | |
627 var y = (xsb === 0) ? x : -x; // y = |x| | |
628 hx &= 0x7fffffff; // High word of |x| | |
629 | |
630 // Filter out huge and non-finite argument | |
631 if (hx >= 0x4043687a) { // if |x| ~=> 56 * ln2 | |
632 if (hx >= 0x40862e42) { // if |x| >= 709.78 | |
633 if (hx >= 0x7ff00000) { | |
634 // expm1(inf) = inf; expm1(-inf) = -1; expm1(nan) = nan; | |
635 return (x === -INFINITY) ? -1 : x; | |
636 } | |
637 if (x > KEXPM1_OVERFLOW) return INFINITY; // Overflow | |
638 } | |
639 if (xsb != 0) return -1; // x < -56 * ln2, return -1. | |
640 } | |
641 | |
642 // Argument reduction | |
643 if (hx > 0x3fd62e42) { // if |x| > 0.5 * ln2 | |
644 if (hx < 0x3ff0a2b2) { // and |x| < 1.5 * ln2 | |
645 if (xsb === 0) { | |
646 hi = x - LN2_HI; | |
647 lo = LN2_LO; | |
648 k = 1; | |
649 } else { | |
650 hi = x + LN2_HI; | |
651 lo = -LN2_LO; | |
652 k = -1; | |
653 } | |
654 } else { | |
655 k = (INVLN2 * x + ((xsb === 0) ? 0.5 : -0.5)) | 0; | |
656 t = k; | |
657 // t * ln2_hi is exact here. | |
658 hi = x - t * LN2_HI; | |
659 lo = t * LN2_LO; | |
660 } | |
661 x = hi - lo; | |
662 c = (hi - x) - lo; | |
663 } else if (hx < 0x3c900000) { | |
664 // When |x| < 2^-54, we can return x. | |
665 return x; | |
666 } else { | |
667 // Fall through. | |
668 k = 0; | |
669 } | |
670 | |
671 // x is now in primary range | |
672 var hfx = 0.5 * x; | |
673 var hxs = x * hfx; | |
674 var r1 = 1 + hxs * (KEXPM1(0) + hxs * (KEXPM1(1) + hxs * | |
675 (KEXPM1(2) + hxs * (KEXPM1(3) + hxs * KEXPM1(4))))); | |
676 t = 3 - r1 * hfx; | |
677 var e = hxs * ((r1 - t) / (6 - x * t)); | |
678 if (k === 0) { // c is 0 | |
679 return x - (x*e - hxs); | |
680 } else { | |
681 e = (x * (e - c) - c); | |
682 e -= hxs; | |
683 if (k === -1) return 0.5 * (x - e) - 0.5; | |
684 if (k === 1) { | |
685 if (x < -0.25) return -2 * (e - (x + 0.5)); | |
686 return 1 + 2 * (x - e); | |
687 } | |
688 | |
689 if (k <= -2 || k > 56) { | |
690 // suffice to return exp(x) + 1 | |
691 y = 1 - (e - x); | |
692 // Add k to y's exponent | |
693 y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y)); | |
694 return y - 1; | |
695 } | |
696 if (k < 20) { | |
697 // t = 1 - 2^k | |
698 t = %_ConstructDouble(0x3ff00000 - (0x200000 >> k), 0); | |
699 y = t - (e - x); | |
700 // Add k to y's exponent | |
701 y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y)); | |
702 } else { | |
703 // t = 2^-k | |
704 t = %_ConstructDouble((0x3ff - k) << 20, 0); | |
705 y = x - (e + t); | |
706 y += 1; | |
707 // Add k to y's exponent | |
708 y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y)); | |
709 } | |
710 } | |
711 return y; | |
712 } | |
OLD | NEW |