Solving quadratic and cubic equations

5 files changed, 310 insertions, 16 deletions

diff --git a/ChangeLog b/ChangeLog
index cc8cd247e..4171e37a1 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,12 @@
+2004-01-05  Sascha Brawer  <brawer@dandelis.ch>
+
+	Fix for Classpath bug #6095
+	Thanks to Brian Gough <bjg@network-theory.com>
+	* java/awt/geom/CubicCurve2D.java (solveCubic): Implemented.
+	* java/awt/geom/QuadCurve2D.java (solveQuadratic): Re-written.
+	* NEWS: Mention the new capability for solving equations.
+	* THANKYOU: Add Brian Gough.
+
 2004-01-04  Michael Koch  <konqueror@gmx.de>
 
 	* java/net/JarURLConnection.java
diff --git a/NEWS b/NEWS
index 951c3d0bd..41b7c1483 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,11 @@
 New since last release
+
 * Split native methods in java.lang.Runtime into java.lang.VMRuntime.
+* java.awt.geom.CubicCurve2D/QuadCurve2D: Can now solve cubic and quadratic
+  equations; implementation adapted from the GNU Scientific Library.
+* Fixed Classpath bugs:
+  #6095 java.awt.geom.QuadCurve2D.solveQuadratic sometimes gives
+        wrong results
 
 New in release 0.07 (2003/30/11)
 
diff --git a/THANKYOU b/THANKYOU
index 1fdf689e6..d191d7c05 100644
--- a/THANKYOU
+++ b/THANKYOU
@@ -18,6 +18,7 @@ Raimar Falke (hawk@hawk.shef.ac.uk)
 Philip Fong (pwlfong@users.sourceforge.net)
 Jeroen Frijters (jeroen@sumatra.nl)
 Etienne M. Gagnon (etienne.gagnon@uqam.ca)
+Brian Gough (bjg@network-theory.com)
 Fred Gray (fegray@npl.uiuc.edu)
 Christian Grothoff <grothoff@cs.purdue.edu>
 David P. Grove (groved@us.ibm.com)
diff --git a/java/awt/geom/CubicCurve2D.java b/java/awt/geom/CubicCurve2D.java
index 1e38d3ada..096e7ad97 100644
--- a/java/awt/geom/CubicCurve2D.java
+++ b/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:
+   *
+   * <blockquote><code>eqn[3]</code> &#xb7; <i>x</i><sup>3</sup>
+   * + <code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
+   * + <code>eqn[1]</code> &#xb7; <i>x</i>
+   * + <code>eqn[0]</code>
+   * = 0
+   * </blockquote>
+   *
+   * <p>For some background about solving cubic equations, see the
+   * article <a
+   * href="http://planetmath.org/encyclopedia/CubicFormula.html"
+   * >&#x201c;Cubic Formula&#x201d;</a> in <a
+   * href="http://planetmath.org/" >PlanetMath</a>.  For an extensive
+   * library of numerical algorithms written in the C programming
+   * language, see the <a href= "http://www.gnu.org/software/gsl/">GNU
+   * Scientific Library</a>, from which this implementation was
+   * adapted.
+   *
+   * @param eqn an array with the coefficients of the equation. When
+   * this procedure has returned, <code>eqn</code> 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 <a href="mailto:bjg@network-theory.com">Brian Gough</a>
+   * (original C implementation in the <a href=
+   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
+   *
+   * @author <a href="mailto:brawer@dandelis.ch">Sascha Brawer</a>
+   * (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:
+   *
+   * <blockquote><code>eqn[3]</code> &#xb7; <i>x</i><sup>3</sup>
+   * + <code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
+   * + <code>eqn[1]</code> &#xb7; <i>x</i>
+   * + <code>eqn[0]</code>
+   * = 0
+   * </blockquote>
+   *
+   * <p>For some background about solving cubic equations, see the
+   * article <a
+   * href="http://planetmath.org/encyclopedia/CubicFormula.html"
+   * >&#x201c;Cubic Formula&#x201d;</a> in <a
+   * href="http://planetmath.org/" >PlanetMath</a>.  For an extensive
+   * library of numerical algorithms written in the C programming
+   * language, see the <a href= "http://www.gnu.org/software/gsl/">GNU
+   * Scientific Library</a>, 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 <code>eqn</code>
+   * and <code>res</code>.
+   *
+   * @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 <a href="mailto:bjg@network-theory.com">Brian Gough</a>
+   * (original C implementation in the <a href=
+   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
+   *
+   * @author <a href="mailto:brawer@dandelis.ch">Sascha Brawer</a>
+   * (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/java/awt/geom/QuadCurve2D.java b/java/awt/geom/QuadCurve2D.java
index 5bc63e6c6..409c4841e 100644
--- a/java/awt/geom/QuadCurve2D.java
+++ b/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:
+   *
+   * <blockquote><code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
+   * + <code>eqn[1]</code> &#xb7; <i>x</i>
+   * + <code>eqn[0]</code>
+   * = 0
+   * </blockquote>
+   *
+   * <p>For some background about solving quadratic equations, see the
+   * article <a href=
+   * "http://planetmath.org/encyclopedia/QuadraticFormula.html"
+   * >&#x201c;Quadratic Formula&#x201d;</a> in <a href=
+   * "http://planetmath.org/">PlanetMath</a>. For an extensive library
+   * of numerical algorithms written in the C programming language,
+   * see the <a href="http://www.gnu.org/software/gsl/">GNU Scientific
+   * Library</a>.
+   *
+   * @see #solveQuadratic(double[], double[])
+   * @see CubicCurve2D#solveCubic(double[], double[])
+   *
+   * @param eqn an array with the coefficients of the equation. When
+   * this procedure has returned, <code>eqn</code> 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 <a href="mailto:bjg@network-theory.com">Brian Gough</a>
+   * (original C implementation in the <a href=
+   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
+   *
+   * @author <a href="mailto:brawer@dandelis.ch">Sascha Brawer</a>
+   * (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:
+   *
+   * <blockquote><code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
+   * + <code>eqn[1]</code> &#xb7; <i>x</i>
+   * + <code>eqn[0]</code>
+   * = 0
+   * </blockquote>
+   *
+   * <p>For some background about solving quadratic equations, see the
+   * article <a href=
+   * "http://planetmath.org/encyclopedia/QuadraticFormula.html"
+   * >&#x201c;Quadratic Formula&#x201d;</a> in <a href=
+   * "http://planetmath.org/">PlanetMath</a>. For an extensive library
+   * of numerical algorithms written in the C programming language,
+   * see the <a href="http://www.gnu.org/software/gsl/">GNU Scientific
+   * Library</a>.
+   *
+   * @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 <code>eqn</code>
+   * and <code>res</code>.
+   *
+   * @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 <a href="mailto:bjg@network-theory.com">Brian Gough</a>
+   * (original C implementation in the <a href=
+   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
+   *
+   * @author <a href="mailto:brawer@dandelis.ch">Sascha Brawer</a>
+   * (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;
   }