Line data Source code
1 : /**
2 : * Safe, platform consistent implementations of some Math functions
3 : *
4 : * These functions are implemented in JS to avoid observed differences
5 : * between results of different floating point libraries, see
6 : * https://bugzilla.mozilla.org/show_bug.cgi?id=531915
7 : *
8 : * They mostly meet the ECMAScript Edition 5 spec, see
9 : * https://www.ecma-international.org/publications/files/ECMA-ST/Ecma-262.pdf
10 : *
11 : * See simulation/components/tests/test_Math.js for tests.
12 : */
13 :
14 : /**
15 : * Approximation of cosine of a (radians)
16 : */
17 72 : Math.cos = function(a)
18 : {
19 : // Bring a into the 0 to +pi range without expensive branching.
20 : // Uses the symmetry that cos is even.
21 1158 : a = (a + Math.PI) % (2*Math.PI);
22 1158 : a = Math.abs((2*Math.PI + a) % (2*Math.PI) - Math.PI);
23 :
24 : // make b = 0 if a < pi/2 and b=1 if a > pi/2
25 1158 : var b = (a-Math.PI/2) + Math.abs(a-Math.PI/2);
26 1158 : b = b/(b+1e-30); // normalize b to one while avoiding divide by zero errors.
27 :
28 : // if a > pi/2 send a to pi-a, otherwise just send a to -a which has no effect
29 : // Using the symmetry cos(x) = -cos(pi-x) to bring a to the 0 to pi/2 range.
30 1158 : a = b*Math.PI - a;
31 :
32 1158 : var c = 1 - 2*b; // sign of the output
33 :
34 : // Taylor expansion about 0 with a correction term in the quadratic to make cos(pi/2)=0
35 1158 : return c * (1 - a*a*(0.5000000025619951 - a*a*(1/24 - a*a*(1/720 - a*a*(1/40320 - a*a*(1/3628800 - a*a/479001600))))));
36 : };
37 :
38 : /**
39 : * Approximation of sine of a (radians)
40 : */
41 72 : Math.sin = function(a)
42 : {
43 579 : return Math.cos(a - Math.PI/2);
44 : };
45 :
46 : /**
47 : * Approximation of arctangent of a, returns angle from -pi/2 to pi/2
48 : */
49 72 : Math.atan = function(a)
50 : {
51 24 : var tanPiBy6 = 0.5773502691896257;
52 24 : var tanPiBy12 = 0.2679491924311227;
53 24 : var sign = 1;
54 24 : var inverted = false;
55 24 : var tanPiBy6Shift = 0;
56 :
57 24 : if (a < 0 || 1/a === -Infinity)
58 : {
59 : // tan(x) = -tan(-x) so remove sign now and put it back at the end
60 2 : sign = -1;
61 2 : a *= -1;
62 : }
63 :
64 24 : if (a > 1)
65 : {
66 : // tan(pi/2 - x) = 1/tan(x)
67 9 : inverted = true;
68 9 : a = 1/a;
69 : }
70 :
71 24 : if (a > tanPiBy12)
72 : {
73 : // tan(x-pi/6) = (tan(x) - tan(pi/6)) / (1 + tan(pi/6)tan(x))
74 8 : tanPiBy6Shift = Math.PI/6;
75 8 : a = (a - tanPiBy6) / (1 + tanPiBy6*a);
76 : }
77 : // Now a will be in the range [-tan(pi/12), tan(pi/12)]
78 :
79 : // Use the taylor expansion around 0 with a correction to the linear term to match the pi/12 boundary
80 : // atan(x) = x - x^3/3 + x^5/5 - ...
81 24 : var r = a*(1.0000000000390272 - a*a*(1/3 - a*a*(1/5 - a*a*(1/7 - a*a*(1/9 - a*a*(1/11 - a*a*(1/13 - a*a/15)))))));
82 :
83 : // shift the result back where necessary
84 24 : r += tanPiBy6Shift;
85 24 : if (inverted)
86 9 : r = Math.PI/2 - r;
87 24 : return sign * r;
88 : };
89 :
90 : /**
91 : * Approximation of arctangent of y/x, returns angle from -pi to pi
92 : */
93 72 : Math.atan2 = function(y,x)
94 : {
95 : // get unsigned x,y for ease of calculation, this means all angles are in the range [0, pi/2]
96 41 : var ux = Math.abs(x);
97 41 : var uy = Math.abs(y);
98 : // holds the result in the upper right quadrant
99 : var r;
100 :
101 : // Handle all edges cases to match the spec
102 41 : if (uy === 0)
103 14 : r = 0;
104 : else
105 : {
106 27 : if (ux === 0)
107 6 : r = Math.PI / 2;
108 :
109 27 : if (uy === Infinity)
110 : {
111 6 : if (ux === Infinity)
112 4 : r = Math.PI / 4;
113 : else
114 2 : r = Math.PI / 2;
115 : }
116 : else
117 : {
118 21 : if (ux === Infinity)
119 4 : r = 0;
120 : else
121 17 : r = Math.atan(uy/ux);
122 : }
123 : }
124 :
125 : // puts the result into the correct quadrant
126 : // 1/(-0) is the only way to determine the sign for a 0 value
127 41 : if (x < 0 || 1/x === -Infinity)
128 : {
129 11 : if (y < 0 || 1/y === -Infinity)
130 5 : return -Math.PI + r;
131 : else
132 6 : return Math.PI - r;
133 : }
134 : else
135 : {
136 30 : if (y < 0 || 1/y === -Infinity)
137 8 : return -r;
138 : else
139 22 : return r;
140 : }
141 : };
142 :
143 72 : Math.acos = function()
144 : {
145 0 : error("Math.acos() does not yet have a synchronization safe implementation");
146 : };
147 :
148 72 : Math.asin = function()
149 : {
150 0 : error("Math.asin() does not yet have a synchronization safe implementation");
151 : };
152 :
153 72 : Math.tan = function()
154 : {
155 0 : error("Math.tan() does not yet have a synchronization safe implementation");
156 : };
157 :
158 : /**
159 : * Approximation of raising x to the power y
160 : */
161 72 : Math.pow = function(x, y)
162 : {
163 185 : if (Math.round(y) === y)
164 : {
165 150 : if (y >= 0)
166 140 : return Math.intPow(x, y);
167 :
168 10 : return 1 / Math.intPow(x, -y);
169 : }
170 : // log has the biggest error when x ~=~ 1
171 : // exp has the biggest error when y*log(x) ~<~ 0
172 : // so the biggest error happens around numbers like pow(0.9999,0.0001),
173 : // that has an error of 10^-17. So I think we're safe
174 35 : return Math.exp(y*Math.log(x));
175 : };
176 :
177 : /**
178 : * Get the square of a number without repeating the value and without calling the slower Math.pow.
179 : */
180 72 : Math.square = function(x)
181 : {
182 83897 : return x * x;
183 : };
184 :
185 : /**
186 : * Approximation of the exponential function, e raised to the power x
187 : */
188 72 : Math.exp = function(x)
189 : {
190 41 : if (x < 0)
191 4 : var iPart = 1/Math.intPow(Math.E, -Math.floor(x));
192 : else
193 37 : var iPart = Math.intPow(Math.E, Math.floor(x));
194 :
195 41 : if (x === Math.floor(x))
196 : // no need to loop if we know the answer
197 11 : return iPart;
198 :
199 : // the integer part is known, work further with the decimal part of x
200 30 : x = x - Math.floor(x); // x \in [0,1)
201 :
202 : // taylor series around 0
203 : // max error ~=~ 10^(-16)
204 30 : var dPart = 1;
205 :
206 30 : for (var i = 22; i > 0; i--)
207 660 : dPart = 1+x*dPart/i;
208 :
209 : // total precision ~=~ 17 decimal digits
210 30 : return iPart*dPart;
211 : };
212 :
213 : /**
214 : * Approximation of the natural logarithm of x
215 : *
216 : * For values very close to 1, the error of 10^-16 could become bigger than the actual value
217 : * But this also happens with the native log function
218 : */
219 72 : Math.log = function(x)
220 : {
221 44 : if (!(x >= 0))
222 3 : return NaN;
223 :
224 41 : if (x === 0)
225 6 : return -Infinity;
226 :
227 35 : if (x === Infinity)
228 3 : return x;
229 :
230 : // start with calculating the binary logarithm
231 : // based on https://en.wikipedia.org/wiki/Binary_logarithm#Real_number
232 :
233 : // calculate to 50 fractional bits -> error ~=~ 10^-16
234 32 : var precisionBits = 50;
235 :
236 : // calculate integer log, rounded down
237 : // when implemented in C, just count the number of bits before the fraction
238 : // without leading zeros. This may be negative.
239 32 : var log = 0;
240 32 : if (x >= 1)
241 : {
242 31 : for (var i = 1; i <= x; i *= 2)
243 64 : log++;
244 :
245 31 : log--;
246 31 : i /= 2;
247 : }
248 : else
249 : {
250 1 : for (var i = 1; i > x; i /= 2)
251 1 : log--;
252 : }
253 : // now lb(x) = log + lb(y) with y = x/i. So y \in [1,2)
254 :
255 32 : var y = x/i;
256 :
257 : // if we're done, or there's a minimal rounding error and we should be done
258 : // convert to natural logarithm
259 32 : if (y <= 1)
260 17 : return log / Math.LOG2E;
261 :
262 15 : var m = 0;
263 15 : var add = 1;
264 15 : while (true)
265 : {
266 1118 : while (m <= precisionBits && y < 4)
267 : {
268 765 : m++;
269 765 : y *= y;
270 765 : add /= 2;
271 : }
272 1118 : if (m > precisionBits)
273 15 : break;
274 1103 : log += add;
275 1103 : y /= 2;
276 : }
277 :
278 : // convert binary logarithm to natural logarithm;
279 15 : return log / Math.LOG2E;
280 :
281 : };
282 :
283 : /**
284 : * Calculate the power for positive integer exponents
285 : */
286 72 : Math.intPow = function(x, y)
287 : {
288 191 : if (Math.abs(y) === Infinity)
289 : {
290 20 : if (Math.abs(x) === 1)
291 4 : return NaN;
292 :
293 16 : if (Math.abs(x) < 1 && y > 0 || Math.abs(x) > 1 && y < 0)
294 4 : return 0;
295 :
296 12 : return Infinity;
297 : }
298 171 : var powers = [x];
299 171 : var binary = [1];
300 171 : var i = 0;
301 171 : for (var e = 2; e <= y; e *= 2)
302 : {
303 : // calculate x^i, using x^(i/2)
304 2093 : powers.push(powers[i]*powers[i]);
305 2093 : binary.push(e);
306 2093 : i++;
307 : }
308 :
309 171 : var result = 1;
310 :
311 171 : var i = binary.length;
312 :
313 171 : while (y > 0)
314 : {
315 159 : if (binary[--i] <= y)
316 : {
317 94 : result *= powers[i];
318 94 : y -= binary[i];
319 : }
320 : }
321 : // error margin = 0 (default JS error)
322 171 : return result;
323 :
324 : };
325 :
326 72 : Math.euclidDistance2DSquared = function(x1, y1, x2, y2)
327 : {
328 41929 : return Math.square(x2 - x1) + Math.square(y2 - y1);
329 : };
330 :
331 : /**
332 : * Can be faster than Math.hypot.
333 : */
334 72 : Math.euclidDistance2D = function(x1, y1, x2, y2)
335 : {
336 9 : return Math.sqrt(Math.euclidDistance2DSquared(x1, y1, x2, y2));
337 : };
338 :
339 72 : Math.euclidDistance3DSquared = function(x1, y1, z1, x2, y2, z2)
340 : {
341 6 : return Math.square(x2 - x1) + Math.square(y2 - y1) + Math.square(z2 - z1);
342 : };
343 :
344 72 : Math.euclidDistance3D = function(x1, y1, z1, x2, y2, z2)
345 : {
346 6 : return Math.sqrt(Math.euclidDistance3DSquared(x1, y1, z1, x2, y2, z2));
347 : };
|
