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