| 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 | 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 | 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. | 21 // from decimal to binary accurately enough to produce the intended values. |
| 22 // kMath is initialized to a Float64Array during genesis and not writable. | 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 |
| 24 // initialized to a two-element Float64Array during genesis. |
| 23 | 25 |
| 24 "use strict"; | 26 "use strict"; |
| 25 | 27 |
| 26 var kMath; | 28 var kMath; |
| 29 var rempio2result; |
| 27 | 30 |
| 28 const INVPIO2 = kMath[0]; | 31 const INVPIO2 = kMath[0]; |
| 29 const PIO2_1 = kMath[1]; | 32 const PIO2_1 = kMath[1]; |
| 30 const PIO2_1T = kMath[2]; | 33 const PIO2_1T = kMath[2]; |
| 31 const PIO2_2 = kMath[3]; | 34 const PIO2_2 = kMath[3]; |
| 32 const PIO2_2T = kMath[4]; | 35 const PIO2_2T = kMath[4]; |
| 33 const PIO2_3 = kMath[5]; | 36 const PIO2_3 = kMath[5]; |
| 34 const PIO2_3T = kMath[6]; | 37 const PIO2_3T = kMath[6]; |
| 35 const PIO4 = kMath[32]; | 38 const PIO4 = kMath[32]; |
| 36 const PIO4LO = kMath[33]; | 39 const PIO4LO = kMath[33]; |
| 37 | 40 |
| 38 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For | 41 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For |
| 39 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 | 42 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 |
| 40 // to more than double precision. | 43 // to more than double precision. |
| 44 |
| 41 macro REMPIO2(X) | 45 macro REMPIO2(X) |
| 42 var n, y0, y1; | 46 var n, y0, y1; |
| 43 var hx = %_DoubleHi(X); | 47 var hx = %_DoubleHi(X); |
| 44 var ix = hx & 0x7fffffff; | 48 var ix = hx & 0x7fffffff; |
| 45 | 49 |
| 46 if (ix < 0x4002d97c) { | 50 if (ix < 0x4002d97c) { |
| 47 // |X| ~< 3*pi/4, special case with n = +/- 1 | 51 // |X| ~< 3*pi/4, special case with n = +/- 1 |
| 48 if (hx > 0) { | 52 if (hx > 0) { |
| 49 var z = X - PIO2_1; | 53 var z = X - PIO2_1; |
| 50 if (ix != 0x3ff921fb) { | 54 if (ix != 0x3ff921fb) { |
| (...skipping 47 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 98 } | 102 } |
| 99 } | 103 } |
| 100 y1 = (r - y0) - w; | 104 y1 = (r - y0) - w; |
| 101 if (hx < 0) { | 105 if (hx < 0) { |
| 102 n = -n; | 106 n = -n; |
| 103 y0 = -y0; | 107 y0 = -y0; |
| 104 y1 = -y1; | 108 y1 = -y1; |
| 105 } | 109 } |
| 106 } else { | 110 } else { |
| 107 // Need to do full Payne-Hanek reduction here. | 111 // Need to do full Payne-Hanek reduction here. |
| 108 var r = %RemPiO2(X); | 112 n = %RemPiO2(X, rempio2result); |
| 109 n = r[0]; | 113 y0 = rempio2result[0]; |
| 110 y0 = r[1]; | 114 y1 = rempio2result[1]; |
| 111 y1 = r[2]; | |
| 112 } | 115 } |
| 113 endmacro | 116 endmacro |
| 114 | 117 |
| 115 | 118 |
| 116 // __kernel_sin(X, Y, IY) | 119 // __kernel_sin(X, Y, IY) |
| 117 // kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 | 120 // kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 |
| 118 // Input X is assumed to be bounded by ~pi/4 in magnitude. | 121 // Input X is assumed to be bounded by ~pi/4 in magnitude. |
| 119 // Input Y is the tail of X so that x = X + Y. | 122 // Input Y is the tail of X so that x = X + Y. |
| 120 // | 123 // |
| 121 // Algorithm | 124 // Algorithm |
| (...skipping 865 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 987 z_l = CP_L * p_h + p_l * CP + dp_l; | 990 z_l = CP_L * p_h + p_l * CP + dp_l; |
| 988 | 991 |
| 989 // log2(ax) = (ss + ...) * 2 / (3 * log(2)) = n + dp_h + z_h + z_l | 992 // log2(ax) = (ss + ...) * 2 / (3 * log(2)) = n + dp_h + z_h + z_l |
| 990 var t = n; | 993 var t = n; |
| 991 var t1 = %_ConstructDouble(%_DoubleHi(((z_h + z_l) + dp_h) + t), 0); | 994 var t1 = %_ConstructDouble(%_DoubleHi(((z_h + z_l) + dp_h) + t), 0); |
| 992 var t2 = z_l - (((t1 - t) - dp_h) - z_h); | 995 var t2 = z_l - (((t1 - t) - dp_h) - z_h); |
| 993 | 996 |
| 994 // t1 + t2 = log2(ax), sum up because we do not care about extra precision. | 997 // t1 + t2 = log2(ax), sum up because we do not care about extra precision. |
| 995 return t1 + t2; | 998 return t1 + t2; |
| 996 } | 999 } |
| OLD | NEW |