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
> // 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
