From bb683c9291592144f1a51745341f1294d14cd553 Mon Sep 17 00:00:00 2001 From: sapics Date: Mon, 20 Jun 2016 09:52:09 +0900 Subject: [PATCH 1/5] Remove same value solution in solveCubic --- src/util/Numerical.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index b81718d8..42457699 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -356,8 +356,8 @@ var Numerical = new function() { } // The cubic has been deflated to a quadratic. var count = Numerical.solveQuadratic(a, b1, c2, roots, min, max); - if (isFinite(x) && count >= 0 - && (count === 0 || x !== roots[count - 1]) + if (isFinite(x) && (count === 0 + || count > 0 && x !== roots[0] && x !== roots[1]) && (min == null || x > min - EPSILON && x < max + EPSILON)) roots[count++] = min == null ? x : clamp(x, min, max); return count; From a9cc938967cff33ac3b04e84175cf326869c5a68 Mon Sep 17 00:00:00 2001 From: sapics Date: Mon, 20 Jun 2016 13:53:39 +0900 Subject: [PATCH 2/5] Fix normalization in solveQuadratic --- src/util/Numerical.js | 1 + 1 file changed, 1 insertion(+) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index 42457699..4aa5f006 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -238,6 +238,7 @@ var Numerical = new function() { abs(Math.floor(Math.log(gmC) * Math.LOG10E))); a *= mult; b *= mult; + B *= mult; c *= mult; // Recalculate the discriminant D = b * b - a * c; From 78f65c9fabf38f896f024b200dfdd2dc41c99d08 Mon Sep 17 00:00:00 2001 From: sapics Date: Mon, 20 Jun 2016 15:41:36 +0900 Subject: [PATCH 3/5] Improve solveQuadratic and solveCubic by hkrish c-code --- src/util/Numerical.js | 74 ++++++++++++++++++++++++++++++++----------- 1 file changed, 55 insertions(+), 19 deletions(-) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index 4aa5f006..50aec173 100644 --- a/src/util/Numerical.js +++ b/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); - if (gmC < 1e-8) { - // We multiply with a factor to normalize the coefficients. - // The factor is just the nearest exponent of 10, big enough - // to raise all the coefficients to nearly [-1, +1] range. - var mult = gmC === 0 ? 0 : pow(10, - abs(Math.floor(Math.log(gmC) * Math.LOG10E))); - a *= mult; - b *= mult; - B *= mult; - c *= mult; - // Recalculate the discriminant - D = b * b - a * c; + var sc = (abs(a) + abs(b) + abs(c)) || MACHINE_EPSILON; + sc = pow(2, -Math.floor(Math.log(sc) / Math.log(2) + 0.5)); + a *= sc; + b *= sc; + c *= sc; + // Recalculate the discriminant + D = b * b - a * c; + E = b * b + a * c; + B = - 2.0 * b; + if (pi * abs(D) < E) { + D = higherPrecisionDiscriminant(a, b, 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 */ solveCubic: function(a, b, c, d, roots, min, max) { - var count = 0, - x, b1, c2; + var count = 0, 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 { From 645e2c2af3d0cc4c446c5537aebd9b8031482973 Mon Sep 17 00:00:00 2001 From: sapics Date: Mon, 20 Jun 2016 17:27:08 +0900 Subject: [PATCH 4/5] Revert EPSILON error in solveCubic --- src/util/Numerical.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index 50aec173..bb3c4df0 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -341,12 +341,12 @@ var Numerical = new function() { } // If a or d is zero, we only need to solve a quadratic, so we set // the coefficients appropriately. - if (a === 0) { + if (abs(a) < EPSILON) { a = b; b1 = c; c2 = d; x = Infinity; - } else if (d === 0) { + } else if (abs(d) < EPSILON) { b1 = b; c2 = c; x = 0; From 4fd120fab8e664a8d6748d10cb895d4a39a968da Mon Sep 17 00:00:00 2001 From: sapics Date: Tue, 21 Jun 2016 08:47:42 +0900 Subject: [PATCH 5/5] Minor optimization in Numerical.js --- src/util/Numerical.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index bb3c4df0..9c339b89 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -255,7 +255,7 @@ var Numerical = new function() { if (D !== 0 && abs(D) < MACHINE_EPSILON) { // If the geometric mean of the coefficients is small enough var sc = (abs(a) + abs(b) + abs(c)) || MACHINE_EPSILON; - sc = pow(2, -Math.floor(Math.log(sc) / Math.log(2) + 0.5)); + sc = pow(2, -Math.floor(Math.log(sc) * Math.LOG2E + 0.5)); a *= sc; b *= sc; c *= sc; @@ -333,7 +333,7 @@ var Numerical = new function() { // 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))); + var p = pow(2, -Math.floor(Math.log(s) * Math.LOG2E)); a *= p; b *= p; c *= p;