Hi Jim.
> Unfortunately, this means that the names here
> and the values assigned to them and the comment above them conflict.
> If the variables could be named "p/3" and "q/2" then all would be clear,
> but I don't know how to do that naming very easily. Perhaps the
> comment could be simply reworded:
>
> // substitute <blah blah blah>
> // x^3 + Px + Q = 0
> // Since we actually need P/3 and Q/2 for all of the
> // calculations that follow, we will calculate
> // p = P/3
> // q = Q/2
> // instead and use those values for simplicity of the code.
Good point. Done.
> Line 1105 - single quotes in comments freaks out my version of
> gnuemacs.
> I usually try to avoid them, except in pairs, but there isn't a better
> way to word this comment. :-(
We could always go with "the method of Cardano", although
that doesn't sound good at all. Or we could use a backtick instead of
the single quote - it looks similar enough.
> Lines 1157-1163 - the old code used to copy the eqn before it got
> clobbered with roots. Here it is too late. You probably need to move
> this code up near line 1135 before the 3 roots are stuffed into the
> res array. (Did you test the eqn==res case?)
You're right. I'm sorry about this.
> I noticed that the "Casus irreducibilis" case isn't in Cordano's
> method.
> He only finds roots for the 1 and 2 root case and punts for 3 roots.
> So, this is someone else's method. It would be nice to figure out who
> or what and list a reference, even though the Graphics Gems and the
> old
> code didn't. The closest reference I can find is unattributed on
> Wikipedia, but you could include it in a comment for reference:
>
> http://en.wikipedia.org/wiki/Cubic_function#Trigonometric_.28and_hyperbolic.29_method
Done.
> Line 1133 - I don't understand why that term has -q in it. The above
> link and the original code both computed essentially the arccos of
> this
> formula without the negation of q. ??? Since acos(-v) == pi - acos(v)
> this would seem to negate the result and bias it by pi/3. Negating it
> won't affect the eventual cosine, but the bias by pi/3 will. Am I
> missing something?
I think you are. What's going on is that our code is computing
ret[0] = t2, ret[1] = t0, ret[2] = t1, where (t0, t1, t2 are the tk's
from the wikipedia link).
Let X = (3q/2p)*sqrt(-3/p) where p and q are the ones from the wikipedia
article, not our code.
So, when we do ret[0] = t * cos(phi), that's:
= t * cos(acos(-X))
= t * cos(pi/3 - acos(X))
= t * cos(acos(X) - pi/3)
= t * cos(acos(X) - pi/3 - pi)
= t * cos(acos(X) - 2*2*pi/3)
= t2
ret[0] and ret[1] are very similar to prove - you just add/subtract pi/3
from the argument to cos.
I unfortunately don't have access to the icedtea servers at this moment,
so I attached a patch. I hope that's ok.
Regards,
Denis.
--- dashing/2d/jdk/src/share/classes/java/awt/geom/CubicCurve2D.java 2010-12-24 10:01:45.378556843 -0500
+++ cc2d/2d/jdk/src/share/classes/java/awt/geom/CubicCurve2D.java 2010-12-24 11:01:42.205614715 -0500
@@ -1083,24 +1083,63 @@
* @since 1.3
*/
public static int solveCubic(double eqn[], double res[]) {
- // From Numerical Recipes, 5.6, Quadratic and Cubic Equations
- double d = eqn[3];
- if (d == 0.0) {
- // The cubic has degenerated to quadratic (or line or ...).
+ // From Graphics Gems:
+ // http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
+ final double d = eqn[3];
+ if (d == 0) {
return 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;
- a = a / 3.0;
- if (R2 < Q3) {
- double theta = Math.acos(R / Math.sqrt(Q3));
- Q = -2.0 * Math.sqrt(Q);
+ /* normal form: x^3 + Ax^2 + Bx + C = 0 */
+ final double A = eqn[2] / d;
+ final double B = eqn[1] / d;
+ final double C = eqn[0] / d;
+
+
+ // substitute x = y - A/3 to eliminate quadratic term:
+ // y^3 +Py + Q = 0
+ //
+ // Since we actually need P/3 and Q/2 for all of the
+ // calculations that follow, we will calculate
+ // p = P/3
+ // q = Q/2
+ // instead and use those values for simplicity of the code.
+
+ double sq_A = A * A;
+ double p = 1.0/3 * (-1.0/3 * sq_A + B);
+ double q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
+
+ /* use Cardano's formula */
+
+ double cb_p = p * p * p;
+ double D = q * q + cb_p;
+
+ int num;
+ // XXX: we consider 0 to be anything within 1e-9 of 0.
+ // Is this really right? Maybe we should use a bound that changes
+ // with the number being tested (math.ulp would do that, but
+ // what input do we give it? And do we scale it, or will
+ // within(D, 0, math.ulp(somenumber)) be good enough).
+ if (iszero(D)) {
+ // XXX: do we really need iszero for q? All we do with it is
+ // take it's cube root, which works fine even for Double.MIN_VALUE
+ // Then again, if we remove it, we will get two extremely close
+ // roots for equations where there is only one zero root.
+ // We probably should use something like iszero, but much more
+ // scrict - so, within(q, 0, 2*Math.ulp(0));
+ if (iszero(q)) { /* one triple solution */
+ res[ 0 ] = 0;
+ num = 1;
+ } else { /* one single and one double solution */
+ final double u = Math.cbrt(-q);
+ res[ 0 ] = 2*u;
+ res[ 1 ] = -u;
+ num = 2;
+ }
+ } else if (D < 0) { /* Casus irreducibilis: three real solutions */
+ // see: http://en.wikipedia.org/wiki/Cubic_function#Trigonometric_.28and_hyperbolic.29_method
+ final double phi = 1.0/3 * Math.acos(-q / Math.sqrt(-cb_p));
+ final double t = 2 * Math.sqrt(-p);
+
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
@@ -1108,24 +1147,41 @@
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;
+
+ res[ 0 ] = ( t * Math.cos(phi));
+ res[ 1 ] = (-t * Math.cos(phi + Math.PI / 3));
+ res[ 2 ] = (-t * Math.cos(phi - Math.PI / 3));
+ num = 3;
+ } else { /* one real solution */
+ final double sqrt_D = Math.sqrt(D);
+ final double u = Math.cbrt(sqrt_D - q);
+ final double v = - Math.cbrt(sqrt_D + q);
+
+ res[ 0 ] = u + v;
+ num = 1;
+ }
+
+ /* resubstitute */
+
+ final double sub = 1.0/3 * A;
+
+ for (int i = 0; i < num; ++i)
+ res[ i ] -= sub;
+
+ if (num == 3) {
fixRoots(res, eqn);
- } else {
- boolean neg = (R < 0.0);
- double S = Math.sqrt(R2 - Q3);
- 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;
+
+ return num;
+ }
+
+ private static boolean iszero(double x) {
+ return within(x, 0, 1e-9);
+ }
+
+ private static boolean within(final double x, final double y, final double err) {
+ final double d = y - x;
+ return (d <= err && d >= -err);
}
/*
