From ab22bc9148e058a649bc50cb0bdc160a9ead763b Mon Sep 17 00:00:00 2001 From: Sascha Brawer Date: Mon, 5 Jan 2004 20:19:29 +0100 Subject: [PATCH] Thanks to Brian Gough 2004-01-05 Sascha Brawer Thanks to Brian Gough * java/awt/geom/CubicCurve2D.java (solveCubic): Implemented. * java/awt/geom/QuadCurve2D.java (solveQuadratic): Re-written. From-SVN: r75437 --- libjava/ChangeLog | 6 ++ libjava/java/awt/geom/CubicCurve2D.java | 166 +++++++++++++++++++++++++++++++- libjava/java/awt/geom/QuadCurve2D.java | 144 ++++++++++++++++++++++++--- 3 files changed, 300 insertions(+), 16 deletions(-) diff --git a/libjava/ChangeLog b/libjava/ChangeLog index b1c9333..f8eadf9 100644 --- a/libjava/ChangeLog +++ b/libjava/ChangeLog @@ -1,3 +1,9 @@ +2004-01-05 Sascha Brawer + + Thanks to Brian Gough + * java/awt/geom/CubicCurve2D.java (solveCubic): Implemented. + * java/awt/geom/QuadCurve2D.java (solveQuadratic): Re-written. + 2004-01-04 Matthias Klose * aclocal.m4: Rebuilt using "aclocal -I .". diff --git a/libjava/java/awt/geom/CubicCurve2D.java b/libjava/java/awt/geom/CubicCurve2D.java index 1e38d3a..096e7ad 100644 --- a/libjava/java/awt/geom/CubicCurve2D.java +++ b/libjava/java/awt/geom/CubicCurve2D.java @@ -624,17 +624,115 @@ public abstract class CubicCurve2D } + /** + * Finds the non-complex roots of a cubic equation, placing the + * results into the same array as the equation coefficients. The + * following equation is being solved: + * + *
eqn[3] · x3 + * + eqn[2] · x2 + * + eqn[1] · x + * + eqn[0] + * = 0 + *
+ * + *

For some background about solving cubic equations, see the + * article “Cubic Formula” in PlanetMath. For an extensive + * library of numerical algorithms written in the C programming + * language, see the GNU + * Scientific Library, from which this implementation was + * adapted. + * + * @param eqn an array with the coefficients of the equation. When + * this procedure has returned, eqn will contain the + * non-complex solutions of the equation, in no particular order. + * + * @return the number of non-complex solutions. A result of 0 + * indicates that the equation has no non-complex solutions. A + * result of -1 indicates that the equation is constant (i.e., + * always or never zero). + * + * @see #solveCubic(double[], double[]) + * @see QuadCurve2D#solveQuadratic(double[],double[]) + * + * @author Brian Gough + * (original C implementation in the GNU Scientific Library) + * + * @author Sascha Brawer + * (adaptation to Java) + */ public static int solveCubic(double[] eqn) { return solveCubic(eqn, eqn); } + /** + * Finds the non-complex roots of a cubic equation. The following + * equation is being solved: + * + *

eqn[3] · x3 + * + eqn[2] · x2 + * + eqn[1] · x + * + eqn[0] + * = 0 + *
+ * + *

For some background about solving cubic equations, see the + * article “Cubic Formula” in PlanetMath. For an extensive + * library of numerical algorithms written in the C programming + * language, see the GNU + * Scientific Library, from which this implementation was + * adapted. + * + * @see QuadCurve2D#solveQuadratic(double[],double[]) + * + * @param eqn an array with the coefficients of the equation. + * + * @param res an array into which the non-complex roots will be + * stored. The results may be in an arbitrary order. It is safe to + * pass the same array object reference for both eqn + * and res. + * + * @return the number of non-complex solutions. A result of 0 + * indicates that the equation has no non-complex solutions. A + * result of -1 indicates that the equation is constant (i.e., + * always or never zero). + * + * @author Brian Gough + * (original C implementation in the GNU Scientific Library) + * + * @author Sascha Brawer + * (adaptation to Java) + */ public static int solveCubic(double[] eqn, double[] res) { + // Adapted from poly/solve_cubic.c in the GNU Scientific Library + // (GSL), revision 1.7 of 2003-07-26. For the original source, see + // http://www.gnu.org/software/gsl/ + // + // Brian Gough, the author of that code, has granted the + // permission to use it in GNU Classpath under the GNU Classpath + // license, and has assigned the copyright to the Free Software + // Foundation. + // + // The Java implementation is very similar to the GSL code, but + // not a strict one-to-one copy. For example, GSL would sort the + // result. + double a, b, c, q, r, Q, R; - - double c3 = eqn[3]; + double c3, Q3, R2, CR2, CQ3; + + // If the cubic coefficient is zero, we have a quadratic equation. + c3 = eqn[3]; if (c3 == 0) return QuadCurve2D.solveQuadratic(eqn, res); @@ -644,7 +742,69 @@ public abstract class CubicCurve2D a = eqn[2] / c3; // We now need to solve x^3 + ax^2 + bx + c = 0. - throw new Error("not implemented"); // FIXME + q = a * a - 3 * b; + r = 2 * a * a * a - 9 * a * b + 27 * c; + + Q = q / 9; + R = r / 54; + + Q3 = Q * Q * Q; + R2 = R * R; + + CR2 = 729 * r * r; + CQ3 = 2916 * q * q * q; + + if (R == 0 && Q == 0) + { + // The GNU Scientific Library would return three identical + // solutions in this case. + res[0] = -a/3; + return 1; + } + + if (CR2 == CQ3) + { + /* this test is actually R2 == Q3, written in a form suitable + for exact computation with integers */ + + /* Due to finite precision some double roots may be missed, and + considered to be a pair of complex roots z = x +/- epsilon i + close to the real axis. */ + + double sqrtQ = Math.sqrt(Q); + + if (R > 0) + { + res[0] = -2 * sqrtQ - a/3; + res[1] = sqrtQ - a/3; + } + else + { + res[0] = -sqrtQ - a/3; + res[1] = 2 * sqrtQ - a/3; + } + return 2; + } + + if (CR2 < CQ3) /* equivalent to R2 < Q3 */ + { + double sqrtQ = Math.sqrt(Q); + double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; + double theta = Math.acos(R / sqrtQ3); + double norm = -2 * sqrtQ; + res[0] = norm * Math.cos(theta / 3) - a / 3; + res[1] = norm * Math.cos((theta + 2.0 * Math.PI) / 3) - a/3; + res[2] = norm * Math.cos((theta - 2.0 * Math.PI) / 3) - a/3; + + // The GNU Scientific Library sorts the results. We don't. + return 3; + } + + double sgnR = (R >= 0 ? 1 : -1); + double A = -sgnR * Math.pow(Math.abs(R) + Math.sqrt(R2 - Q3), 1.0/3.0); + double B = Q / A ; + res[0] = A + B - a/3; + return 1; } diff --git a/libjava/java/awt/geom/QuadCurve2D.java b/libjava/java/awt/geom/QuadCurve2D.java index 5bc63e6..409c484 100644 --- a/libjava/java/awt/geom/QuadCurve2D.java +++ b/libjava/java/awt/geom/QuadCurve2D.java @@ -550,39 +550,157 @@ public abstract class QuadCurve2D } + /** + * Finds the non-complex roots of a quadratic equation, placing the + * results into the same array as the equation coefficients. The + * following equation is being solved: + * + *

eqn[2] · x2 + * + eqn[1] · x + * + eqn[0] + * = 0 + *
+ * + *

For some background about solving quadratic equations, see the + * article “Quadratic Formula” in PlanetMath. For an extensive library + * of numerical algorithms written in the C programming language, + * see the GNU Scientific + * Library. + * + * @see #solveQuadratic(double[], double[]) + * @see CubicCurve2D#solveCubic(double[], double[]) + * + * @param eqn an array with the coefficients of the equation. When + * this procedure has returned, eqn will contain the + * non-complex solutions of the equation, in no particular order. + * + * @return the number of non-complex solutions. A result of 0 + * indicates that the equation has no non-complex solutions. A + * result of -1 indicates that the equation is constant (i.e., + * always or never zero). + * + * @author Brian Gough + * (original C implementation in the GNU Scientific Library) + * + * @author Sascha Brawer + * (adaptation to Java) + */ public static int solveQuadratic(double[] eqn) { return solveQuadratic(eqn, eqn); } + /** + * Finds the non-complex roots of a quadratic equation. The + * following equation is being solved: + * + *

eqn[2] · x2 + * + eqn[1] · x + * + eqn[0] + * = 0 + *
+ * + *

For some background about solving quadratic equations, see the + * article “Quadratic Formula” in PlanetMath. For an extensive library + * of numerical algorithms written in the C programming language, + * see the GNU Scientific + * Library. + * + * @see CubicCurve2D#solveCubic(double[],double[]) + * + * @param eqn an array with the coefficients of the equation. + * + * @param res an array into which the non-complex roots will be + * stored. The results may be in an arbitrary order. It is safe to + * pass the same array object reference for both eqn + * and res. + * + * @return the number of non-complex solutions. A result of 0 + * indicates that the equation has no non-complex solutions. A + * result of -1 indicates that the equation is constant (i.e., + * always or never zero). + * + * @author Brian Gough + * (original C implementation in the GNU Scientific Library) + * + * @author Sascha Brawer + * (adaptation to Java) + */ public static int solveQuadratic(double[] eqn, double[] res) { - double c = eqn[0]; - double b = eqn[1]; - double a = eqn[2]; + // Taken from poly/solve_quadratic.c in the GNU Scientific Library + // (GSL), cvs revision 1.7 of 2003-07-26. For the original source, + // see http://www.gnu.org/software/gsl/ + // + // Brian Gough, the author of that code, has granted the + // permission to use it in GNU Classpath under the GNU Classpath + // license, and has assigned the copyright to the Free Software + // Foundation. + // + // The Java implementation is very similar to the GSL code, but + // not a strict one-to-one copy. For example, GSL would sort the + // result. + + double a, b, c, disc; + + c = eqn[0]; + b = eqn[1]; + a = eqn[2]; + + // Check for linear or constant functions. This is not done by the + // GNU Scientific Library. Without this special check, we + // wouldn't return -1 for constant functions, and 2 instead of 1 + // for linear functions. if (a == 0) { if (b == 0) return -1; + res[0] = -c / b; return 1; } - c /= a; - b /= a * 2; - double det = Math.sqrt(b * b - c); - if (det != det) + + disc = b * b - 4 * a * c; + + if (disc < 0) return 0; - // For fewer rounding errors, we calculate the two roots differently. - if (b > 0) + + if (disc == 0) + { + // The GNU Scientific Library returns two identical results here. + // We just return one. + res[0] = -0.5 * b / a ; + return 1; + } + + // disc > 0 + if (b == 0) { - res[0] = -b - det; - res[1] = -c / (b + det); + double r; + + r = Math.abs(0.5 * Math.sqrt(disc) / a); + res[0] = -r; + res[1] = r; } else { - res[0] = -c / (b - det); - res[1] = -b + det; + double sgnb, temp; + + sgnb = (b > 0 ? 1 : -1); + temp = -0.5 * (b + sgnb * Math.sqrt(disc)); + + // The GNU Scientific Library sorts the result here. We don't. + res[0] = temp / a; + res[1] = c / temp; } return 2; } -- 2.7.4