/*
GeoGebra - Dynamic Mathematics for Everyone
http://www.geogebra.org
This file is part of GeoGebra.
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation.
*/
package org.geogebra.common.kernel;
import java.util.Arrays;
import java.util.Comparator;
import java.util.Iterator;
import java.util.TreeSet;
import org.apache.commons.math3.analysis.solvers.BrentSolver;
import org.apache.commons.math3.analysis.solvers.LaguerreSolver;
import org.apache.commons.math3.analysis.solvers.NewtonSolver;
import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
import org.apache.commons.math3.complex.Complex;
import org.geogebra.common.kernel.algos.AlgoRootNewton;
import org.geogebra.common.kernel.arithmetic.MyDouble;
import org.geogebra.common.kernel.arithmetic.PolyFunction;
import org.geogebra.common.util.debug.Log;
/**
* Class for solving polynomial equations
*/
@SuppressWarnings("deprecation")
public class EquationSolver implements EquationSolverInterface {
private static final double LAGUERRE_START = -1;
private static final int LAGUERRE_MAX_EVAL = 100;
private static final double LAGUERRE_EPS = 1E-5;
private LaguerreSolver laguerreSolver;
private UnivariateSolver rootFinderBrent;
private NewtonSolver rootFinderNewton;
/**
* Creates new equation solver
*
*/
public EquationSolver() {
// we need someone to polish our roots
// rootPolisher = new RealRoot();
// extrFinder = kernel.getExtremumFinder();
}
@Override
final public int polynomialRoots(double[] roots, boolean multiple) {
int realRoots;
/*
* update: if roots[n], roots[n-1], ..., we need to determine "real"
* degree Darko Drakulic, 28.6.2011.
*/
int degree = roots.length - 1;
for (int i = degree; i >= 0 && roots[i] == 0; i--) {
degree--;
}
switch (degree) { // degree of polynomial
case 0:
realRoots = 0;
break;
case 1:
roots[0] = -roots[0] / roots[1];
realRoots = 1;
break;
case 2:
realRoots = solveQuadratic(roots, roots, Kernel.STANDARD_PRECISION);
if (multiple && realRoots == 1) {
realRoots = 2;
roots[1] = roots[0];
}
break;
case 3:
realRoots = solveCubic(roots, roots, Kernel.STANDARD_PRECISION);
break;
default:
realRoots = laguerreAll(roots);
}
// solveQuadratic and solveCubic may return -1
return Math.max(0, realRoots);
}
@Override
final public int polynomialComplexRoots(double[] real, double complex[]) {
int ret = -1;
switch (real.length - 1) { // degree of polynomial
case 0:
ret = 0;
break;
case 1:
real[0] = -real[0] / real[1];
complex[0] = 0;
ret = 1;
break;
case 2:
ret = solveQuadraticComplex(real, complex);
break;
default:
ret = laguerreAllComplex(real, complex);
}
return ret;
}
/**
* Solves the quadratic whose coefficients are in the <code>eqn</code> array
* and places the non-complex roots back into the same array, returning the
* number of roots. The quadratic solved is represented by the equation:
*
* <pre>
* eqn = {C, B, A};
* ax^2 + bx + c = 0
* </pre>
*
* A return value of <code>-1</code> is used to distinguish a constant
* equation, which might be always 0 or never 0, from an equation that has
* no zeroes.
*
* @param eqn
* the array that contains the quadratic coefficients
* @return the number of roots, or <code>-1</code> if the equation is a
* constant
*/
final public int solveQuadratic(double eqn[]) {
return solveQuadratic(eqn, eqn, Kernel.STANDARD_PRECISION);
}
@Override
final public int solveQuadratic(double eqn[], double res[], double eps) {
return solveQuadraticS(eqn, res, eps);
}
final public static int solveQuadraticS(double eqn[], double res[],
double eps) {
double a = eqn[2];
double b = eqn[1];
double c = eqn[0];
int roots = 0;
if (Math.abs(a) < eps) {
// The quadratic parabola has degenerated to a line.
if (Math.abs(b) < eps) {
// The line has degenerated to a constant.
return -1;
}
res[roots++] = -c / b;
} else if (Math.abs(b) < eps * Math.abs(a)) { // a*x^2 + c = 0
double x2 = -c / a;
if (Kernel.isZero(x2, eps)) {
res[roots++] = 0;
} else if (x2 < 0) { // no roots
return 0;
} else {
res[roots++] = Math.sqrt(x2);
res[roots++] = -Math.sqrt(x2);
}
} else {
// From Numerical Recipes, 5.6, Quadratic and Cubic Equations
double d = b * b - 4.0 * a * c;
if (Math.abs(d) < eps * b * b) {
res[roots++] = -b / (2.0 * a);
} else {
if (d < 0.0) {
// If d < 0.0, then there are no roots
return 0;
}
d = Math.sqrt(d);
// For accuracy, calculate one root using:
// (-b +/- d) / 2a
// and the other using:
// 2c / (-b +/- d)
// Choose the sign of the +/- so that b+d gets larger in
// magnitude
if (b < 0.0) {
d = -d;
}
double q = (b + d) / -2.0;
// We already tested a for being 0 above
res[roots++] = q / a;
res[roots++] = c / q;
}
}
return roots;
}
/**
* @param real
* put coefficients here on input, receive real parts on output
* @param complex
* imaginary parts on output
* @return number of roots
*/
final public static int solveQuadraticComplex(double real[],
double complex[]) {
double a = real[2];
double b = real[1];
double c = real[0];
int roots = 0;
if (Math.abs(a) < Kernel.STANDARD_PRECISION) {
// The quadratic parabola has degenerated to a line.
if (Math.abs(b) < Kernel.STANDARD_PRECISION) {
// The line has degenerated to a constant.
return -1;
}
complex[roots] = 0;
real[roots++] = -c / b;
} else {
// From Numerical Recipes, 5.6, Quadratic and Cubic Equations
double d = b * b - 4.0 * a * c;
if (Math.abs(d) < Kernel.STANDARD_PRECISION) {
complex[roots] = 0;
real[roots++] = -b / (2.0 * a);
} else {
if (d < 0.0) {
// If d < 0.0, then there are 2 complex roots
d = Math.sqrt(-d);
complex[0] = d / (2.0 * a);
real[0] = -b / (2.0 * a);
complex[1] = -complex[0];
real[1] = real[0];
roots = 2;
} else {
d = Math.sqrt(d);
// For accuracy, calculate one root using:
// (-b +/- d) / 2a
// and the other using:
// 2c / (-b +/- d)
// Choose the sign of the +/- so that b+d gets larger in
// magnitude
if (b < 0.0) {
d = -d;
}
double q = (b + d) / -2.0;
// We already tested a for being 0 above
complex[roots] = 0;
real[roots++] = q / a;
complex[roots] = 0;
real[roots++] = c / q;
}
}
}
return roots;
}
/**
* Solves the cubic whose coefficients are in the <code>eqn</code> array and
* places the non-complex roots back into the same array, returning the
* number of roots. The solved cubic is represented by the equation:
*
* <pre>
* eqn = {c, b, a, d}
* dx^3 + ax^2 + bx + c = 0
* </pre>
*
* A return value of -1 is used to distinguish a constant equation that
* might be always 0 or never 0 from an equation that has no zeroes.
*
* @param eqn
* an array containing coefficients for a cubic
* @return the number of roots, or -1 if the equation is a constant.
*/
final public int solveCubic(double eqn[]) {
return solveCubic(eqn, eqn, Kernel.STANDARD_PRECISION);
}
/*
* adapted from gsl/poly/solve_cubic.c
*
* added solveQuadratic() added fixRoots() removed sorting of roots
*
* solve_cubic.c - finds the real roots of x^3 + a x^2 + b x + c = 0
*/
@Override
final public int solveCubic(double eqn[], double res[], double eps) {
return solveCubicS(eqn, res, eps);
}
final static public int solveCubicS(double eqn[], double res[],
double eps) {
int roots = 0;
double d = eqn[3];
if (Math.abs(d) < eps) {
// The cubic has degenerated to quadratic (or line or ...).
return solveQuadraticS(eqn, res, eps);
}
double a = eqn[2] / d;
double b = eqn[1] / d;
double c = eqn[0] / d;
double q = (a * a - 3 * b); // D(q) = 2aD(a) + 3D(b)
double r = (2 * a * a * a - 9 * a * b + 27 * c); // D(r) = (3aa-9b)D(a)
// - 9aD(b) + 27D(c)
double Q = q / 9; // D(Q) = D(q)/9
double R = r / 54;
double Q3 = Q * Q * Q;
double R2 = R * R;
double CR2 = 729 * r * r; // D(CR2) = 729*2rD(r) ( |1458r(3aa-9b) -
// 8748q*2a| + |-9a*1458r -8748q*3| +
// |27*1458r| )*kernel.STANDARD_PRECISION
double CQ3 = 2916 * q * q * q; // D(CQ3) = 2916*3qD(q) (D(root) ~
// D(2sqrt(Q))= -1/sqrt(Q) D(Q),
// |D(root)| |2a+3|/sqrt(9q)
// *kernel.STANDARD_PRECISION
// Application.debug("Q = "+Q+" R = "+R+" Q3 = "+Q3+" R2 = "+R2+" CR2 =
// "+CR2+" CQ3 = "+CQ3);
// Application.debug(Math.abs(CR2 - CQ3)+"");
// changed back to original GGB-1725
// if (Math.abs(R) < Kernel.STANDARD_PRECISION
// && Math.abs(Q) < Kernel.STANDARD_PRECISION)
if (R == 0 && Q == 0)
{
res[roots++] = -a / 3;
res[roots++] = -a / 3;
res[roots++] = -a / 3;
return 3;
}
// Michael Borcherds changed to check CR2 equal to CQ3 to first 8
// significant digits
// important for eg y=(ax-b)^2(cx-d)
// |D(CR2-CQ3)|< (|r(aa-3b) - 4qa| + |-3ar -6q| + |9r|)*13122*sqrt(q) /
// |2a+3| *kernel.STANDARD_PRECISION
// for simplicity, it (may be) about 10* max(CR2,CR3)/|2a+3| *
// kernel.STANDARD_PRECISION
// else if (Math.abs(CR2 - CQ3) < Math.max(CR2, CQ3) *
// kernel.STANDARD_PRECISION)
// // else if (CR2 == CQ3)
else if (Math.abs(CR2 - CQ3) < Math.max(CR2, CQ3) * 10
/ Math.max(1, Math.abs(2 * a + 3)) * Kernel.STANDARD_PRECISION) // else
// 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 +/-
// kernel.STANDARD_PRECISION
// i
// close to the real axis.
double sqrtQ = Math.sqrt(Q);
if (R > 0) {
res[roots++] = -2 * sqrtQ - a / 3;
res[roots++] = sqrtQ - a / 3;
res[roots++] = sqrtQ - a / 3;
} else {
res[roots++] = -sqrtQ - a / 3;
res[roots++] = -sqrtQ - a / 3;
res[roots++] = 2 * sqrtQ - a / 3;
}
return 3;
} else 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[roots++] = norm * Math.cos(theta / 3) - a / 3;
res[roots++] = norm * Math.cos((theta + 2.0 * Math.PI) / 3) - a / 3;
res[roots++] = norm * Math.cos((theta - 2.0 * Math.PI) / 3) - a / 3;
// GeoGebra addition
// TODO: find a better way to deal with this
if (res != eqn) {
// possible to fix roots
fixRoots(res, eqn);
}
return 3;
} else {
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[roots++] = A + B - a / 3;
return 1;
}
}// */
/*
* This pruning step is necessary since solveCubic uses the cosine function
* to calculate the roots when there are 3 of them. Since the cosine method
* can have an error of +/- 1E-14 we need to make sure that we don't make
* any bad decisions due to an error.
*
* If the root is not near one of the endpoints, then we will only have a
* slight inaccuracy in calculating the x intercept which will only cause a
* slightly wrong answer for some points very close to the curve. While the
* results in that case are not as accurate as they could be, they are not
* disastrously inaccurate either.
*
* On the other hand, if the error happens near one end of the curve, then
* our processing to reject values outside of the t=[0,1] range will fail
* and the results of that failure will be disastrous since for an entire
* horizontal range of test points, we will either overcount or undercount
* the crossings and get a wrong answer for all of them, even when they are
* clearly and obviously inside or outside the curve.
*
* To work around this problem, we try a couple of Newton-Raphson iterations
* to see if the true root is closer to the endpoint or further away. If it
* is further away, then we can stop since we know we are on the right side
* of the endpoint. If we change direction, then either we are now being
* dragged away from the endpoint in which case the first condition will
* cause us to stop, or we have passed the endpoint and are headed back. In
* the second case, we simply evaluate the slope at the endpoint itself and
* place ourselves on the appropriate side of it or on it depending on that
* result.
*/
private static void fixRoots(double res[], double eqn[]) {
double myepsilon = 1E-5;
for (int i = 0; i < 3; i++) {
double t = res[i];
if (Math.abs(t) < myepsilon) {
res[i] = findZero(t, 0, eqn);
} else if (Math.abs(t - 1) < myepsilon) {
res[i] = findZero(t, 1, eqn);
}
}
}
private static double solveEqn(double eqn[], int order, double t) {
double v = eqn[order];
int counter = order;
while (--counter >= 0) {
v = v * t + eqn[counter];
}
return v;
}
private static double findZero(double init, double target, double eqn[]) {
double t = init;
double slopeqn[] = { eqn[1], 2 * eqn[2], 3 * eqn[3] };
double slope;
double origdelta = 0;
double origt = t;
while (true) {
slope = solveEqn(slopeqn, 2, t);
if (slope == 0.0) {
// At a local minima - must return
return t;
}
double y = solveEqn(eqn, 3, t);
if (y == 0.0) {
// Found it! - return it
return t;
}
// assert(slope != 0 && y != 0);
double delta = -(y / slope);
// assert(delta != 0);
if (origdelta == 0.0) {
origdelta = delta;
}
if (t < target) {
if (delta < 0) {
return t;
}
} else if (t > target) {
if (delta > 0) {
return t;
}
} else {
return (delta > 0 ? (target + java.lang.Double.MIN_VALUE)
: (target - java.lang.Double.MIN_VALUE));
}
double newt = t + delta;
if (MyDouble.exactEqual(t, newt)) {
// The deltas are so small that we aren't moving...
return t;
}
if (delta * origdelta < 0) {
// We have reversed our path.
int tag = (origt < t ? getTag(target, origt, t)
: getTag(target, t, origt));
if (tag != INSIDE) {
// Local minima found away from target - return the middle
return (origt + t) / 2;
}
// Local minima somewhere near target - move to target
// and let the slope determine the resulting t.
t = target;
} else {
t = newt;
}
}
}
private static final int BELOW = -2;
private static final int LOWEDGE = -1;
private static final int INSIDE = 0;
private static final int HIGHEDGE = 1;
private static final int ABOVE = 2;
/*
* Determine where coord lies with respect to the range from low to high. It
* is assumed that low <= high. The return value is one of the 5 values
* BELOW, LOWEDGE, INSIDE, HIGHEDGE, or ABOVE.
*/
private static int getTag(double coord, double low, double high) {
if (coord <= low) {
return (coord < low ? BELOW : LOWEDGE);
}
if (coord >= high) {
return (coord > high ? ABOVE : HIGHEDGE);
}
return INSIDE;
}
/* ************************************************* */
/**
* Calculates all roots of a polynomial given by eqn using Laguerres method.
* Polishes roots found. The roots are stored in eqn again.
*
* @param eqn
* coefficients of polynomial
*/
private int laguerreAll(double[] eqn) {
// for fast evaluation of polynomial (used for root polishing)
PolyFunction polyFunc = new PolyFunction(eqn);
PolyFunction derivFunc = polyFunc.getDerivative();
Complex[] complexRoots = null;
try {
if (laguerreSolver == null) {
laguerreSolver = new LaguerreSolver();
}
complexRoots = laguerreSolver.solveAllComplex(eqn, LAGUERRE_START,
LAGUERRE_MAX_EVAL);
} catch (ArithmeticException e) {
Log.warn("Too many iterations. Degree: " + eqn.length);
} catch (RuntimeException e) {
Log.error("EquationSolver.LaguerreSolver: "
+ e.getLocalizedMessage());
}
if (complexRoots == null) {
complexRoots = new Complex[0];
}
// sort complexRoots by real part into laguerreRoots
double[] laguerreRoots = new double[complexRoots.length];
for (int i = 0; i < laguerreRoots.length; i++) {
laguerreRoots[i] = complexRoots[i].getReal();
}
Arrays.sort(laguerreRoots);
// get the roots from Laguerre method
int realRoots = 0;
double root;
for (int i = 0; i < laguerreRoots.length; i++) {
// System.out.println("laguerreRoots[" + i + "] = " +
// laguerreRoots[i]);
// let's polish all complex roots to get all real roots
root = laguerreRoots[i];
// check if root is bounded in intervall [root-eps, root+eps]
double left = i == 0 ? root - 1 : (root + laguerreRoots[i - 1]) / 2;
double right = i == laguerreRoots.length - 1 ? root + 1
: (root + laguerreRoots[i + 1]) / 2;
double f_left = polyFunc.value(left);
double f_right = polyFunc.value(right);
boolean bounded = f_left * f_right < 0.0;
try {
if (bounded) {
if (rootFinderBrent == null) {
rootFinderBrent = new BrentSolver();
}
// small f'(root): don't go too far from our laguerre root !
double brentRoot = rootFinderBrent.solve(
AlgoRootNewton.MAX_ITERATIONS, polyFunc,
left, right/* , root */);
if (Math.abs(polyFunc.value(brentRoot)) < Math
.abs(polyFunc.value(root))) {
root = brentRoot;
}
// System.out.println("Polish bisectNewtonRaphson: " +
// root);
} else {
if (rootFinderNewton == null) {
rootFinderNewton = new NewtonSolver();
}
// the root is not bounded: give Mr. Newton a chance
double newtonRoot = rootFinderNewton.solve(
AlgoRootNewton.MAX_ITERATIONS, polyFunc,
left, right,
root);
if (Math.abs(polyFunc.value(newtonRoot)) < Math
.abs(polyFunc.value(root))) {
root = newtonRoot;
}
// System.out.println("Polish newtonRaphson: " + root);
}
} catch (RuntimeException e) {
// Application.debug("Polish FAILED: ");
// polishing failed: maybe we have an extremum here
// try to find a local extremum
try {
if (rootFinderBrent == null) {
rootFinderBrent = new BrentSolver();
}
if (left < right) {
double brentRoot = rootFinderBrent.solve(
AlgoRootNewton.MAX_ITERATIONS, derivFunc,
left, right);
if (Math.abs(polyFunc.value(brentRoot)) < Math
.abs(polyFunc.value(root))) {
root = brentRoot;
}
}
// System.out.println(" find extremum successfull: " +
// root);
} catch (RuntimeException ex) {
//Log.debug(ex.getMessage());
}
}
// check if the found root is really ok
double[] val = polyFunc.evaluateDerivFunc(root); // get f(root) and
// f'(root)
double error = Math.abs(val[0]); // | f(root) |
double slope = Math.abs(val[1]);
boolean success;
if (slope < 1) {
success = error < LAGUERRE_EPS;
} else {
success = error < LAGUERRE_EPS * slope;
}
if (success) {
// Application.debug("FOUND ROOT: " + root);
eqn[realRoots] = root;
realRoots++;
} else {
// Application.debug("no root: " + root + ", error " + error);
}
}
return realRoots;
}
/**
* Calculates all roots of a polynomial given by eqn using Laguerres method.
* Polishes roots found. The roots are stored in eqn again.
*/
private int laguerreAllComplex(double[] real, double[] complex) {
Complex[] complexRoots = null;
try {
if (laguerreSolver == null) {
laguerreSolver = new LaguerreSolver();
}
complexRoots = laguerreSolver.solveAllComplex(real, LAGUERRE_START);
} catch (Exception e) {
Log.debug("Problem solving with LaguerreSolver"
+ e.getLocalizedMessage());
return 0;
}
// sort by real part & remove duplicates
TreeSet<Complex> sortedSet = new TreeSet<Complex>(getComparatorReal());
for (int i = 0; i < complexRoots.length; i++) {
sortedSet.add(complexRoots[i]);
}
int roots = 0;
Complex temp;
Iterator<Complex> iterator = sortedSet.iterator();
while (iterator.hasNext()) {
temp = iterator.next();
real[roots] = temp.getReal();
complex[roots] = temp.getImaginary();
roots++;
}
return roots;
}
// avoid too big polynomial coefficients in laguerreAll
// private static final int MAX_POLYNOMIAL_COEFFICIENT = 1000;
//
// // avoid huge coefficients
// private void checkCoefficients(double [] eqn) {
// double max = Math.abs(eqn[0]);
// double temp;
// for (int i=1; i < eqn.length; i++) {
// temp = Math.abs(eqn[i]);
// if (temp > max) max = temp;
// }
// if (max > MAX_POLYNOMIAL_COEFFICIENT) {
// Application.debug("changed coefficients");
// double factor = MAX_POLYNOMIAL_COEFFICIENT/max;
// for (int i=0; i < eqn.length; i++) {
// eqn[i] *= factor;
// }
// }
// }
/*
* solve_quartic.c - finds the real roots of x^4 + a x^3 + b x^2 + c x + d =
* 0
*
* modified from GSL Quartic extension
* http://www.network-theory.co.uk/download/gslextras/Quartic/
*/
@Override
public int solveQuartic(double eqn[], double res[], double eps) {
if (Math.abs(eqn[4]) < 0) {
return solveCubic(eqn, res, Kernel.STANDARD_PRECISION);
}
double a = eqn[3] / eqn[4], b = eqn[2] / eqn[4], c = eqn[1] / eqn[4],
d = eqn[0] / eqn[4];
/*
* This code is based on a simplification of the algorithm from
* zsolve_quartic.c for real roots
*/
double u[] = new double[3], v[] = new double[3], zarr[] = new double[4];
double aa, pp, qq, rr, rc, sc, tc, mt;
double w1r, w1i, w2r, w2i, w3r;
double v1, v2, arg, theta;
double disc, h;
int k1 = 0, k2 = 0;
int roots = 0;
/*
* Deal easily with the cases where the quartic is degenerate. The
* ordering of solutions is done explicitly.
*/
if (0 == b && 0 == c) {
if (0 == d) {
if (a > 0) {
res[roots++] = -a;
res[roots++] = 0.0;
res[roots++] = 0.0;
res[roots++] = 0.0;
} else {
res[roots++] = 0.0;
res[roots++] = 0.0;
res[roots++] = 0.0;
res[roots++] = -a;
}
return 4;
} else if (0 == a) {
if (d > 0) {
return 0;
}
res[roots++] = Math.sqrt(Math.sqrt(-d));
res[roots] = -res[roots - 1];
roots++;
return 2;
}
}
if (0.0 == c && 0.0 == d) {
res[roots++] = 0.0;
res[roots++] = 0.0;
double[] res2 = new double[3];
res2[2] = 1.0;
res2[1] = a;
res2[0] = b;
int n = solveQuadratic(res2, res2, eps);
res[roots++] = res2[0];
res[roots++] = res2[1];
// if (gsl_poly_solve_quadratic(1.0,a,b,x2,x3)==0) {
if (n == 0) {
mt = 3;
} else {
mt = 1;
}
} else {
/*
* For non-degenerate solutions, proceed by constructing and solving
* the resolvent cubic
*/
aa = a * a;
pp = b - (3.0 / 8.0) * aa;
qq = c - (1.0 / 2.0) * a * (b - (1.0 / 4.0) * aa);
rr = d - (1.0 / 4.0)
* (a * c - (1.0 / 4.0) * aa * (b - (3.0 / 16.0) * aa));
rc = (1.0 / 2.0) * pp;
sc = (1.0 / 4.0) * ((1.0 / 4.0) * pp * pp - rr);
tc = -((1.0 / 8.0) * qq * (1.0 / 8.0) * qq);
/*
* This code solves the resolvent cubic in a convenient fashion for
* this implementation of the quartic. If there are three real
* roots, then they are placed directly into u[]. If two are
* complex, then the real root is put into u[0] and the real and
* imaginary part of the complex roots are placed into u[1] and
* u[2], respectively. Additionally, this calculates the
* discriminant of the cubic and puts it into the variable disc.
*/
{
double qcub = (rc * rc - 3 * sc);
double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc);
double Q = qcub / 9;
double R = rcub / 54;
double Q3 = Q * Q * Q;
double R2 = R * R;
double CR2 = 729 * rcub * rcub;
double CQ3 = 2916 * qcub * qcub * qcub;
disc = (CR2 - CQ3) / 2125764.0;
if (0 == R && 0 == Q) {
u[0] = -rc / 3;
u[1] = -rc / 3;
u[2] = -rc / 3;
} else if (MyDouble.exactEqual(CR2, CQ3)) {
double sqrtQ = Math.sqrt(Q);
if (R > 0) {
u[0] = -2 * sqrtQ - rc / 3;
u[1] = sqrtQ - rc / 3;
u[2] = sqrtQ - rc / 3;
} else {
u[0] = -sqrtQ - rc / 3;
u[1] = -sqrtQ - rc / 3;
u[2] = 2 * sqrtQ - rc / 3;
}
} else if (CR2 < CQ3) {
double sqrtQ = Math.sqrt(Q);
double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
theta = Math.acos(R / sqrtQ3);
if (R / sqrtQ3 >= 1.0) {
theta = 0.0;
}
{
double norm = -2 * sqrtQ;
u[0] = norm * Math.cos(theta / 3) - rc / 3;
u[1] = norm * Math.cos((theta + 2.0 * Math.PI) / 3)
- rc / 3;
u[2] = norm * Math.cos((theta - 2.0 * Math.PI) / 3)
- rc / 3;
}
} else {
double sgnR = (R >= 0 ? 1 : -1);
double modR = Math.abs(R);
double sqrt_disc = Math.sqrt(R2 - Q3);
double A = -sgnR * Math.pow(modR + sqrt_disc, 1.0 / 3.0);
double B = Q / A;
double mod_diffAB = Math.abs(A - B);
u[0] = A + B - rc / 3;
u[1] = -0.5 * (A + B) - rc / 3;
u[2] = -(Math.sqrt(3.0) / 2.0) * mod_diffAB;
}
}
/* End of solution to resolvent cubic */
/*
* Combine the square roots of the roots of the cubic resolvent
* appropriately. Also, calculate 'mt' which designates the nature
* of the roots: mt=1 : 4 real roots (disc == 0) mt=2 : 0 real roots
* (disc < 0) mt=3 : 2 real roots (disc > 0)
*/
if (0.0 == disc) {
u[2] = u[1];
}
if (0 >= disc) {
mt = 2;
/*
* One would think that we could return 0 here and exit, since
* mt=2. However, this assignment is temporary and changes to
* mt=1 under certain conditions below.
*/
v[0] = Math.abs(u[0]);
v[1] = Math.abs(u[1]);
v[2] = Math.abs(u[2]);
v1 = Math.max(Math.max(v[0], v[1]), v[2]);
/* Work out which two roots have the largest moduli */
k1 = 0;
k2 = 0;
if (v1 == v[0]) {
k1 = 0;
v2 = Math.max(v[1], v[2]);
} else if (v1 == v[1]) {
k1 = 1;
v2 = Math.max(v[0], v[2]);
} else {
k1 = 2;
v2 = Math.max(v[0], v[1]);
}
if (v2 == v[0]) {
k2 = 0;
} else if (v2 == v[1]) {
k2 = 1;
} else {
k2 = 2;
}
if (0.0 <= u[k1]) {
w1r = Math.sqrt(u[k1]);
w1i = 0.0;
} else {
w1r = 0.0;
w1i = Math.sqrt(-u[k1]);
}
if (0.0 <= u[k2]) {
w2r = Math.sqrt(u[k2]);
w2i = 0.0;
} else {
w2r = 0.0;
w2i = Math.sqrt(-u[k2]);
}
} else {
mt = 3;
if (0.0 == u[1] && 0.0 == u[2]) {
arg = 0.0;
} else {
arg = Math.sqrt(Math.sqrt(u[1] * u[1] + u[2] * u[2]));
}
theta = Math.atan2(u[2], u[1]);
w1r = arg * Math.cos(theta / 2.0);
w1i = arg * Math.sin(theta / 2.0);
w2r = w1r;
w2i = -w1i;
}
/* Solve the quadratic to obtain the roots to the quartic */
w3r = qq / 8.0 * (w1i * w2i - w1r * w2r) / (w1i * w1i + w1r * w1r)
/ (w2i * w2i + w2r * w2r);
h = a / 4.0;
zarr[0] = w1r + w2r + w3r - h;
zarr[1] = -w1r - w2r + w3r - h;
zarr[2] = -w1r + w2r - w3r - h;
zarr[3] = w1r - w2r - w3r - h;
/* Arrange the roots into the variables z0, z1, z2, z3 */
if (2 == mt) {
if (u[k1] >= 0 && u[k2] >= 0) {
mt = 1;
res[roots++] = zarr[0];
res[roots++] = zarr[1];
res[roots++] = zarr[2];
res[roots++] = zarr[3];
} else {
return 0;
}
} else {
res[roots++] = zarr[0];
res[roots++] = zarr[1];
}
}
/* Sort the roots as usual */
if (1 == mt) {
/*
* Roots are all real, sort them by the real part if (*x0 > *x1)
* SWAPD (*x0, *x1); if (*x0 > *x2) SWAPD (*x0, *x2); if (*x0 > *x3)
* SWAPD (*x0, *x3);
*
* if (*x1 > *x2) SWAPD (*x1, *x2); if (*x2 > *x3) { SWAPD (*x2,
* *x3); if (*x1 > *x2) SWAPD (*x1, *x2); }
*/
return 4;
}
/*
* 2 real roots if (*x0 > *x1) SWAPD (*x0, *x1);
*/
return 2;
}
/**
* Returns a comparator for Complex objects. (sorts on real part coordinate)
* If equal does return zero (deletes duplicates!)
*
* @return comparator for complex numbers
*/
public static Comparator<Complex> getComparatorReal() {
if (comparatorReal == null) {
comparatorReal = new Comparator<Complex>() {
@Override
public int compare(Complex itemA, Complex itemB) {
double compReal = itemA.getReal() - itemB.getReal();
if (Kernel.isZero(compReal)) {
double compImaginary = itemA.getImaginary()
- itemB.getImaginary();
// if real parts equal, sort on imaginary
if (!Kernel.isZero(compImaginary)) {
return compImaginary < 0 ? -1 : +1;
}
// return 0 -> remove duplicates!
return 0;
}
return compReal < 0 ? -1 : +1;
}
};
}
return comparatorReal;
}
private static Comparator<Complex> comparatorReal;
}