From 02658c9e743685973b6779d60116adc28008df04 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrg=20Lehni?= Date: Sat, 9 Jul 2016 01:01:19 +0200 Subject: [PATCH] Clean-up code from PR #1087 Closes #1085 --- src/util/Numerical.js | 93 ++++++++++++++++++++++------------------- test/tests/Numerical.js | 45 ++++++++++++++++++++ test/tests/load.js | 2 + 3 files changed, 96 insertions(+), 44 deletions(-) create mode 100644 test/tests/Numerical.js diff --git a/src/util/Numerical.js b/src/util/Numerical.js index 9c339b89..d08a93dc 100644 --- a/src/util/Numerical.js +++ b/src/util/Numerical.js @@ -60,6 +60,7 @@ var Numerical = new function() { var abs = Math.abs, sqrt = Math.sqrt, pow = Math.pow, + // Constants EPSILON = 1e-12, MACHINE_EPSILON = 1.12e-16; @@ -67,23 +68,37 @@ 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 getDiscriminant(a, b, c) { + // Ported from @hkrish's polysolve.c + function split(a) { + var x = a * 134217729, + y = a - x, + hi = y + x, // Don't optimize y away! + lo = a - hi; + return [hi, lo]; + } + + var D = b * b - a * c, + E = b * b + a * c; + if (abs(D) * 3 < E) { + var ad = split(a), + bd = split(b), + cd = split(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]; + D = (p - q) + (dp - dq); // Don’t omit parentheses! + } + return D; } - 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); + function getNormalizationFactor(x) { + // 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.floor( + Math.log(x || MACHINE_EPSILON) * Math.LOG2E + 0.5)); } return /** @lends Numerical */{ @@ -241,31 +256,21 @@ var Numerical = new function() { eMax = max + EPSILON, x1, x2 = Infinity, B = b, - D, E, pi = 3; + D, E; // 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); - } + 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 !== 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.LOG2E + 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 (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); } if (abs(a) < EPSILON) { // This could just be a linear equation @@ -326,18 +331,18 @@ var Numerical = new function() { * @author Harikrishnan Gopalakrishnan */ solveCubic: function(a, b, c, d, roots, min, max) { - var count = 0, x, b1, c2, + var 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) { + // 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 p = pow(2, -Math.floor(Math.log(s) * Math.LOG2E)); - a *= p; - b *= p; - c *= p; - d *= p; + var f = getNormalizationFactor(s); + a *= f; + b *= f; + c *= f; + d *= f; } // If a or d is zero, we only need to solve a quadratic, so we set // the coefficients appropriately. diff --git a/test/tests/Numerical.js b/test/tests/Numerical.js new file mode 100644 index 00000000..01a5cb08 --- /dev/null +++ b/test/tests/Numerical.js @@ -0,0 +1,45 @@ +/* + * Paper.js - The Swiss Army Knife of Vector Graphics Scripting. + * http://paperjs.org/ + * + * Copyright (c) 2011 - 2016, Juerg Lehni & Jonathan Puckey + * http://scratchdisk.com/ & http://jonathanpuckey.com/ + * + * Distributed under the MIT license. See LICENSE file for details. + * + * All rights reserved. + */ + +QUnit.module('Numerical'); + +test('Numerical.solveQuadratic()', function() { + function solve(s) { + var roots = [], + count = Numerical.solveQuadratic(s, 0, -s, roots); + return roots; + } + + var expected = [1, -1]; + + equals(solve(1), expected, + 'Numerical.solveQuadratic().'); + equals(solve(Numerical.EPSILON), expected, + 'Numerical.solveQuadratic() with an identical set of' + + 'coefficients at different scale.'); +}); + +test('Numerical.solveCubic()', function() { + function solve(s) { + var roots = [], + count = Numerical.solveCubic(0.5 * s, -s, -s, -s, roots); + return roots; + } + + var expected = [2.919639565839418]; + + equals(solve(1), expected, + 'Numerical.solveCubic().'); + equals(solve(Numerical.EPSILON), expected, + 'Numerical.solveCubic() with an identical set of' + + 'coefficients at different scale.'); +}); diff --git a/test/tests/load.js b/test/tests/load.js index 6514919d..c1f9f6aa 100644 --- a/test/tests/load.js +++ b/test/tests/load.js @@ -58,3 +58,5 @@ /*#*/ include('SvgImport.js'); /*#*/ include('SvgExport.js'); + +/*#*/ include('Numerical.js');