Improve solveQuadratic and solveCubic by hkrish c-code

src/util/Numerical.js

@@ -67,6 +67,25 @@ var Numerical = new function() {
         return value < min ? min : value > max ? max : value;
     }
 
+    function splitDouble(X) {
+        var bigX = X * 134217729, // X*(2^27 + 1)
+            Y = X - bigX,
+            Xh = Y + bigX; // Don't optimize Y away!
+        return [Xh, X - Xh];
+    }
+
+    function higherPrecisionDiscriminant(a, b, c) {
+        var ad = splitDouble(a),
+            bd = splitDouble(b),
+            cd = splitDouble(c),
+            p = b * b,
+            dp = (bd[0] * bd[0] - p + 2 * bd[0] * bd[1]) + bd[1] * bd[1],
+            q = a * c,
+            dq = (ad[0] * cd[0] - q + ad[0] * cd[1] + ad[1] * cd[0])
+                    + ad[1] * cd[1];
+        return (p - q) + (dp - dq);
+    }
+
     return /** @lends Numerical */{
         TOLERANCE: 1e-6,
         /**
@@ -202,6 +221,8 @@ var Numerical = new function() {
          *  Kahan W. - "To Solve a Real Cubic Equation"
          *  http://www.cs.berkeley.edu/~wkahan/Math128/Cubic.pdf
          *  Blinn J. - "How to solve a Quadratic Equation"
+         *  Harikrishnan G.
+         *  https://gist.github.com/hkrish/9e0de1f121971ee0fbab281f5c986de9
          *
          * @param {Number} a the quadratic term
          * @param {Number} b the linear term
@@ -220,28 +241,30 @@ var Numerical = new function() {
                 eMax = max + EPSILON,
                 x1, x2 = Infinity,
                 B = b,
-                D;
+                D, E, pi = 3;
             // a, b, c are expected to be the coefficients of the equation:
             // Ax² - 2Bx + C == 0, so we take b = -B/2:
             b /= -2;
             D = b * b - a * c; // Discriminant
+            E = b * b + a * c;
+            if (pi * abs(D) < E) {
+                D = higherPrecisionDiscriminant(a, b, c);
+            }
             // If the discriminant is very small, we can try to pre-condition
             // the coefficients, so that we may get better accuracy
             if (D !== 0 && abs(D) < MACHINE_EPSILON) {
                 // If the geometric mean of the coefficients is small enough
-                var gmC = pow(abs(a * b * c), 1 / 3);
+                var sc = (abs(a) + abs(b) + abs(c)) || MACHINE_EPSILON;
-                if (gmC < 1e-8) {
+                sc = pow(2, -Math.floor(Math.log(sc) / Math.log(2) + 0.5));
-                    // We multiply with a factor to normalize the coefficients.
+                a *= sc;
-                    // The factor is just the nearest exponent of 10, big enough
+                b *= sc;
-                    // to raise all the coefficients to nearly [-1, +1] range.
+                c *= sc;
-                    var mult = gmC === 0 ? 0 : pow(10,
+                // Recalculate the discriminant
-                        abs(Math.floor(Math.log(gmC) * Math.LOG10E)));
+                D = b * b - a * c;
-                    a *= mult;
+                E = b * b + a * c;
-                    b *= mult;
+                B = - 2.0 * b;
-                    B *= mult;
+                if (pi * abs(D) < E) {
-                    c *= mult;
+                    D = higherPrecisionDiscriminant(a, b, c);
-                    // Recalculate the discriminant
-                    D = b * b - a * c;
                 }
             }
             if (abs(a) < EPSILON) {
@@ -284,6 +307,8 @@ var Numerical = new function() {
          * References:
          *  Kahan W. - "To Solve a Real Cubic Equation"
          *   http://www.cs.berkeley.edu/~wkahan/Math128/Cubic.pdf
+         *  Harikrishnan G.
+         *  https://gist.github.com/hkrish/9e0de1f121971ee0fbab281f5c986de9
          *
          * W. Kahan's paper contains inferences on accuracy of cubic
          * zero-finding methods. Also testing methods for robustness.
@@ -301,16 +326,27 @@ var Numerical = new function() {
          * @author Harikrishnan Gopalakrishnan <hari.exeption@gmail.com>
          */
         solveCubic: function(a, b, c, d, roots, min, max) {
-            var count = 0,
+            var count = 0, x, b1, c2,
-                x, b1, c2;
+                s = Math.max(abs(a), abs(b), abs(c), abs(d));
+            // Normalise coefficients a la Jenkins & Traub's RPOLY
+            if ((s < 1e-7 && s > 0) || s > 1e7) {
+                // Scale the coefficients by a multiple of the exponent of
+                // coefficients so that all the bits in the mantissa are
+                // preserved.
+                var p = pow(2, - Math.floor(Math.log(s) / Math.log(2)));
+                a *= p;
+                b *= p;
+                c *= p;
+                d *= p;
+            }
             // If a or d is zero, we only need to solve a quadratic, so we set
             // the coefficients appropriately.
-            if (abs(a) < EPSILON) {
+            if (a === 0) {
                 a = b;
                 b1 = c;
                 c2 = d;
                 x = Infinity;
-            } else if (abs(d) < EPSILON) {
+            } else if (d === 0) {
                 b1 = b;
                 c2 = c;
                 x = 0;
@@ -333,7 +369,7 @@ var Numerical = new function() {
                 s = t < 0 ? -1 : 1;
                 t = -qd / a;
                 // See Kahan's notes on why 1.324718*... works.
-                r = t > 0 ? 1.3247179572 * Math.max(r, sqrt(t)) : r;
+                r = t > 0 ? 1.324717957244746 * Math.max(r, sqrt(t)) : r;
                 x0 = x - s * r;
                 if (x0 !== x) {
                     do {