OpenE2K
/
gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Thanks to Brian Gough <bjg@network-theory.com> 

2004-01-05  Sascha Brawer  <brawer@dandelis.ch>

	Thanks to Brian Gough <bjg@network-theory.com>
	* java/awt/geom/CubicCurve2D.java (solveCubic): Implemented.
	* java/awt/geom/QuadCurve2D.java (solveQuadratic): Re-written.

From-SVN: r75437
This commit is contained in:
Sascha Brawer 2004-01-05 20:19:29 +01:00 committed by Michael Koch
parent 60b799fd29
commit ab22bc9148
3 changed files with 300 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
libjava/ChangeLog
View File

@ -1,3 +1,9 @@
2004-01-05 Sascha Brawer <brawer@dandelis.ch>
Thanks to Brian Gough <bjg@network-theory.com>
* java/awt/geom/CubicCurve2D.java (solveCubic): Implemented.
* java/awt/geom/QuadCurve2D.java (solveQuadratic): Re-written.
2004-01-04 Matthias Klose <doko@debian.org>
* aclocal.m4: Rebuilt using "aclocal -I .".

166
libjava/java/awt/geom/CubicCurve2D.java
View File

@ -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;
}

144
libjava/java/awt/geom/QuadCurve2D.java
View File

@ -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)
{
res[0] = -b - det;
res[1] = -c / (b + det);
// 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)
{
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;
}