Streamline handling of getNormalizationFactor() to share more code.

Based on comments by @hkrish in https://github.com/paperjs/paper.js/pull/1087#issuecomment-231529361
This commit is contained in:
Jürg Lehni 2016-07-09 13:28:50 +02:00
parent 2532f205a7
commit 9d6aab3802

View file

@ -74,11 +74,11 @@ var Numerical = new function() {
function getDiscriminant(a, b, c) { function getDiscriminant(a, b, c) {
// Ported from @hkrish's polysolve.c // Ported from @hkrish's polysolve.c
function split(a) { function split(v) {
var x = a * 134217729, var x = v * 134217729,
y = a - x, y = v - x,
hi = y + x, // Don't optimize y away! hi = y + x, // Don't optimize y away!
lo = a - hi; lo = v - hi;
return [hi, lo]; return [hi, lo];
} }
@ -98,10 +98,14 @@ var Numerical = new function() {
return D; 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 // Normalization is done by scaling coefficients with a power of 2, so
// that all the bits in the mantissa remain unchanged. // 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 */{ return /** @lends Numerical */{
@ -254,49 +258,52 @@ var Numerical = new function() {
* @author Harikrishnan Gopalakrishnan <hari.exeption@gmail.com> * @author Harikrishnan Gopalakrishnan <hari.exeption@gmail.com>
*/ */
solveQuadratic: function(a, b, c, roots, min, max) { solveQuadratic: function(a, b, c, roots, min, max) {
var count = 0, var x1, x2 = Infinity;
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);
}
if (abs(a) < EPSILON) { if (abs(a) < EPSILON) {
// This could just be a linear equation // This could just be a linear equation
if (abs(b) < EPSILON) if (abs(b) < EPSILON)
return abs(c) < EPSILON ? -1 : 0; return abs(c) < EPSILON ? -1 : 0;
x1 = -c / b; x1 = -c / b;
} else if (D >= -MACHINE_EPSILON) { // No real roots if D < 0 } else {
var Q = D < 0 ? 0 : sqrt(D), // a, b, c are expected to be the coefficients of the equation:
R = B + (B < 0 ? -Q : Q); // Ax² - 2Bx + C == 0, so we take b = -b/2:
// Try to minimize floating point noise. b *= -0.5;
if (R === 0) { var D = getDiscriminant(a, b, c);
x1 = c / a; // If the discriminant is very small, we can try to normalize
x2 = -x1; // the coefficients, so that we may get better accuracy.
} else { if (D && abs(D) < MACHINE_EPSILON) {
x1 = R / a; var f = getNormalizationFactor(abs(a), abs(b), abs(c));
x2 = c / R; 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, // We need to include EPSILON in the comparisons with min / max,
// as some solutions are ever so lightly out of bounds. // as some solutions are ever so lightly out of bounds.
if (isFinite(x1) && (min == null || x1 > eMin && x1 < eMax)) if (isFinite(x1) && (boundless || x1 > minB && x1 < maxB))
roots[count++] = min == null ? x1 : clamp(x1, min, max); roots[count++] = boundless ? x1 : clamp(x1, min, max);
if (x2 !== x1 if (x2 !== x1
&& isFinite(x2) && (min == null || x2 > eMin && x2 < eMax)) && isFinite(x2) && (boundless || x2 > minB && x2 < maxB))
roots[count++] = min == null ? x2 : clamp(x2, min, max); roots[count++] = boundless ? x2 : clamp(x2, min, max);
return count; return count;
}, },
@ -332,14 +339,9 @@ var Numerical = new function() {
* @author Harikrishnan Gopalakrishnan <hari.exeption@gmail.com> * @author Harikrishnan Gopalakrishnan <hari.exeption@gmail.com>
*/ */
solveCubic: function(a, b, c, d, roots, min, max) { solveCubic: function(a, b, c, d, roots, min, max) {
var x, b1, c2, var f = getNormalizationFactor(abs(a), abs(b), abs(c), abs(d)),
s = Math.max(abs(a), abs(b), abs(c), abs(d)); x, b1, c2;
// Normalize coefficients à la Jenkins & Traub's RPOLY if (f) {
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);
a *= f; a *= f;
b *= f; b *= f;
c *= f; c *= f;
@ -398,11 +400,12 @@ var Numerical = new function() {
} }
} }
// The cubic has been deflated to a quadratic. // 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 if (isFinite(x) && (count === 0
|| count > 0 && x !== roots[0] && x !== roots[1]) || count > 0 && x !== roots[0] && x !== roots[1])
&& (min == null || x > min - EPSILON && x < max + EPSILON)) && (boundless || x > min - EPSILON && x < max + EPSILON))
roots[count++] = min == null ? x : clamp(x, min, max); roots[count++] = boundless ? x : clamp(x, min, max);
return count; return count;
} }
}; };