From 9d6aab3802082439f0fe933a22fddf14dfe34173 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrg=20Lehni?= Date: Sat, 9 Jul 2016 13:28:50 +0200 Subject: [PATCH] Streamline handling of getNormalizationFactor() to share more code. Based on comments by @hkrish in https://github.com/paperjs/paper.js/pull/1087#issuecomment-231529361 --- src/util/Numerical.js | 103 ++++++++++++++++++++++-------------------- 1 file changed, 53 insertions(+), 50 deletions(-) diff --git a/src/util/Numerical.js b/src/util/Numerical.js index a621b134..7d3d0a5d 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -74,11 +74,11 @@ var Numerical = new function() { function getDiscriminant(a, b, c) { // Ported from @hkrish's polysolve.c - function split(a) { - var x = a * 134217729, - y = a - x, + function split(v) { + var x = v * 134217729, + y = v - x, hi = y + x, // Don't optimize y away! - lo = a - hi; + lo = v - hi; return [hi, lo]; } @@ -98,10 +98,14 @@ var Numerical = new function() { return D; } - function getNormalizationFactor(x) { + function getNormalizationFactor() { + var max = Math.max.apply(Math, arguments); + // Normalize coefficients à la Jenkins & Traub's RPOLY. // Normalization is done by scaling coefficients with a power of 2, so // that all the bits in the mantissa remain unchanged. - return pow(2, -Math.round(log2(x || MACHINE_EPSILON))); + return max && (max < 1e-8 || max > 1e8) + ? pow(2, -Math.round(log2(max))) + : 0; } return /** @lends Numerical */{ @@ -254,49 +258,52 @@ var Numerical = new function() { * @author Harikrishnan Gopalakrishnan */ solveQuadratic: function(a, b, c, roots, min, max) { - var count = 0, - eMin = min - EPSILON, - eMax = max + EPSILON, - x1, x2 = Infinity, - // a, b, c are expected to be the coefficients of the equation: - // Ax² - 2Bx + C == 0, so we take B = -b/2: - B = b * -0.5, - D = getDiscriminant(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 && abs(D) < MACHINE_EPSILON) { - // Normalize coefficients. - var f = getNormalizationFactor(abs(a) + abs(B) + abs(c)); - a *= f; - b *= f; - c *= f; - B *= f; - D = getDiscriminant(a, B, c); - } + var x1, x2 = Infinity; if (abs(a) < EPSILON) { // This could just be a linear equation if (abs(b) < EPSILON) return abs(c) < EPSILON ? -1 : 0; x1 = -c / b; - } else if (D >= -MACHINE_EPSILON) { // No real roots if D < 0 - var Q = D < 0 ? 0 : sqrt(D), - R = B + (B < 0 ? -Q : Q); - // Try to minimize floating point noise. - if (R === 0) { - x1 = c / a; - x2 = -x1; - } else { - x1 = R / a; - x2 = c / R; + } else { + // a, b, c are expected to be the coefficients of the equation: + // Ax² - 2Bx + C == 0, so we take b = -b/2: + b *= -0.5; + var D = getDiscriminant(a, b, c); + // If the discriminant is very small, we can try to normalize + // the coefficients, so that we may get better accuracy. + if (D && abs(D) < MACHINE_EPSILON) { + var f = getNormalizationFactor(abs(a), abs(b), abs(c)); + if (f) { + a *= f; + b *= f; + c *= f; + D = getDiscriminant(a, b, c); + } + } + if (D >= -MACHINE_EPSILON) { // No real roots if D < 0 + var Q = D < 0 ? 0 : sqrt(D), + R = b + (b < 0 ? -Q : Q); + // Try to minimize floating point noise. + if (R === 0) { + x1 = c / a; + x2 = -x1; + } else { + x1 = R / a; + x2 = c / R; + } } } + var count = 0, + boundless = min == null, + minB = min - EPSILON, + maxB = max + EPSILON; // We need to include EPSILON in the comparisons with min / max, // as some solutions are ever so lightly out of bounds. - if (isFinite(x1) && (min == null || x1 > eMin && x1 < eMax)) - roots[count++] = min == null ? x1 : clamp(x1, min, max); + if (isFinite(x1) && (boundless || x1 > minB && x1 < maxB)) + roots[count++] = boundless ? x1 : clamp(x1, min, max); if (x2 !== x1 - && isFinite(x2) && (min == null || x2 > eMin && x2 < eMax)) - roots[count++] = min == null ? x2 : clamp(x2, min, max); + && isFinite(x2) && (boundless || x2 > minB && x2 < maxB)) + roots[count++] = boundless ? x2 : clamp(x2, min, max); return count; }, @@ -332,14 +339,9 @@ var Numerical = new function() { * @author Harikrishnan Gopalakrishnan */ solveCubic: function(a, b, c, d, roots, min, max) { - var x, b1, c2, - s = Math.max(abs(a), abs(b), abs(c), abs(d)); - // Normalize coefficients à la Jenkins & Traub's RPOLY - if (s < 1e-8 || s > 1e8) { - // Scale the coefficients by a multiple of the exponent of - // coefficients so that all the bits in the mantissa are - // preserved. - var f = getNormalizationFactor(s); + var f = getNormalizationFactor(abs(a), abs(b), abs(c), abs(d)), + x, b1, c2; + if (f) { a *= f; b *= f; c *= f; @@ -398,11 +400,12 @@ var Numerical = new function() { } } // The cubic has been deflated to a quadratic. - var count = Numerical.solveQuadratic(a, b1, c2, roots, min, max); + var count = Numerical.solveQuadratic(a, b1, c2, roots, min, max), + boundless = min == null; 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); + && (boundless || x > min - EPSILON && x < max + EPSILON)) + roots[count++] = boundless ? x : clamp(x, min, max); return count; } };