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 |