| 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 // ==================================================== |
| (...skipping 18 matching lines...) Expand all Loading... |
| 29 (function() { | 29 (function() { |
| 30 | 30 |
| 31 "use strict"; | 31 "use strict"; |
| 32 | 32 |
| 33 %CheckIsBootstrapping(); | 33 %CheckIsBootstrapping(); |
| 34 | 34 |
| 35 var GlobalMath = global.Math; | 35 var GlobalMath = global.Math; |
| 36 | 36 |
| 37 //------------------------------------------------------------------- | 37 //------------------------------------------------------------------- |
| 38 | 38 |
| 39 const INVPIO2 = kMath[0]; | 39 define INVPIO2 = kMath[0]; |
| 40 const PIO2_1 = kMath[1]; | 40 define PIO2_1 = kMath[1]; |
| 41 const PIO2_1T = kMath[2]; | 41 define PIO2_1T = kMath[2]; |
| 42 const PIO2_2 = kMath[3]; | 42 define PIO2_2 = kMath[3]; |
| 43 const PIO2_2T = kMath[4]; | 43 define PIO2_2T = kMath[4]; |
| 44 const PIO2_3 = kMath[5]; | 44 define PIO2_3 = kMath[5]; |
| 45 const PIO2_3T = kMath[6]; | 45 define PIO2_3T = kMath[6]; |
| 46 const PIO4 = kMath[32]; | 46 define PIO4 = kMath[32]; |
| 47 const PIO4LO = kMath[33]; | 47 define PIO4LO = kMath[33]; |
| 48 | 48 |
| 49 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For | 49 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For |
| 50 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 | 50 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 |
| 51 // to more than double precision. | 51 // to more than double precision. |
| 52 | 52 |
| 53 macro REMPIO2(X) | 53 macro REMPIO2(X) |
| 54 var n, y0, y1; | 54 var n, y0, y1; |
| 55 var hx = %_DoubleHi(X); | 55 var hx = %_DoubleHi(X); |
| 56 var ix = hx & 0x7fffffff; | 56 var ix = hx & 0x7fffffff; |
| 57 | 57 |
| (...skipping 84 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 142 // | x | | 142 // | x | |
| 143 // | 143 // |
| 144 // 3. ieee_sin(X+Y) = ieee_sin(X) + sin'(X')*Y | 144 // 3. ieee_sin(X+Y) = ieee_sin(X) + sin'(X')*Y |
| 145 // ~ ieee_sin(X) + (1-X*X/2)*Y | 145 // ~ ieee_sin(X) + (1-X*X/2)*Y |
| 146 // For better accuracy, let | 146 // For better accuracy, let |
| 147 // 3 2 2 2 2 | 147 // 3 2 2 2 2 |
| 148 // r = X *(S2+X *(S3+X *(S4+X *(S5+X *S6)))) | 148 // r = X *(S2+X *(S3+X *(S4+X *(S5+X *S6)))) |
| 149 // then 3 2 | 149 // then 3 2 |
| 150 // sin(x) = X + (S1*X + (X *(r-Y/2)+Y)) | 150 // sin(x) = X + (S1*X + (X *(r-Y/2)+Y)) |
| 151 // | 151 // |
| 152 const S1 = -1.66666666666666324348e-01; | 152 define S1 = -1.66666666666666324348e-01; |
| 153 const S2 = 8.33333333332248946124e-03; | 153 define S2 = 8.33333333332248946124e-03; |
| 154 const S3 = -1.98412698298579493134e-04; | 154 define S3 = -1.98412698298579493134e-04; |
| 155 const S4 = 2.75573137070700676789e-06; | 155 define S4 = 2.75573137070700676789e-06; |
| 156 const S5 = -2.50507602534068634195e-08; | 156 define S5 = -2.50507602534068634195e-08; |
| 157 const S6 = 1.58969099521155010221e-10; | 157 define S6 = 1.58969099521155010221e-10; |
| 158 | 158 |
| 159 macro RETURN_KERNELSIN(X, Y, SIGN) | 159 macro RETURN_KERNELSIN(X, Y, SIGN) |
| 160 var z = X * X; | 160 var z = X * X; |
| 161 var v = z * X; | 161 var v = z * X; |
| 162 var r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); | 162 var r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); |
| 163 return (X - ((z * (0.5 * Y - v * r) - Y) - v * S1)) SIGN; | 163 return (X - ((z * (0.5 * Y - v * r) - Y) - v * S1)) SIGN; |
| 164 endmacro | 164 endmacro |
| 165 | 165 |
| 166 // __kernel_cos(X, Y) | 166 // __kernel_cos(X, Y) |
| 167 // kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 | 167 // kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 |
| (...skipping 20 matching lines...) Expand all Loading... |
| 188 // a correction term is necessary in ieee_cos(x) and hence | 188 // a correction term is necessary in ieee_cos(x) and hence |
| 189 // cos(X+Y) = 1 - (X*X/2 - (r - X*Y)) | 189 // cos(X+Y) = 1 - (X*X/2 - (r - X*Y)) |
| 190 // For better accuracy when x > 0.3, let qx = |x|/4 with | 190 // For better accuracy when x > 0.3, let qx = |x|/4 with |
| 191 // the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. | 191 // the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. |
| 192 // Then | 192 // Then |
| 193 // cos(X+Y) = (1-qx) - ((X*X/2-qx) - (r-X*Y)). | 193 // cos(X+Y) = (1-qx) - ((X*X/2-qx) - (r-X*Y)). |
| 194 // Note that 1-qx and (X*X/2-qx) is EXACT here, and the | 194 // Note that 1-qx and (X*X/2-qx) is EXACT here, and the |
| 195 // magnitude of the latter is at least a quarter of X*X/2, | 195 // magnitude of the latter is at least a quarter of X*X/2, |
| 196 // thus, reducing the rounding error in the subtraction. | 196 // thus, reducing the rounding error in the subtraction. |
| 197 // | 197 // |
| 198 const C1 = 4.16666666666666019037e-02; | 198 define C1 = 4.16666666666666019037e-02; |
| 199 const C2 = -1.38888888888741095749e-03; | 199 define C2 = -1.38888888888741095749e-03; |
| 200 const C3 = 2.48015872894767294178e-05; | 200 define C3 = 2.48015872894767294178e-05; |
| 201 const C4 = -2.75573143513906633035e-07; | 201 define C4 = -2.75573143513906633035e-07; |
| 202 const C5 = 2.08757232129817482790e-09; | 202 define C5 = 2.08757232129817482790e-09; |
| 203 const C6 = -1.13596475577881948265e-11; | 203 define C6 = -1.13596475577881948265e-11; |
| 204 | 204 |
| 205 macro RETURN_KERNELCOS(X, Y, SIGN) | 205 macro RETURN_KERNELCOS(X, Y, SIGN) |
| 206 var ix = %_DoubleHi(X) & 0x7fffffff; | 206 var ix = %_DoubleHi(X) & 0x7fffffff; |
| 207 var z = X * X; | 207 var z = X * X; |
| 208 var r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); | 208 var r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); |
| 209 if (ix < 0x3fd33333) { // |x| ~< 0.3 | 209 if (ix < 0x3fd33333) { // |x| ~< 0.3 |
| 210 return (1 - (0.5 * z - (z * r - X * Y))) SIGN; | 210 return (1 - (0.5 * z - (z * r - X * Y))) SIGN; |
| 211 } else { | 211 } else { |
| 212 var qx; | 212 var qx; |
| 213 if (ix > 0x3fe90000) { // |x| > 0.78125 | 213 if (ix > 0x3fe90000) { // |x| > 0.78125 |
| (...skipping 219 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 433 // | 433 // |
| 434 // Note: Assuming log() return accurate answer, the following | 434 // Note: Assuming log() return accurate answer, the following |
| 435 // algorithm can be used to compute log1p(x) to within a few ULP: | 435 // algorithm can be used to compute log1p(x) to within a few ULP: |
| 436 // | 436 // |
| 437 // u = 1+x; | 437 // u = 1+x; |
| 438 // if (u==1.0) return x ; else | 438 // if (u==1.0) return x ; else |
| 439 // return log(u)*(x/(u-1.0)); | 439 // return log(u)*(x/(u-1.0)); |
| 440 // | 440 // |
| 441 // See HP-15C Advanced Functions Handbook, p.193. | 441 // See HP-15C Advanced Functions Handbook, p.193. |
| 442 // | 442 // |
| 443 const LN2_HI = kMath[34]; | 443 define LN2_HI = kMath[34]; |
| 444 const LN2_LO = kMath[35]; | 444 define LN2_LO = kMath[35]; |
| 445 const TWO_THIRD = kMath[36]; | 445 define TWO_THIRD = kMath[36]; |
| 446 macro KLOG1P(x) | 446 macro KLOG1P(x) |
| 447 (kMath[37+x]) | 447 (kMath[37+x]) |
| 448 endmacro | 448 endmacro |
| 449 // 2^54 | 449 // 2^54 |
| 450 const TWO54 = 18014398509481984; | 450 define TWO54 = 18014398509481984; |
| 451 | 451 |
| 452 function MathLog1p(x) { | 452 function MathLog1p(x) { |
| 453 x = x * 1; // Convert to number. | 453 x = x * 1; // Convert to number. |
| 454 var hx = %_DoubleHi(x); | 454 var hx = %_DoubleHi(x); |
| 455 var ax = hx & 0x7fffffff; | 455 var ax = hx & 0x7fffffff; |
| 456 var k = 1; | 456 var k = 1; |
| 457 var f = x; | 457 var f = x; |
| 458 var hu = 1; | 458 var hu = 1; |
| 459 var c = 0; | 459 var c = 0; |
| 460 var u = x; | 460 var u = x; |
| (...skipping 158 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 619 // for finite argument, only expm1(0)=0 is exact. | 619 // for finite argument, only expm1(0)=0 is exact. |
| 620 // | 620 // |
| 621 // Accuracy: | 621 // Accuracy: |
| 622 // according to an error analysis, the error is always less than | 622 // according to an error analysis, the error is always less than |
| 623 // 1 ulp (unit in the last place). | 623 // 1 ulp (unit in the last place). |
| 624 // | 624 // |
| 625 // Misc. info. | 625 // Misc. info. |
| 626 // For IEEE double | 626 // For IEEE double |
| 627 // if x > 7.09782712893383973096e+02 then expm1(x) overflow | 627 // if x > 7.09782712893383973096e+02 then expm1(x) overflow |
| 628 // | 628 // |
| 629 const KEXPM1_OVERFLOW = kMath[44]; | 629 define KEXPM1_OVERFLOW = kMath[44]; |
| 630 const INVLN2 = kMath[45]; | 630 define INVLN2 = kMath[45]; |
| 631 macro KEXPM1(x) | 631 macro KEXPM1(x) |
| 632 (kMath[46+x]) | 632 (kMath[46+x]) |
| 633 endmacro | 633 endmacro |
| 634 | 634 |
| 635 function MathExpm1(x) { | 635 function MathExpm1(x) { |
| 636 x = x * 1; // Convert to number. | 636 x = x * 1; // Convert to number. |
| 637 var y; | 637 var y; |
| 638 var hi; | 638 var hi; |
| 639 var lo; | 639 var lo; |
| 640 var k; | 640 var k; |
| (...skipping 101 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 742 // 2 | 742 // 2 |
| 743 // | 743 // |
| 744 // 22 <= x <= lnovft : sinh(x) := exp(x)/2 | 744 // 22 <= x <= lnovft : sinh(x) := exp(x)/2 |
| 745 // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) | 745 // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) |
| 746 // ln2ovft < x : sinh(x) := x*shuge (overflow) | 746 // ln2ovft < x : sinh(x) := x*shuge (overflow) |
| 747 // | 747 // |
| 748 // Special cases: | 748 // Special cases: |
| 749 // sinh(x) is |x| if x is +Infinity, -Infinity, or NaN. | 749 // sinh(x) is |x| if x is +Infinity, -Infinity, or NaN. |
| 750 // only sinh(0)=0 is exact for finite x. | 750 // only sinh(0)=0 is exact for finite x. |
| 751 // | 751 // |
| 752 const KSINH_OVERFLOW = kMath[51]; | 752 define KSINH_OVERFLOW = kMath[51]; |
| 753 const TWO_M28 = 3.725290298461914e-9; // 2^-28, empty lower half | 753 define TWO_M28 = 3.725290298461914e-9; // 2^-28, empty lower half |
| 754 const LOG_MAXD = 709.7822265625; // 0x40862e42 00000000, empty lower half | 754 define LOG_MAXD = 709.7822265625; // 0x40862e42 00000000, empty lower half |
| 755 | 755 |
| 756 function MathSinh(x) { | 756 function MathSinh(x) { |
| 757 x = x * 1; // Convert to number. | 757 x = x * 1; // Convert to number. |
| 758 var h = (x < 0) ? -0.5 : 0.5; | 758 var h = (x < 0) ? -0.5 : 0.5; |
| 759 // |x| in [0, 22]. return sign(x)*0.5*(E+E/(E+1)) | 759 // |x| in [0, 22]. return sign(x)*0.5*(E+E/(E+1)) |
| 760 var ax = $abs(x); | 760 var ax = $abs(x); |
| 761 if (ax < 22) { | 761 if (ax < 22) { |
| 762 // For |x| < 2^-28, sinh(x) = x | 762 // For |x| < 2^-28, sinh(x) = x |
| 763 if (ax < TWO_M28) return x; | 763 if (ax < TWO_M28) return x; |
| 764 var t = MathExpm1(ax); | 764 var t = MathExpm1(ax); |
| (...skipping 29 matching lines...) Expand all Loading... |
| 794 // ln2/2 <= x <= 22 : cosh(x) := ------------------- | 794 // ln2/2 <= x <= 22 : cosh(x) := ------------------- |
| 795 // 2 | 795 // 2 |
| 796 // 22 <= x <= lnovft : cosh(x) := exp(x)/2 | 796 // 22 <= x <= lnovft : cosh(x) := exp(x)/2 |
| 797 // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) | 797 // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) |
| 798 // ln2ovft < x : cosh(x) := huge*huge (overflow) | 798 // ln2ovft < x : cosh(x) := huge*huge (overflow) |
| 799 // | 799 // |
| 800 // Special cases: | 800 // Special cases: |
| 801 // cosh(x) is |x| if x is +INF, -INF, or NaN. | 801 // cosh(x) is |x| if x is +INF, -INF, or NaN. |
| 802 // only cosh(0)=1 is exact for finite x. | 802 // only cosh(0)=1 is exact for finite x. |
| 803 // | 803 // |
| 804 const KCOSH_OVERFLOW = kMath[51]; | 804 define KCOSH_OVERFLOW = kMath[51]; |
| 805 | 805 |
| 806 function MathCosh(x) { | 806 function MathCosh(x) { |
| 807 x = x * 1; // Convert to number. | 807 x = x * 1; // Convert to number. |
| 808 var ix = %_DoubleHi(x) & 0x7fffffff; | 808 var ix = %_DoubleHi(x) & 0x7fffffff; |
| 809 // |x| in [0,0.5*log2], return 1+expm1(|x|)^2/(2*exp(|x|)) | 809 // |x| in [0,0.5*log2], return 1+expm1(|x|)^2/(2*exp(|x|)) |
| 810 if (ix < 0x3fd62e43) { | 810 if (ix < 0x3fd62e43) { |
| 811 var t = MathExpm1($abs(x)); | 811 var t = MathExpm1($abs(x)); |
| 812 var w = 1 + t; | 812 var w = 1 + t; |
| 813 // For |x| < 2^-55, cosh(x) = 1 | 813 // For |x| < 2^-55, cosh(x) = 1 |
| 814 if (ix < 0x3c800000) return w; | 814 if (ix < 0x3c800000) return w; |
| (...skipping 37 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 852 // [1/log(10)] rounded to 53 bits has error .198 ulps; | 852 // [1/log(10)] rounded to 53 bits has error .198 ulps; |
| 853 // log10 is monotonic at all binary break points. | 853 // log10 is monotonic at all binary break points. |
| 854 // | 854 // |
| 855 // Special cases: | 855 // Special cases: |
| 856 // log10(x) is NaN if x < 0; | 856 // log10(x) is NaN if x < 0; |
| 857 // log10(+INF) is +INF; log10(0) is -INF; | 857 // log10(+INF) is +INF; log10(0) is -INF; |
| 858 // log10(NaN) is that NaN; | 858 // log10(NaN) is that NaN; |
| 859 // log10(10**N) = N for N=0,1,...,22. | 859 // log10(10**N) = N for N=0,1,...,22. |
| 860 // | 860 // |
| 861 | 861 |
| 862 const IVLN10 = kMath[52]; | 862 define IVLN10 = kMath[52]; |
| 863 const LOG10_2HI = kMath[53]; | 863 define LOG10_2HI = kMath[53]; |
| 864 const LOG10_2LO = kMath[54]; | 864 define LOG10_2LO = kMath[54]; |
| 865 | 865 |
| 866 function MathLog10(x) { | 866 function MathLog10(x) { |
| 867 x = x * 1; // Convert to number. | 867 x = x * 1; // Convert to number. |
| 868 var hx = %_DoubleHi(x); | 868 var hx = %_DoubleHi(x); |
| 869 var lx = %_DoubleLo(x); | 869 var lx = %_DoubleLo(x); |
| 870 var k = 0; | 870 var k = 0; |
| 871 | 871 |
| 872 if (hx < 0x00100000) { | 872 if (hx < 0x00100000) { |
| 873 // x < 2^-1022 | 873 // x < 2^-1022 |
| 874 // log10(+/- 0) = -Infinity. | 874 // log10(+/- 0) = -Infinity. |
| (...skipping 27 matching lines...) Expand all Loading... |
| 902 // fdlibm does not have an explicit log2 function, but fdlibm's pow | 902 // fdlibm does not have an explicit log2 function, but fdlibm's pow |
| 903 // function does implement an accurate log2 function as part of the | 903 // function does implement an accurate log2 function as part of the |
| 904 // pow implementation. This extracts the core parts of that as a | 904 // pow implementation. This extracts the core parts of that as a |
| 905 // separate log2 function. | 905 // separate log2 function. |
| 906 | 906 |
| 907 // Method: | 907 // Method: |
| 908 // Compute log2(x) in two pieces: | 908 // Compute log2(x) in two pieces: |
| 909 // log2(x) = w1 + w2 | 909 // log2(x) = w1 + w2 |
| 910 // where w1 has 53-24 = 29 bits of trailing zeroes. | 910 // where w1 has 53-24 = 29 bits of trailing zeroes. |
| 911 | 911 |
| 912 const DP_H = kMath[64]; | 912 define DP_H = kMath[64]; |
| 913 const DP_L = kMath[65]; | 913 define DP_L = kMath[65]; |
| 914 | 914 |
| 915 // Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3) | 915 // Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3) |
| 916 macro KLOG2(x) | 916 macro KLOG2(x) |
| 917 (kMath[55+x]) | 917 (kMath[55+x]) |
| 918 endmacro | 918 endmacro |
| 919 | 919 |
| 920 // cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy. | 920 // cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy. |
| 921 const CP = kMath[61]; | 921 define CP = kMath[61]; |
| 922 const CP_H = kMath[62]; | 922 define CP_H = kMath[62]; |
| 923 const CP_L = kMath[63]; | 923 define CP_L = kMath[63]; |
| 924 // 2^53 | 924 // 2^53 |
| 925 const TWO53 = 9007199254740992; | 925 define TWO53 = 9007199254740992; |
| 926 | 926 |
| 927 function MathLog2(x) { | 927 function MathLog2(x) { |
| 928 x = x * 1; // Convert to number. | 928 x = x * 1; // Convert to number. |
| 929 var ax = $abs(x); | 929 var ax = $abs(x); |
| 930 var hx = %_DoubleHi(x); | 930 var hx = %_DoubleHi(x); |
| 931 var lx = %_DoubleLo(x); | 931 var lx = %_DoubleLo(x); |
| 932 var ix = hx & 0x7fffffff; | 932 var ix = hx & 0x7fffffff; |
| 933 | 933 |
| 934 // Handle special cases. | 934 // Handle special cases. |
| 935 // log2(+/- 0) = -Infinity | 935 // log2(+/- 0) = -Infinity |
| (...skipping 85 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 1021 "log10", MathLog10, | 1021 "log10", MathLog10, |
| 1022 "log2", MathLog2, | 1022 "log2", MathLog2, |
| 1023 "log1p", MathLog1p, | 1023 "log1p", MathLog1p, |
| 1024 "expm1", MathExpm1 | 1024 "expm1", MathExpm1 |
| 1025 ]); | 1025 ]); |
| 1026 | 1026 |
| 1027 %SetInlineBuiltinFlag(MathSin); | 1027 %SetInlineBuiltinFlag(MathSin); |
| 1028 %SetInlineBuiltinFlag(MathCos); | 1028 %SetInlineBuiltinFlag(MathCos); |
| 1029 | 1029 |
| 1030 })(); | 1030 })(); |
| OLD | NEW |