Update Numerical.solveQuadratic() / solveCubic() to optionally filter results to be in a given range.

This commit is contained in:
Jürg Lehni 2013-10-18 13:52:01 +02:00
parent 523b9ea592
commit eae526f38c

View file

@ -126,15 +126,22 @@ var Numerical = new function() {
*
* a*x^2 + b*x + c = 0
*/
solveQuadratic: function(a, b, c, roots) {
solveQuadratic: function(a, b, c, roots, min, max) {
var unbound = min === undefined,
count = 0;
function add(root) {
if (unbound || root >= min && root <= max)
roots[count++] = root;
return count;
}
// Code ported over and adapted from Uintah library (MIT license).
var epsilon = this.EPSILON;
// If a is 0, equation is actually linear, return 0 or 1 easy roots.
if (abs(a) < epsilon) {
if (abs(b) >= epsilon) {
roots[0] = -c / b;
return 1;
}
if (abs(b) >= epsilon)
return add(-c / b);
// If all the coefficients are 0, we have infinite solutions!
return abs(c) < epsilon ? -1 : 0; // Infinite or 0 solutions
}
@ -145,10 +152,10 @@ var Numerical = new function() {
if (p2 < q - epsilon)
return 0;
var s = p2 > q ? sqrt(p2 - q) : 0;
roots[0] = s - p;
add (s - p);
if (s > 0)
roots[1] = -s - p;
return s > 0 ? 2 : 1;
add(-s - p);
return count;
},
/**
@ -157,12 +164,22 @@ var Numerical = new function() {
*
* a*x^3 + b*x^2 + c*x + d = 0
*/
solveCubic: function(a, b, c, d, roots) {
// Code ported over and adapted from Uintah library (MIT license).
solveCubic: function(a, b, c, d, roots, min, max) {
var epsilon = this.EPSILON;
// If a is 0, equation is actually quadratic.
if (abs(a) < epsilon)
return Numerical.solveQuadratic(b, c, d, roots);
return Numerical.solveQuadratic(b, c, d, roots, min, max);
var unbound = min === undefined,
count = 0;
function add(root) {
if (unbound || root >= min && root <= max)
roots[count++] = root;
return count;
}
// Code ported over and adapted from Uintah library (MIT license).
// Normalize to form: x^3 + b x^2 + c x + d = 0:
b /= a;
c /= a;
@ -177,31 +194,26 @@ var Numerical = new function() {
// Substitute x = y - b/3 to eliminate quadric term: x^3 +px + q = 0
b /= 3;
if (abs(D) < epsilon) {
if (abs(q) < epsilon) { // One triple solution.
roots[0] = - b;
return 1;
}
if (abs(q) < epsilon) // One triple solution.
return add(-b);
// One single and one double solution.
var sqp = sqrt(p),
snq = q > 0 ? 1 : -1;
roots[0] = -snq * 2 * sqp - b;
roots[1] = snq * sqp - b;
return 2;
add(-snq * 2 * sqp - b);
return add(snq * sqp - b);
}
if (D < 0) { // Casus irreducibilis: three real solutions
var sqp = sqrt(p),
phi = Math.acos(q / (sqp * sqp * sqp)) / 3,
t = -2 * sqp,
o = 2 * PI / 3;
roots[0] = t * cos(phi) - b;
roots[1] = t * cos(phi + o) - b;
roots[2] = t * cos(phi - o) - b;
return 3;
add(t * cos(phi) - b);
add(t * cos(phi + o) - b);
return add(t * cos(phi - o) - b);
}
// One real solution
var A = (q > 0 ? -1 : 1) * pow(abs(q) + sqrt(D), 1 / 3);
roots[0] = A + p / A - b;
return 1;
return add(A + p / A - b);
}
};
};