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 // rempio2result is used as a container for return values of %RemPiO2. It is | |
20 // initialized to a two-element Float64Array during genesis. | |
21 | |
22 (function(global, utils) { | 19 (function(global, utils) { |
23 | 20 |
24 "use strict"; | 21 "use strict"; |
25 | 22 |
26 %CheckIsBootstrapping(); | 23 %CheckIsBootstrapping(); |
27 | 24 |
28 // ------------------------------------------------------------------- | 25 // ------------------------------------------------------------------- |
29 // Imports | 26 // Imports |
30 | 27 |
31 var GlobalFloat64Array = global.Float64Array; | |
32 var GlobalMath = global.Math; | 28 var GlobalMath = global.Math; |
33 var MathAbs; | 29 var MathAbs; |
34 var MathExpm1; | 30 var MathExpm1; |
35 var NaN = %GetRootNaN(); | |
36 var rempio2result; | |
37 | 31 |
38 utils.Import(function(from) { | 32 utils.Import(function(from) { |
39 MathAbs = from.MathAbs; | 33 MathAbs = from.MathAbs; |
40 MathExpm1 = from.MathExpm1; | 34 MathExpm1 = from.MathExpm1; |
41 }); | 35 }); |
42 | 36 |
43 utils.CreateDoubleResultArray = function(global) { | |
44 rempio2result = new GlobalFloat64Array(2); | |
45 }; | |
46 | |
47 // ------------------------------------------------------------------- | |
48 | |
49 define INVPIO2 = 6.36619772367581382433e-01; | |
50 define PIO2_1 = 1.57079632673412561417; | |
51 define PIO2_1T = 6.07710050650619224932e-11; | |
52 define PIO2_2 = 6.07710050630396597660e-11; | |
53 define PIO2_2T = 2.02226624879595063154e-21; | |
54 define PIO2_3 = 2.02226624871116645580e-21; | |
55 define PIO2_3T = 8.47842766036889956997e-32; | |
56 define PIO4 = 7.85398163397448278999e-01; | |
57 define PIO4LO = 3.06161699786838301793e-17; | |
58 | |
59 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For | |
60 // precision, r is returned as two values y0 and y1 such that r = y0 + y1 | |
61 // to more than double precision. | |
62 | |
63 macro REMPIO2(X) | |
64 var n, y0, y1; | |
65 var hx = %_DoubleHi(X); | |
66 var ix = hx & 0x7fffffff; | |
67 | |
68 if (ix < 0x4002d97c) { | |
69 // |X| ~< 3*pi/4, special case with n = +/- 1 | |
70 if (hx > 0) { | |
71 var z = X - PIO2_1; | |
72 if (ix != 0x3ff921fb) { | |
73 // 33+53 bit pi is good enough | |
74 y0 = z - PIO2_1T; | |
75 y1 = (z - y0) - PIO2_1T; | |
76 } else { | |
77 // near pi/2, use 33+33+53 bit pi | |
78 z -= PIO2_2; | |
79 y0 = z - PIO2_2T; | |
80 y1 = (z - y0) - PIO2_2T; | |
81 } | |
82 n = 1; | |
83 } else { | |
84 // Negative X | |
85 var z = X + PIO2_1; | |
86 if (ix != 0x3ff921fb) { | |
87 // 33+53 bit pi is good enough | |
88 y0 = z + PIO2_1T; | |
89 y1 = (z - y0) + PIO2_1T; | |
90 } else { | |
91 // near pi/2, use 33+33+53 bit pi | |
92 z += PIO2_2; | |
93 y0 = z + PIO2_2T; | |
94 y1 = (z - y0) + PIO2_2T; | |
95 } | |
96 n = -1; | |
97 } | |
98 } else if (ix <= 0x413921fb) { | |
99 // |X| ~<= 2^19*(pi/2), medium size | |
100 var t = MathAbs(X); | |
101 n = (t * INVPIO2 + 0.5) | 0; | |
102 var r = t - n * PIO2_1; | |
103 var w = n * PIO2_1T; | |
104 // First round good to 85 bit | |
105 y0 = r - w; | |
106 if (ix - (%_DoubleHi(y0) & 0x7ff00000) > 0x1000000) { | |
107 // 2nd iteration needed, good to 118 | |
108 t = r; | |
109 w = n * PIO2_2; | |
110 r = t - w; | |
111 w = n * PIO2_2T - ((t - r) - w); | |
112 y0 = r - w; | |
113 if (ix - (%_DoubleHi(y0) & 0x7ff00000) > 0x3100000) { | |
114 // 3rd iteration needed. 151 bits accuracy | |
115 t = r; | |
116 w = n * PIO2_3; | |
117 r = t - w; | |
118 w = n * PIO2_3T - ((t - r) - w); | |
119 y0 = r - w; | |
120 } | |
121 } | |
122 y1 = (r - y0) - w; | |
123 if (hx < 0) { | |
124 n = -n; | |
125 y0 = -y0; | |
126 y1 = -y1; | |
127 } | |
128 } else { | |
129 // Need to do full Payne-Hanek reduction here. | |
130 n = %RemPiO2(X, rempio2result); | |
131 y0 = rempio2result[0]; | |
132 y1 = rempio2result[1]; | |
133 } | |
134 endmacro | |
135 | |
136 | |
137 // kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 | |
138 // Input x is assumed to be bounded by ~pi/4 in magnitude. | |
139 // Input y is the tail of x. | |
140 // Input k indicates whether ieee_tan (if k = 1) or -1/tan (if k = -1) | |
141 // is returned. | |
142 // | |
143 // Algorithm | |
144 // 1. Since ieee_tan(-x) = -ieee_tan(x), we need only to consider positive x. | |
145 // 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0. | |
146 // 3. ieee_tan(x) is approximated by a odd polynomial of degree 27 on | |
147 // [0,0.67434] | |
148 // 3 27 | |
149 // tan(x) ~ x + T1*x + ... + T13*x | |
150 // where | |
151 // | |
152 // |ieee_tan(x) 2 4 26 | -59.2 | |
153 // |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 | |
154 // | x | | |
155 // | |
156 // Note: ieee_tan(x+y) = ieee_tan(x) + tan'(x)*y | |
157 // ~ ieee_tan(x) + (1+x*x)*y | |
158 // Therefore, for better accuracy in computing ieee_tan(x+y), let | |
159 // 3 2 2 2 2 | |
160 // r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) | |
161 // then | |
162 // 3 2 | |
163 // tan(x+y) = x + (T1*x + (x *(r+y)+y)) | |
164 // | |
165 // 4. For x in [0.67434,pi/4], let y = pi/4 - x, then | |
166 // tan(x) = ieee_tan(pi/4-y) = (1-ieee_tan(y))/(1+ieee_tan(y)) | |
167 // = 1 - 2*(ieee_tan(y) - (ieee_tan(y)^2)/(1+ieee_tan(y))) | |
168 // | |
169 // Set returnTan to 1 for tan; -1 for cot. Anything else is illegal | |
170 // and will cause incorrect results. | |
171 // | |
172 define T00 = 3.33333333333334091986e-01; | |
173 define T01 = 1.33333333333201242699e-01; | |
174 define T02 = 5.39682539762260521377e-02; | |
175 define T03 = 2.18694882948595424599e-02; | |
176 define T04 = 8.86323982359930005737e-03; | |
177 define T05 = 3.59207910759131235356e-03; | |
178 define T06 = 1.45620945432529025516e-03; | |
179 define T07 = 5.88041240820264096874e-04; | |
180 define T08 = 2.46463134818469906812e-04; | |
181 define T09 = 7.81794442939557092300e-05; | |
182 define T10 = 7.14072491382608190305e-05; | |
183 define T11 = -1.85586374855275456654e-05; | |
184 define T12 = 2.59073051863633712884e-05; | |
185 | |
186 function KernelTan(x, y, returnTan) { | |
187 var z; | |
188 var w; | |
189 var hx = %_DoubleHi(x); | |
190 var ix = hx & 0x7fffffff; | |
191 | |
192 if (ix < 0x3e300000) { // |x| < 2^-28 | |
193 if (((ix | %_DoubleLo(x)) | (returnTan + 1)) == 0) { | |
194 // x == 0 && returnTan = -1 | |
195 return 1 / MathAbs(x); | |
196 } else { | |
197 if (returnTan == 1) { | |
198 return x; | |
199 } else { | |
200 // Compute -1/(x + y) carefully | |
201 var w = x + y; | |
202 var z = %_ConstructDouble(%_DoubleHi(w), 0); | |
203 var v = y - (z - x); | |
204 var a = -1 / w; | |
205 var t = %_ConstructDouble(%_DoubleHi(a), 0); | |
206 var s = 1 + t * z; | |
207 return t + a * (s + t * v); | |
208 } | |
209 } | |
210 } | |
211 if (ix >= 0x3fe59428) { // |x| > .6744 | |
212 if (x < 0) { | |
213 x = -x; | |
214 y = -y; | |
215 } | |
216 z = PIO4 - x; | |
217 w = PIO4LO - y; | |
218 x = z + w; | |
219 y = 0; | |
220 } | |
221 z = x * x; | |
222 w = z * z; | |
223 | |
224 // Break x^5 * (T1 + x^2*T2 + ...) into | |
225 // x^5 * (T1 + x^4*T3 + ... + x^20*T11) + | |
226 // x^5 * (x^2 * (T2 + x^4*T4 + ... + x^22*T12)) | |
227 var r = T01 + w * (T03 + w * (T05 + | |
228 w * (T07 + w * (T09 + w * T11)))); | |
229 var v = z * (T02 + w * (T04 + w * (T06 + | |
230 w * (T08 + w * (T10 + w * T12))))); | |
231 var s = z * x; | |
232 r = y + z * (s * (r + v) + y); | |
233 r = r + T00 * s; | |
234 w = x + r; | |
235 if (ix >= 0x3fe59428) { | |
236 return (1 - ((hx >> 30) & 2)) * | |
237 (returnTan - 2.0 * (x - (w * w / (w + returnTan) - r))); | |
238 } | |
239 if (returnTan == 1) { | |
240 return w; | |
241 } else { | |
242 z = %_ConstructDouble(%_DoubleHi(w), 0); | |
243 v = r - (z - x); | |
244 var a = -1 / w; | |
245 var t = %_ConstructDouble(%_DoubleHi(a), 0); | |
246 s = 1 + t * z; | |
247 return t + a * (s + t * v); | |
248 } | |
249 } | |
250 | |
251 // ECMA 262 - 15.8.2.18 | |
252 function MathTan(x) { | |
253 x = x * 1; // Convert to number. | |
254 if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) { | |
255 // |x| < pi/4, approximately. No reduction needed. | |
256 return KernelTan(x, 0, 1); | |
257 } | |
258 REMPIO2(x); | |
259 return KernelTan(y0, y1, (n & 1) ? -1 : 1); | |
260 } | |
261 | |
262 define LN2_HI = 6.93147180369123816490e-01; | |
263 define LN2_LO = 1.90821492927058770002e-10; | |
264 | |
265 // 2^54 | |
266 define TWO54 = 18014398509481984; | |
267 | |
268 // ES6 draft 09-27-13, section 20.2.2.30. | 37 // ES6 draft 09-27-13, section 20.2.2.30. |
269 // Math.sinh | 38 // Math.sinh |
270 // Method : | 39 // Method : |
271 // mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 | 40 // mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 |
272 // 1. Replace x by |x| (sinh(-x) = -sinh(x)). | 41 // 1. Replace x by |x| (sinh(-x) = -sinh(x)). |
273 // 2. | 42 // 2. |
274 // E + E/(E+1) | 43 // E + E/(E+1) |
275 // 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) | 44 // 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) |
276 // 2 | 45 // 2 |
277 // | 46 // |
(...skipping 141 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
419 } else { | 188 } else { |
420 // |x| > 22, return +/- 1 | 189 // |x| > 22, return +/- 1 |
421 z = 1; | 190 z = 1; |
422 } | 191 } |
423 return (x >= 0) ? z : -z; | 192 return (x >= 0) ? z : -z; |
424 } | 193 } |
425 | 194 |
426 //------------------------------------------------------------------- | 195 //------------------------------------------------------------------- |
427 | 196 |
428 utils.InstallFunctions(GlobalMath, DONT_ENUM, [ | 197 utils.InstallFunctions(GlobalMath, DONT_ENUM, [ |
429 "tan", MathTan, | |
430 "sinh", MathSinh, | 198 "sinh", MathSinh, |
431 "cosh", MathCosh, | 199 "cosh", MathCosh, |
432 "tanh", MathTanh | 200 "tanh", MathTanh |
433 ]); | 201 ]); |
434 | 202 |
435 }) | 203 }) |
OLD | NEW |