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 |