LCOV - code coverage report
 Current view: top level - globalscripts - Math.js (source / functions) Hit Total Coverage Test: lcov.info Lines: 133 136 97.8 % Date: 2023-04-02 12:52:40 Functions: 13 16 81.2 %
  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 : }; 

 Generated by: LCOV version 1.14