private PolynomialRoots()
{
// Do not instantiate
}
/**
* Emulate the F77 sign function.
* @param a double; the value to optionally sign invert
* @param b double; the sign of which determines what to do
* @return double; if b >= 0 then a; else -a
*/
private static double sign(final double a, final double b)
{
return b >= 0 ? a : -a;
}
/**
* LINEAR POLYNOMIAL ROOT SOLVER.
* <p>
* Calculates the root of the linear polynomial:<br>
* q1 * x + q0<br>
* Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
* @param q1 double; coefficient of the x term
* @param q0 double; independent coefficient
* @return Complex[]; the roots of the equation
*/
public static Complex[] linearRoots(final double q1, final double q0)
{
if (q1 == 0)
{
return new Complex[] {}; // No roots; return empty array
}
return linearRoots(q0 / q1);
}
/**
* LINEAR POLYNOMIAL ROOT SOLVER.
* <p>
* Calculates the root of the linear polynomial:<br>
* x + q0<br>
* Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
* @param q0 double; independent coefficient
* @return Complex[]; the roots of the equation
*/
public static Complex[] linearRoots(final double q0)
{
return new Complex[] { new Complex(-q0, 0) };
}
/**
* QUADRATIC POLYNOMIAL ROOT SOLVER
* <p>
* Calculates all real + complex roots of the quadratic polynomial:<br>
* q2 * x^2 + q1 * x + q0<br>
* The code checks internally if rescaling of the coefficients is needed to avoid overflow.
* <p>
* The order of the roots is as follows:<br>
* 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
* negative last).<br>
* 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
* q1 : coefficient of x term q0 : independent coefficient
* @param q2 double; coefficient of the quadratic term
* @param q1 double; coefficient of the x term
* @param q0 double; independent coefficient
* @return Complex[]; the roots of the equation
*/
public static Complex[] quadraticRoots(final double q2, final double q1, final double q0)
{
if (q2 == 0)
{
return linearRoots(q1, q0);
}
return quadraticRoots(q1 / q2, q0 / q2);
}
/**
* QUADRATIC POLYNOMIAL ROOT SOLVER
* <p>
* Calculates all real + complex roots of the quadratic polynomial:<br>
* x^2 + q1 * x + q0<br>
* The code checks internally if rescaling of the coefficients is needed to avoid overflow.
* <p>
* The order of the roots is as follows:<br>
* 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
* negative last).<br>
* 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
* q1 : coefficient of x term q0 : independent coefficient
* @param q1 double; coefficient of the x term
* @param q0 double; independent coefficient
* @return Complex[]; the roots of the equation
*/
public static Complex[] quadraticRoots(final double q1, final double q0)
{
boolean rescale;
double a0, a1;
double k = 0, x, y, z;
// Handle special cases.
if (q0 == 0.0 && q1 == 0.0)
{
// Two real roots at 0,0
return new Complex[] { Complex.ZERO, Complex.ZERO };
}
else if (q0 == 0.0)
{
// Two real roots; one of these is 0,0
// x^2 + q1 * x == x * (x + q1)
Complex nonZeroRoot = new Complex(-q1);
return new Complex[] { q1 > 0 ? Complex.ZERO : nonZeroRoot, q1 <= 0 ? nonZeroRoot : Complex.ZERO };
}
else if (q1 == 0.0)
{
x = Math.sqrt(Math.abs(q0));
if (q0 < 0.0)
{
// Two real roots, symmetrically around 0
return new Complex[] { new Complex(x, 0), new Complex(-x, 0) };
}
else
{
// Two complex roots, symmetrically around 0
return new Complex[] { new Complex(0, x), new Complex(0, -x) };
}
}
else
{
// The general case. Do rescaling, if either squaring of q1/2 or evaluation of
// (q1/2)^2 - q0 will lead to overflow. This is better than to have the solver
// crashed. Note, that rescaling might lead to loss of accuracy, so we only
// invoke it when absolutely necessary.
final double sqrtLPN = Math.sqrt(Double.MAX_VALUE); // Square root of the Largest Positive Number
rescale = (q1 > sqrtLPN + sqrtLPN); // this detects overflow of (q1/2)^2
if (!rescale)
{
x = q1 * 0.5; // we are sure here that x*x will not overflow
rescale = (q0 < x * x - Double.MAX_VALUE); // this detects overflow of (q1/2)^2 - q0
}
if (rescale)
{
x = Math.abs(q1);
y = Math.sqrt(Math.abs(q0));
if (x > y)
{
k = x;
z = 1.0 / x;
a1 = sign(1.0, q1);
a0 = (q0 * z) * z;
}
else
{
k = y;
a1 = q1 / y;
a0 = sign(1.0, q0);
}
}
else
{
a1 = q1;
a0 = q0;
}
// Determine the roots of the quadratic. Note, that either a1 or a0 might
// have become equal to zero due to underflow. But both cannot be zero.
x = a1 * 0.5;
y = x * x - a0;
if (y >= 0.0)
{
// Two real roots
y = Math.sqrt(y);
if (x > 0.0)
{
y = -x - y;
}
else
{
y = -x + y;
}
if (rescale)
{
y = y * k; // very important to convert to original
z = q0 / y; // root first, otherwise complete loss of
}
else // root due to possible a0 = 0 underflow
{
z = a0 / y;
}
return new Complex[] { new Complex(Math.max(y, z), 0), new Complex(Math.min(y, z), 0) };
}
else
{
// Two complex roots (zero real roots)
y = Math.sqrt(-y);
if (rescale)
{
x *= k;
y *= k;
}
return new Complex[] { new Complex(-x, y), new Complex(-x, -y) };
}
}
}
/**
* CUBIC POLYNOMIAL ROOT SOLVER.
* <p>
* Calculates all (real and complex) roots of the cubic polynomial:<br>
* c3 * x^3 + c2 * x^2 + c1 * x + c0<br>
* The first real root (which always exists) is obtained using an optimized Newton-Raphson scheme. The other remaining roots
* are obtained through composite deflation into a quadratic.
* <P>
* The cubic root solver can handle any size of cubic coefficients and there is no danger of overflow due to proper
* rescaling of the cubic polynomial. The order of the roots is as follows: 1) For real roots, the order is according to
* their algebraic value on the number scale (largest positive first, largest negative last). 2) Since there can be only one
* complex conjugate pair root, no order is necessary. 3) All real roots precede the complex ones.
* @param c3 double; coefficient of the cubic term
* @param c2 double; coefficient of the quadratic term
* @param c1 double; coefficient of the linear term
* @param c0 double; coefficient of the independent term
* @return Complex[]; array of Complex with all the roots
*/
public static Complex[] cubicRoots(final double c3, final double c2, final double c1, final double c0) 