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 |