Name: gm110360 Date: 03/01/2002
FULL PRODUCT VERSION :
java version "1.4.0-rc"
Java(TM) 2 Runtime Environment, Standard Edition (build 1.4.
Java HotSpot(TM) Client VM (build 1.4.0-rc-b91, mixed mode)
FULL OPERATING SYSTEM VERSION :
Microsoft Windows 2000 [Version 5.00.2195]
A DESCRIPTION OF THE PROBLEM :
The method CubicCurve2D.cubicSolve() in java.awt.geom does
not find ALL real roots of a cubic equation.
For example: x^3 + x^2 = 0 has the two solutions 0 and -1,
but cubicSolve() only finds -1.
STEPS TO FOLLOW TO REPRODUCE THE PROBLEM :
1. use method CubicCurve2D.solveCubic()
(see my code example)
EXPECTED VERSUS ACTUAL BEHAVIOR :
As stated above, cubicSolve() should return ALL real
solutions in any situation.
ERROR MESSAGES/STACK TRACES THAT OCCUR :
No Errors or Exceptions are thrown.
This bug can be reproduced always.
---------- BEGIN SOURCE ----------
class SolveCubicBug {
public static void main (String args[]) {
double [] eqn = {0,0,1,1};
double [] res = new double[3];
// JDK function
int num = java.awt.geom.CubicCurve2D.solveCubic(eqn, res);
System.out.println("solveCubic() says:");
System.out.println("The cubic equation " +
eqn[3] + " x^3 + " +
eqn[2] + " x^2 + " +
eqn[1] + " x + " +
eqn[0] + " = 0 has " + num + " root(s):");
for(int i = 0; i < num; i++) {
System.out.println("x" + (i+1) + " = " + res[i]);
}
// my function
num = mySolveCubic(eqn, res);
System.out.println("\nmySolveCubic() says:");
System.out.println("The cubic equation " +
eqn[3] + " x^3 + " +
eqn[2] + " x^2 + " +
eqn[1] + " x + " +
eqn[0] + " = 0 has " + num + " root(s):");
for(int i = 0; i < num; i++) {
System.out.println("x" + (i+1) + " = " + res[i]);
}
}
public static int mySolveCubic(double eqn[], double res[]) {
// From Numerical Recipes, 5.6, Quadratic and Cubic Equations
// case discriminant == 0 added by Markus Hohenwarter, 20.1.2002
double d = eqn[3];
if (d == 0.0) {
// The cubic has degenerated to quadratic (or line or ...).
// CHANGE:
// replace "return java.awt.geom.QuadCurve2D.solveQuadratic(eqn, res);"
// by "return solveQuadratic(eqn, res);"
return java.awt.geom.QuadCurve2D.solveQuadratic(eqn, res);
}
double a = eqn[2] / d;
double b = eqn[1] / d;
double c = eqn[0] / d;
int roots = 0;
double Q = (a * a - 3.0 * b) / 9.0;
double R = (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) / 54.0;
double R2 = R * R;
double Q3 = Q * Q * Q;
double discriminant = R2 - Q3;
a = a / 3.0;
if (discriminant == 0.0) {
if (Q == 0.0) {
// one real solution
res[roots++] = -a;
} else {
// two real solutions
Q = Math.sqrt(Q);
res[roots++] = -2.0 * Q - a;
res[roots++] = Q - a;
}
} else {
if (R2 < Q3) { // => Q3 > 0.0
double theta = Math.acos(R / Math.sqrt(Q3));
Q = -2.0 * Math.sqrt(Q);
if (res == eqn) {
// Copy the eqn so that we don't clobber it with the
// roots. This is needed so that fixRoots can do its
// work with the original equation.
eqn = new double[4];
System.arraycopy(res, 0, eqn, 0, 4);
}
res[roots++] = Q * Math.cos(theta / 3.0) - a;
res[roots++] = Q * Math.cos((theta - Math.PI * 2.0)/ 3.0) - a;
res[roots++] = Q * Math.cos((theta + Math.PI * 2.0)/ 3.0) - a;
// CHANGE: "fixRoots(res, eqn);" should be invoked now
} else {
boolean neg = (R < 0.0);
double S = Math.sqrt(discriminant);
if (neg) {
R = -R;
}
double A = Math.pow(R + S, 1.0 / 3.0);
if (!neg) {
A = -A;
}
double B = (A == 0.0) ? 0.0 : (Q / A);
res[roots++] = (A + B) - a;
}
}
return roots;
}
}
---------- END SOURCE ----------
CUSTOMER WORKAROUND :
Reason of the problem: The Java implementation does not
deal with an important special case:
For R^2 == Q^3 = 0 there is only one real solution (this
one is found correctly by cubicSolve()).
For R^2 == Q^3 != 0 there are TWO different real solutions,
but cubicSolve() only finds one of them. Here the two
complex conjugate solutions become one real solution. The
Numerical Recipes (5.6) say that there is only one real
solution in the case of R^2 >= Q^3. This is wrong, as there
are two for R^2 == Q^3 != 0 because the complex conjugates
become real.
Solution to the problem:
Include the special case R^2 == Q^3. The repaired code is
in the source code field above (see function mySolveCubic
()).
(Review ID: 139073)
======================================================================