| 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-2004 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 // Double constants that do not have empty lower 32 bits are found in fdlibm.cc | |
| 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 // rempio2result is used as a container for return values of %RemPiO2. It is | 19 // rempio2result is used as a container for return values of %RemPiO2. It is |
| 24 // initialized to a two-element Float64Array during genesis. | 20 // initialized to a two-element Float64Array during genesis. |
| 25 | 21 |
| 26 (function(global, utils) { | 22 (function(global, utils) { |
| 27 | 23 |
| 28 "use strict"; | 24 "use strict"; |
| 29 | 25 |
| 30 %CheckIsBootstrapping(); | 26 %CheckIsBootstrapping(); |
| 31 | 27 |
| 32 // ------------------------------------------------------------------- | 28 // ------------------------------------------------------------------- |
| 33 // Imports | 29 // Imports |
| 34 | 30 |
| 35 var GlobalMath = global.Math; | 31 var GlobalMath = global.Math; |
| 36 var kMath; | |
| 37 var MathAbs; | 32 var MathAbs; |
| 38 var MathExp; | 33 var MathExp; |
| 39 var NaN = %GetRootNaN(); | 34 var NaN = %GetRootNaN(); |
| 40 var rempio2result; | 35 var rempio2result; |
| 41 | 36 |
| 42 utils.Import(function(from) { | 37 utils.Import(function(from) { |
| 43 MathAbs = from.MathAbs; | 38 MathAbs = from.MathAbs; |
| 44 MathExp = from.MathExp; | 39 MathExp = from.MathExp; |
| 45 }); | 40 }); |
| 46 | 41 |
| 47 utils.SetupTypedArray(function(arg1, arg2, arg3) { | 42 utils.SetupTypedArray(function(arg1, arg2) { |
| 48 kMath = arg2; | 43 rempio2result = arg2; |
| 49 rempio2result = arg3; | |
| 50 }); | 44 }); |
| 51 | 45 |
| 52 // ------------------------------------------------------------------- | 46 // ------------------------------------------------------------------- |
| 53 | 47 |
| 54 define INVPIO2 = kMath[0]; | 48 define INVPIO2 = 6.36619772367581382433e-01; |
| 55 define PIO2_1 = kMath[1]; | 49 define PIO2_1 = 1.57079632673412561417; |
| 56 define PIO2_1T = kMath[2]; | 50 define PIO2_1T = 6.07710050650619224932e-11; |
| 57 define PIO2_2 = kMath[3]; | 51 define PIO2_2 = 6.07710050630396597660e-11; |
| 58 define PIO2_2T = kMath[4]; | 52 define PIO2_2T = 2.02226624879595063154e-21; |
| 59 define PIO2_3 = kMath[5]; | 53 define PIO2_3 = 2.02226624871116645580e-21; |
| 60 define PIO2_3T = kMath[6]; | 54 define PIO2_3T = 8.47842766036889956997e-32; |
| 61 define PIO4 = kMath[32]; | 55 define PIO4 = 7.85398163397448278999e-01; |
| 62 define PIO4LO = kMath[33]; | 56 define PIO4LO = 3.06161699786838301793e-17; |
| 63 | 57 |
| 64 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For | 58 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For |
| 65 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 | 59 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 |
| 66 // to more than double precision. | 60 // to more than double precision. |
| 67 | 61 |
| 68 macro REMPIO2(X) | 62 macro REMPIO2(X) |
| 69 var n, y0, y1; | 63 var n, y0, y1; |
| 70 var hx = %_DoubleHi(X); | 64 var hx = %_DoubleHi(X); |
| 71 var ix = hx & 0x7fffffff; | 65 var ix = hx & 0x7fffffff; |
| 72 | 66 |
| (...skipping 191 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 264 // 3 2 | 258 // 3 2 |
| 265 // tan(x+y) = x + (T1*x + (x *(r+y)+y)) | 259 // tan(x+y) = x + (T1*x + (x *(r+y)+y)) |
| 266 // | 260 // |
| 267 // 4. For x in [0.67434,pi/4], let y = pi/4 - x, then | 261 // 4. For x in [0.67434,pi/4], let y = pi/4 - x, then |
| 268 // tan(x) = ieee_tan(pi/4-y) = (1-ieee_tan(y))/(1+ieee_tan(y)) | 262 // tan(x) = ieee_tan(pi/4-y) = (1-ieee_tan(y))/(1+ieee_tan(y)) |
| 269 // = 1 - 2*(ieee_tan(y) - (ieee_tan(y)^2)/(1+ieee_tan(y))) | 263 // = 1 - 2*(ieee_tan(y) - (ieee_tan(y)^2)/(1+ieee_tan(y))) |
| 270 // | 264 // |
| 271 // Set returnTan to 1 for tan; -1 for cot. Anything else is illegal | 265 // Set returnTan to 1 for tan; -1 for cot. Anything else is illegal |
| 272 // and will cause incorrect results. | 266 // and will cause incorrect results. |
| 273 // | 267 // |
| 274 macro KTAN(x) | 268 define T00 = 3.33333333333334091986e-01; |
| 275 kMath[19+x] | 269 define T01 = 1.33333333333201242699e-01; |
| 276 endmacro | 270 define T02 = 5.39682539762260521377e-02; |
| 271 define T03 = 2.18694882948595424599e-02; |
| 272 define T04 = 8.86323982359930005737e-03; |
| 273 define T05 = 3.59207910759131235356e-03; |
| 274 define T06 = 1.45620945432529025516e-03; |
| 275 define T07 = 5.88041240820264096874e-04; |
| 276 define T08 = 2.46463134818469906812e-04; |
| 277 define T09 = 7.81794442939557092300e-05; |
| 278 define T10 = 7.14072491382608190305e-05; |
| 279 define T11 = -1.85586374855275456654e-05; |
| 280 define T12 = 2.59073051863633712884e-05; |
| 277 | 281 |
| 278 function KernelTan(x, y, returnTan) { | 282 function KernelTan(x, y, returnTan) { |
| 279 var z; | 283 var z; |
| 280 var w; | 284 var w; |
| 281 var hx = %_DoubleHi(x); | 285 var hx = %_DoubleHi(x); |
| 282 var ix = hx & 0x7fffffff; | 286 var ix = hx & 0x7fffffff; |
| 283 | 287 |
| 284 if (ix < 0x3e300000) { // |x| < 2^-28 | 288 if (ix < 0x3e300000) { // |x| < 2^-28 |
| 285 if (((ix | %_DoubleLo(x)) | (returnTan + 1)) == 0) { | 289 if (((ix | %_DoubleLo(x)) | (returnTan + 1)) == 0) { |
| 286 // x == 0 && returnTan = -1 | 290 // x == 0 && returnTan = -1 |
| (...skipping 22 matching lines...) Expand all Loading... |
| 309 w = PIO4LO - y; | 313 w = PIO4LO - y; |
| 310 x = z + w; | 314 x = z + w; |
| 311 y = 0; | 315 y = 0; |
| 312 } | 316 } |
| 313 z = x * x; | 317 z = x * x; |
| 314 w = z * z; | 318 w = z * z; |
| 315 | 319 |
| 316 // Break x^5 * (T1 + x^2*T2 + ...) into | 320 // Break x^5 * (T1 + x^2*T2 + ...) into |
| 317 // x^5 * (T1 + x^4*T3 + ... + x^20*T11) + | 321 // x^5 * (T1 + x^4*T3 + ... + x^20*T11) + |
| 318 // x^5 * (x^2 * (T2 + x^4*T4 + ... + x^22*T12)) | 322 // x^5 * (x^2 * (T2 + x^4*T4 + ... + x^22*T12)) |
| 319 var r = KTAN(1) + w * (KTAN(3) + w * (KTAN(5) + | 323 var r = T01 + w * (T03 + w * (T05 + |
| 320 w * (KTAN(7) + w * (KTAN(9) + w * KTAN(11))))); | 324 w * (T07 + w * (T09 + w * T11)))); |
| 321 var v = z * (KTAN(2) + w * (KTAN(4) + w * (KTAN(6) + | 325 var v = z * (T02 + w * (T04 + w * (T06 + |
| 322 w * (KTAN(8) + w * (KTAN(10) + w * KTAN(12)))))); | 326 w * (T08 + w * (T10 + w * T12))))); |
| 323 var s = z * x; | 327 var s = z * x; |
| 324 r = y + z * (s * (r + v) + y); | 328 r = y + z * (s * (r + v) + y); |
| 325 r = r + KTAN(0) * s; | 329 r = r + T00 * s; |
| 326 w = x + r; | 330 w = x + r; |
| 327 if (ix >= 0x3fe59428) { | 331 if (ix >= 0x3fe59428) { |
| 328 return (1 - ((hx >> 30) & 2)) * | 332 return (1 - ((hx >> 30) & 2)) * |
| 329 (returnTan - 2.0 * (x - (w * w / (w + returnTan) - r))); | 333 (returnTan - 2.0 * (x - (w * w / (w + returnTan) - r))); |
| 330 } | 334 } |
| 331 if (returnTan == 1) { | 335 if (returnTan == 1) { |
| 332 return w; | 336 return w; |
| 333 } else { | 337 } else { |
| 334 z = %_ConstructDouble(%_DoubleHi(w), 0); | 338 z = %_ConstructDouble(%_DoubleHi(w), 0); |
| 335 v = r - (z - x); | 339 v = r - (z - x); |
| (...skipping 112 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 448 // | 452 // |
| 449 // Note: Assuming log() return accurate answer, the following | 453 // Note: Assuming log() return accurate answer, the following |
| 450 // algorithm can be used to compute log1p(x) to within a few ULP: | 454 // algorithm can be used to compute log1p(x) to within a few ULP: |
| 451 // | 455 // |
| 452 // u = 1+x; | 456 // u = 1+x; |
| 453 // if (u==1.0) return x ; else | 457 // if (u==1.0) return x ; else |
| 454 // return log(u)*(x/(u-1.0)); | 458 // return log(u)*(x/(u-1.0)); |
| 455 // | 459 // |
| 456 // See HP-15C Advanced Functions Handbook, p.193. | 460 // See HP-15C Advanced Functions Handbook, p.193. |
| 457 // | 461 // |
| 458 define LN2_HI = kMath[34]; | 462 define LN2_HI = 6.93147180369123816490e-01; |
| 459 define LN2_LO = kMath[35]; | 463 define LN2_LO = 1.90821492927058770002e-10; |
| 460 define TWO_THIRD = kMath[36]; | 464 define TWO_THIRD = 6.666666666666666666e-01; |
| 461 macro KLOG1P(x) | 465 define LP1 = 6.666666666666735130e-01; |
| 462 (kMath[37+x]) | 466 define LP2 = 3.999999999940941908e-01; |
| 463 endmacro | 467 define LP3 = 2.857142874366239149e-01; |
| 468 define LP4 = 2.222219843214978396e-01; |
| 469 define LP5 = 1.818357216161805012e-01; |
| 470 define LP6 = 1.531383769920937332e-01; |
| 471 define LP7 = 1.479819860511658591e-01; |
| 472 |
| 464 // 2^54 | 473 // 2^54 |
| 465 define TWO54 = 18014398509481984; | 474 define TWO54 = 18014398509481984; |
| 466 | 475 |
| 467 function MathLog1p(x) { | 476 function MathLog1p(x) { |
| 468 x = x * 1; // Convert to number. | 477 x = x * 1; // Convert to number. |
| 469 var hx = %_DoubleHi(x); | 478 var hx = %_DoubleHi(x); |
| 470 var ax = hx & 0x7fffffff; | 479 var ax = hx & 0x7fffffff; |
| 471 var k = 1; | 480 var k = 1; |
| 472 var f = x; | 481 var f = x; |
| 473 var hu = 1; | 482 var hu = 1; |
| (...skipping 61 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 535 var R = hfsq * (1 - TWO_THIRD * f); | 544 var R = hfsq * (1 - TWO_THIRD * f); |
| 536 if (k === 0) { | 545 if (k === 0) { |
| 537 return f - R; | 546 return f - R; |
| 538 } else { | 547 } else { |
| 539 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); | 548 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); |
| 540 } | 549 } |
| 541 } | 550 } |
| 542 | 551 |
| 543 var s = f / (2 + f); | 552 var s = f / (2 + f); |
| 544 var z = s * s; | 553 var z = s * s; |
| 545 var R = z * (KLOG1P(0) + z * (KLOG1P(1) + z * | 554 var R = z * (LP1 + z * (LP2 + z * (LP3 + z * (LP4 + |
| 546 (KLOG1P(2) + z * (KLOG1P(3) + z * | 555 z * (LP5 + z * (LP6 + z * LP7)))))); |
| 547 (KLOG1P(4) + z * (KLOG1P(5) + z * KLOG1P(6))))))); | |
| 548 if (k === 0) { | 556 if (k === 0) { |
| 549 return f - (hfsq - s * (hfsq + R)); | 557 return f - (hfsq - s * (hfsq + R)); |
| 550 } else { | 558 } else { |
| 551 return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); | 559 return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); |
| 552 } | 560 } |
| 553 } | 561 } |
| 554 | 562 |
| 555 // ES6 draft 09-27-13, section 20.2.2.14. | 563 // ES6 draft 09-27-13, section 20.2.2.14. |
| 556 // Math.expm1 | 564 // Math.expm1 |
| 557 // Returns exp(x)-1, the exponential of x minus 1. | 565 // Returns exp(x)-1, the exponential of x minus 1. |
| (...skipping 76 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 634 // for finite argument, only expm1(0)=0 is exact. | 642 // for finite argument, only expm1(0)=0 is exact. |
| 635 // | 643 // |
| 636 // Accuracy: | 644 // Accuracy: |
| 637 // according to an error analysis, the error is always less than | 645 // according to an error analysis, the error is always less than |
| 638 // 1 ulp (unit in the last place). | 646 // 1 ulp (unit in the last place). |
| 639 // | 647 // |
| 640 // Misc. info. | 648 // Misc. info. |
| 641 // For IEEE double | 649 // For IEEE double |
| 642 // if x > 7.09782712893383973096e+02 then expm1(x) overflow | 650 // if x > 7.09782712893383973096e+02 then expm1(x) overflow |
| 643 // | 651 // |
| 644 define KEXPM1_OVERFLOW = kMath[44]; | 652 define KEXPM1_OVERFLOW = 7.09782712893383973096e+02; |
| 645 define INVLN2 = kMath[45]; | 653 define INVLN2 = 1.44269504088896338700; |
| 646 macro KEXPM1(x) | 654 define EXPM1_1 = -3.33333333333331316428e-02; |
| 647 (kMath[46+x]) | 655 define EXPM1_2 = 1.58730158725481460165e-03; |
| 648 endmacro | 656 define EXPM1_3 = -7.93650757867487942473e-05; |
| 657 define EXPM1_4 = 4.00821782732936239552e-06; |
| 658 define EXPM1_5 = -2.01099218183624371326e-07; |
| 649 | 659 |
| 650 function MathExpm1(x) { | 660 function MathExpm1(x) { |
| 651 x = x * 1; // Convert to number. | 661 x = x * 1; // Convert to number. |
| 652 var y; | 662 var y; |
| 653 var hi; | 663 var hi; |
| 654 var lo; | 664 var lo; |
| 655 var k; | 665 var k; |
| 656 var t; | 666 var t; |
| 657 var c; | 667 var c; |
| 658 | 668 |
| (...skipping 39 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 698 // When |x| < 2^-54, we can return x. | 708 // When |x| < 2^-54, we can return x. |
| 699 return x; | 709 return x; |
| 700 } else { | 710 } else { |
| 701 // Fall through. | 711 // Fall through. |
| 702 k = 0; | 712 k = 0; |
| 703 } | 713 } |
| 704 | 714 |
| 705 // x is now in primary range | 715 // x is now in primary range |
| 706 var hfx = 0.5 * x; | 716 var hfx = 0.5 * x; |
| 707 var hxs = x * hfx; | 717 var hxs = x * hfx; |
| 708 var r1 = 1 + hxs * (KEXPM1(0) + hxs * (KEXPM1(1) + hxs * | 718 var r1 = 1 + hxs * (EXPM1_1 + hxs * (EXPM1_2 + hxs * |
| 709 (KEXPM1(2) + hxs * (KEXPM1(3) + hxs * KEXPM1(4))))); | 719 (EXPM1_3 + hxs * (EXPM1_4 + hxs * EXPM1_5)))); |
| 710 t = 3 - r1 * hfx; | 720 t = 3 - r1 * hfx; |
| 711 var e = hxs * ((r1 - t) / (6 - x * t)); | 721 var e = hxs * ((r1 - t) / (6 - x * t)); |
| 712 if (k === 0) { // c is 0 | 722 if (k === 0) { // c is 0 |
| 713 return x - (x*e - hxs); | 723 return x - (x*e - hxs); |
| 714 } else { | 724 } else { |
| 715 e = (x * (e - c) - c); | 725 e = (x * (e - c) - c); |
| 716 e -= hxs; | 726 e -= hxs; |
| 717 if (k === -1) return 0.5 * (x - e) - 0.5; | 727 if (k === -1) return 0.5 * (x - e) - 0.5; |
| 718 if (k === 1) { | 728 if (k === 1) { |
| 719 if (x < -0.25) return -2 * (e - (x + 0.5)); | 729 if (x < -0.25) return -2 * (e - (x + 0.5)); |
| (...skipping 37 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 757 // 2 | 767 // 2 |
| 758 // | 768 // |
| 759 // 22 <= x <= lnovft : sinh(x) := exp(x)/2 | 769 // 22 <= x <= lnovft : sinh(x) := exp(x)/2 |
| 760 // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) | 770 // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) |
| 761 // ln2ovft < x : sinh(x) := x*shuge (overflow) | 771 // ln2ovft < x : sinh(x) := x*shuge (overflow) |
| 762 // | 772 // |
| 763 // Special cases: | 773 // Special cases: |
| 764 // sinh(x) is |x| if x is +Infinity, -Infinity, or NaN. | 774 // sinh(x) is |x| if x is +Infinity, -Infinity, or NaN. |
| 765 // only sinh(0)=0 is exact for finite x. | 775 // only sinh(0)=0 is exact for finite x. |
| 766 // | 776 // |
| 767 define KSINH_OVERFLOW = kMath[51]; | 777 define KSINH_OVERFLOW = 710.4758600739439; |
| 768 define TWO_M28 = 3.725290298461914e-9; // 2^-28, empty lower half | 778 define TWO_M28 = 3.725290298461914e-9; // 2^-28, empty lower half |
| 769 define LOG_MAXD = 709.7822265625; // 0x40862e42 00000000, empty lower half | 779 define LOG_MAXD = 709.7822265625; // 0x40862e42 00000000, empty lower half |
| 770 | 780 |
| 771 function MathSinh(x) { | 781 function MathSinh(x) { |
| 772 x = x * 1; // Convert to number. | 782 x = x * 1; // Convert to number. |
| 773 var h = (x < 0) ? -0.5 : 0.5; | 783 var h = (x < 0) ? -0.5 : 0.5; |
| 774 // |x| in [0, 22]. return sign(x)*0.5*(E+E/(E+1)) | 784 // |x| in [0, 22]. return sign(x)*0.5*(E+E/(E+1)) |
| 775 var ax = MathAbs(x); | 785 var ax = MathAbs(x); |
| 776 if (ax < 22) { | 786 if (ax < 22) { |
| 777 // For |x| < 2^-28, sinh(x) = x | 787 // For |x| < 2^-28, sinh(x) = x |
| (...skipping 31 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 809 // ln2/2 <= x <= 22 : cosh(x) := ------------------- | 819 // ln2/2 <= x <= 22 : cosh(x) := ------------------- |
| 810 // 2 | 820 // 2 |
| 811 // 22 <= x <= lnovft : cosh(x) := exp(x)/2 | 821 // 22 <= x <= lnovft : cosh(x) := exp(x)/2 |
| 812 // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) | 822 // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) |
| 813 // ln2ovft < x : cosh(x) := huge*huge (overflow) | 823 // ln2ovft < x : cosh(x) := huge*huge (overflow) |
| 814 // | 824 // |
| 815 // Special cases: | 825 // Special cases: |
| 816 // cosh(x) is |x| if x is +INF, -INF, or NaN. | 826 // cosh(x) is |x| if x is +INF, -INF, or NaN. |
| 817 // only cosh(0)=1 is exact for finite x. | 827 // only cosh(0)=1 is exact for finite x. |
| 818 // | 828 // |
| 819 define KCOSH_OVERFLOW = kMath[51]; | 829 define KCOSH_OVERFLOW = 710.4758600739439; |
| 820 | 830 |
| 821 function MathCosh(x) { | 831 function MathCosh(x) { |
| 822 x = x * 1; // Convert to number. | 832 x = x * 1; // Convert to number. |
| 823 var ix = %_DoubleHi(x) & 0x7fffffff; | 833 var ix = %_DoubleHi(x) & 0x7fffffff; |
| 824 // |x| in [0,0.5*log2], return 1+expm1(|x|)^2/(2*exp(|x|)) | 834 // |x| in [0,0.5*log2], return 1+expm1(|x|)^2/(2*exp(|x|)) |
| 825 if (ix < 0x3fd62e43) { | 835 if (ix < 0x3fd62e43) { |
| 826 var t = MathExpm1(MathAbs(x)); | 836 var t = MathExpm1(MathAbs(x)); |
| 827 var w = 1 + t; | 837 var w = 1 + t; |
| 828 // For |x| < 2^-55, cosh(x) = 1 | 838 // For |x| < 2^-55, cosh(x) = 1 |
| 829 if (ix < 0x3c800000) return w; | 839 if (ix < 0x3c800000) return w; |
| (...skipping 94 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 924 // [1/log(10)] rounded to 53 bits has error .198 ulps; | 934 // [1/log(10)] rounded to 53 bits has error .198 ulps; |
| 925 // log10 is monotonic at all binary break points. | 935 // log10 is monotonic at all binary break points. |
| 926 // | 936 // |
| 927 // Special cases: | 937 // Special cases: |
| 928 // log10(x) is NaN if x < 0; | 938 // log10(x) is NaN if x < 0; |
| 929 // log10(+INF) is +INF; log10(0) is -INF; | 939 // log10(+INF) is +INF; log10(0) is -INF; |
| 930 // log10(NaN) is that NaN; | 940 // log10(NaN) is that NaN; |
| 931 // log10(10**N) = N for N=0,1,...,22. | 941 // log10(10**N) = N for N=0,1,...,22. |
| 932 // | 942 // |
| 933 | 943 |
| 934 define IVLN10 = kMath[52]; | 944 define IVLN10 = 4.34294481903251816668e-01; |
| 935 define LOG10_2HI = kMath[53]; | 945 define LOG10_2HI = 3.01029995663611771306e-01; |
| 936 define LOG10_2LO = kMath[54]; | 946 define LOG10_2LO = 3.69423907715893078616e-13; |
| 937 | 947 |
| 938 function MathLog10(x) { | 948 function MathLog10(x) { |
| 939 x = x * 1; // Convert to number. | 949 x = x * 1; // Convert to number. |
| 940 var hx = %_DoubleHi(x); | 950 var hx = %_DoubleHi(x); |
| 941 var lx = %_DoubleLo(x); | 951 var lx = %_DoubleLo(x); |
| 942 var k = 0; | 952 var k = 0; |
| 943 | 953 |
| 944 if (hx < 0x00100000) { | 954 if (hx < 0x00100000) { |
| 945 // x < 2^-1022 | 955 // x < 2^-1022 |
| 946 // log10(+/- 0) = -Infinity. | 956 // log10(+/- 0) = -Infinity. |
| (...skipping 27 matching lines...) Expand all Loading... |
| 974 // fdlibm does not have an explicit log2 function, but fdlibm's pow | 984 // fdlibm does not have an explicit log2 function, but fdlibm's pow |
| 975 // function does implement an accurate log2 function as part of the | 985 // function does implement an accurate log2 function as part of the |
| 976 // pow implementation. This extracts the core parts of that as a | 986 // pow implementation. This extracts the core parts of that as a |
| 977 // separate log2 function. | 987 // separate log2 function. |
| 978 | 988 |
| 979 // Method: | 989 // Method: |
| 980 // Compute log2(x) in two pieces: | 990 // Compute log2(x) in two pieces: |
| 981 // log2(x) = w1 + w2 | 991 // log2(x) = w1 + w2 |
| 982 // where w1 has 53-24 = 29 bits of trailing zeroes. | 992 // where w1 has 53-24 = 29 bits of trailing zeroes. |
| 983 | 993 |
| 984 define DP_H = kMath[64]; | 994 define DP_H = 5.84962487220764160156e-01; |
| 985 define DP_L = kMath[65]; | 995 define DP_L = 1.35003920212974897128e-08; |
| 986 | 996 |
| 987 // Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3) | 997 // Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3) |
| 988 macro KLOG2(x) | 998 define LOG2_1 = 5.99999999999994648725e-01; |
| 989 (kMath[55+x]) | 999 define LOG2_2 = 4.28571428578550184252e-01; |
| 990 endmacro | 1000 define LOG2_3 = 3.33333329818377432918e-01; |
| 1001 define LOG2_4 = 2.72728123808534006489e-01; |
| 1002 define LOG2_5 = 2.30660745775561754067e-01; |
| 1003 define LOG2_6 = 2.06975017800338417784e-01; |
| 991 | 1004 |
| 992 // cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy. | 1005 // cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy. |
| 993 define CP = kMath[61]; | 1006 define CP = 9.61796693925975554329e-01; |
| 994 define CP_H = kMath[62]; | 1007 define CP_H = 9.61796700954437255859e-01; |
| 995 define CP_L = kMath[63]; | 1008 define CP_L = -7.02846165095275826516e-09; |
| 996 // 2^53 | 1009 // 2^53 |
| 997 define TWO53 = 9007199254740992; | 1010 define TWO53 = 9007199254740992; |
| 998 | 1011 |
| 999 function MathLog2(x) { | 1012 function MathLog2(x) { |
| 1000 x = x * 1; // Convert to number. | 1013 x = x * 1; // Convert to number. |
| 1001 var ax = MathAbs(x); | 1014 var ax = MathAbs(x); |
| 1002 var hx = %_DoubleHi(x); | 1015 var hx = %_DoubleHi(x); |
| 1003 var lx = %_DoubleLo(x); | 1016 var lx = %_DoubleLo(x); |
| 1004 var ix = hx & 0x7fffffff; | 1017 var ix = hx & 0x7fffffff; |
| 1005 | 1018 |
| (...skipping 44 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 1050 var ss = u * v; | 1063 var ss = u * v; |
| 1051 var s_h = %_ConstructDouble(%_DoubleHi(ss), 0); | 1064 var s_h = %_ConstructDouble(%_DoubleHi(ss), 0); |
| 1052 | 1065 |
| 1053 // t_h = ax + bp[k] High | 1066 // t_h = ax + bp[k] High |
| 1054 var t_h = %_ConstructDouble(%_DoubleHi(ax + bp), 0) | 1067 var t_h = %_ConstructDouble(%_DoubleHi(ax + bp), 0) |
| 1055 var t_l = ax - (t_h - bp); | 1068 var t_l = ax - (t_h - bp); |
| 1056 var s_l = v * ((u - s_h * t_h) - s_h * t_l); | 1069 var s_l = v * ((u - s_h * t_h) - s_h * t_l); |
| 1057 | 1070 |
| 1058 // Compute log2(ax) | 1071 // Compute log2(ax) |
| 1059 var s2 = ss * ss; | 1072 var s2 = ss * ss; |
| 1060 var r = s2 * s2 * (KLOG2(0) + s2 * (KLOG2(1) + s2 * (KLOG2(2) + s2 * ( | 1073 var r = s2 * s2 * (LOG2_1 + s2 * (LOG2_2 + s2 * (LOG2_3 + s2 * ( |
| 1061 KLOG2(3) + s2 * (KLOG2(4) + s2 * KLOG2(5)))))); | 1074 LOG2_4 + s2 * (LOG2_5 + s2 * LOG2_6))))); |
| 1062 r += s_l * (s_h + ss); | 1075 r += s_l * (s_h + ss); |
| 1063 s2 = s_h * s_h; | 1076 s2 = s_h * s_h; |
| 1064 t_h = %_ConstructDouble(%_DoubleHi(3.0 + s2 + r), 0); | 1077 t_h = %_ConstructDouble(%_DoubleHi(3.0 + s2 + r), 0); |
| 1065 t_l = r - ((t_h - 3.0) - s2); | 1078 t_l = r - ((t_h - 3.0) - s2); |
| 1066 // u + v = ss * (1 + ...) | 1079 // u + v = ss * (1 + ...) |
| 1067 u = s_h * t_h; | 1080 u = s_h * t_h; |
| 1068 v = s_l * t_h + t_l * ss; | 1081 v = s_l * t_h + t_l * ss; |
| 1069 | 1082 |
| 1070 // 2 / (3 * log(2)) * (ss + ...) | 1083 // 2 / (3 * log(2)) * (ss + ...) |
| 1071 p_h = %_ConstructDouble(%_DoubleHi(u + v), 0); | 1084 p_h = %_ConstructDouble(%_DoubleHi(u + v), 0); |
| (...skipping 22 matching lines...) Expand all Loading... |
| 1094 "log10", MathLog10, | 1107 "log10", MathLog10, |
| 1095 "log2", MathLog2, | 1108 "log2", MathLog2, |
| 1096 "log1p", MathLog1p, | 1109 "log1p", MathLog1p, |
| 1097 "expm1", MathExpm1 | 1110 "expm1", MathExpm1 |
| 1098 ]); | 1111 ]); |
| 1099 | 1112 |
| 1100 %SetForceInlineFlag(MathSin); | 1113 %SetForceInlineFlag(MathSin); |
| 1101 %SetForceInlineFlag(MathCos); | 1114 %SetForceInlineFlag(MathCos); |
| 1102 | 1115 |
| 1103 }) | 1116 }) |
| OLD | NEW |