From f73717a7e774942675bda4308a4c535943eba0e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrg=20Lehni?= Date: Mon, 22 Oct 2012 18:21:33 -0400 Subject: [PATCH] Fix issues in Numerical.solveQuadratic(), solveCubic() and Path#contains(). Closes #71. --- src/path/Curve.js | 4 +- src/path/Path.js | 2 +- src/util/Numerical.js | 85 +++++++++++++++++++++++-------------------- 3 files changed, 48 insertions(+), 43 deletions(-) diff --git a/src/path/Curve.js b/src/path/Curve.js index 886390c5..dcc34934 100644 --- a/src/path/Curve.js +++ b/src/path/Curve.js @@ -319,7 +319,7 @@ var Curve = this.Curve = Base.extend(/** @lends Curve# */{ crossings = 0; for (var i = 0; i < num; i++) { var t = roots[i]; - if (t >= 0 && t < 1 && Curve.evaluate(vals, t, 0).x > point.x) { + if (t >= 0 && t <= 1 && Curve.evaluate(vals, t, 0).x > point.x) { // If we're close to 0 and are not changing y-direction from the // previous curve, do not count this root, as we're merely // touching a tip. Passing 1 for Curve.evaluate()'s type means @@ -665,7 +665,7 @@ var Curve = this.Curve = Base.extend(/** @lends Curve# */{ */ function toBezierForm(v, point) { var n = 3, // degree of B(t) - degree = 5, // degree of B(t) . P + degree = 5, // degree of B(t) . P c = [], d = [], cd = [], diff --git a/src/path/Path.js b/src/path/Path.js index ca4be3cc..696dcee1 100644 --- a/src/path/Path.js +++ b/src/path/Path.js @@ -558,7 +558,7 @@ var Path = this.Path = PathItem.extend(/** @lends Path# */{ */ removeSegments: function(from, to) { from = from || 0; - to = Base.pick(to, this._segments.length); + to = Base.pick(to, this._segments.length); var segments = this._segments, curves = this._curves, last = to >= segments.length, diff --git a/src/util/Numerical.js b/src/util/Numerical.js index 03af4fa5..bda03428 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -58,9 +58,15 @@ var Numerical = new function() { // Math short-cuts for often used methods and values var abs = Math.abs, sqrt = Math.sqrt, + pow = Math.pow, cos = Math.cos, PI = Math.PI; + // Define the missing Math.cbrt() + function cbrt(x) { + return x > 0 ? pow(x, 1 / 3) : x < 0 ? -pow(-x, 1 / 3) : 0; + } + return { TOLERANCE: 10e-6, // Precision when comparing against 0 @@ -121,33 +127,26 @@ var Numerical = new function() { * a*x^2 + b*x + c = 0 */ solveQuadratic: function(a, b, c, roots, tolerance) { - // After Numerical Recipes in C, 2nd edition, Press et al., - // 5.6, Quadratic and Cubic Equations + // Code ported over and adapted from Uintah library (MIT license). // If problem is actually linear, return 0 or 1 easy roots if (abs(a) < tolerance) { if (abs(b) >= tolerance) { roots[0] = -c / b; return 1; } - // If all the coefficients are 0, infinite values are - // possible! - if (abs(c) < tolerance) - return -1; // Infinite solutions - return 0; // 0 solutions + // If all the coefficients are 0, we have infinite solutions! + return abs(c) < tolerance ? -1 : 0; // Infinite or 0 solutions } var q = b * b - 4 * a * c; if (q < 0) return 0; // 0 solutions q = sqrt(q); - if (b < 0) - q = -q; - q = (b + q) * -0.5; + a *= 2; // Prepare division by (2 * a) var n = 0; - if (abs(q) >= tolerance) - roots[n++] = c / q; - if (abs(a) >= tolerance) - roots[n++] = q / a; - return n; // 0, 1 or 2 solutions + roots[n++] = (-b - q) / a; + if (q > 0) + roots[n++] = (-b + q) / a; + return n; // 1 or 2 solutions }, /** @@ -157,37 +156,43 @@ var Numerical = new function() { * a*x^3 + b*x^2 + c*x + d = 0 */ solveCubic: function(a, b, c, d, roots, tolerance) { - // After Numerical Recipes in C, 2nd edition, Press et al., - // 5.6, Quadratic and Cubic Equations + // Code ported over and adapted from Uintah library (MIT license). if (abs(a) < tolerance) return Numerical.solveQuadratic(b, c, d, roots, tolerance); - // Normalize + // Normalize to form: x^3 + b x^2 + c x + d = 0: b /= a; c /= a; d /= a; // Compute discriminants - var Q = (b * b - 3 * c) / 9, - R = (2 * b * b * b - 9 * b * c + 27 * d) / 54, - Q3 = Q * Q * Q, - R2 = R * R; - b /= 3; // Divide by 3 as that's required below - if (R2 < Q3) { // Three real roots - // This sqrt and division is safe, since R2 >= 0, so Q3 > R2, - // so Q3 > 0. The acos is also safe, since R2/Q3 < 1, and - // thus R/sqrt(Q3) < 1. - var theta = Math.acos(R / sqrt(Q3)), - // This sqrt is safe, since Q3 >= 0, and thus Q >= 0 - q = -2 * sqrt(Q); - roots[0] = q * cos(theta / 3) - b; - roots[1] = q * cos((theta + 2 * PI) / 3) - b; - roots[2] = q * cos((theta - 2 * PI) / 3) - b; - return 3; - } else { // One real root - var A = -Math.pow(abs(R) + sqrt(R2 - Q3), 1 / 3); - if (R < 0) A = -A; - var B = (abs(A) < tolerance) ? 0 : Q / A; - roots[0] = (A + B) - b; - return 1; + var bb = b * b, + p = 1 / 3 * (-1 / 3 * bb + c), + q = 1 / 2 * (2 / 27 * b * bb - 1 / 3 * b * c + d), + // Use Cardano's formula + ppp = p * p * p, + D = q * q + ppp; + // Substitute x = y - b/3 to eliminate quadric term: x^3 +px + q = 0 + b /= 3; + if (abs(D) < tolerance) { + if (abs(q) < tolerance) { // One triple solution. + roots[0] = - b; + return 1; + } else { // One single and one double solution. + var u = cbrt(-q); + roots[0] = 2 * u - b; + roots[1] = - u - b; + return 2; + } + } else if (D < 0) { // Casus irreducibilis: three real solutions + var phi = 1 / 3 * Math.acos(-q / sqrt(-ppp)); + var t = 2 * sqrt(-p); + roots[0] = t * cos(phi) - b; + roots[1] = - t * cos(phi + PI / 3) - b; + roots[2] = - t * cos(phi - PI / 3) - b; + return 3; + } else { // One real solution + D = sqrt(D); + roots[0] = cbrt(D - q) - cbrt(D + q) - b; + return 1; } } };