CaRMtl/rene/zirkel/expression/Quartic.java
2018-09-04 22:51:42 -04:00

278 lines
7.8 KiB
Java

/*
Copyright 2006 Rene Grothmann, modified by Eric Hakenholz
This file is part of C.a.R. software.
C.a.R. is a 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, version 3 of the License.
C.a.R. is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package rene.zirkel.expression;
public class Quartic {
static public int solve(final double c[], final double sol[]) {
double sum = 0;
for (final double element : c)
sum += Math.abs(element);
if (sum < 1e-15) {
sol[0] = 0;
return 1;
}
for (int k = 0; k < c.length; k++)
c[k] /= sum;
int n = 4;
while (Math.abs(c[n]) < 1e-15)
n--;
if (n == 4) {
final double soli[] = new double[4];
return quartic(c, sol, soli);
} else if (n == 3) {
return cubic(c, sol);
} else if (n == 2) {
final double a = -c[1] / (2 * c[2]);
final double r = a * a - c[0] / c[2];
if (r < -1e-10)
return 0;
else if (Math.abs(r) < 1e-10) {
sol[0] = a;
return 1;
} else {
sol[0] = a + Math.sqrt(r);
sol[1] = a - Math.sqrt(r);
return 2;
}
} else if (n == 1) {
sol[0] = -c[0] / c[1];
return 1;
} else {
if (Math.abs(c[0]) < 1e-10) {
sol[0] = 0;
return 1;
}
}
return 0;
}
/*-------------------- Global Function Description Block ----------------------
*
* ***QUARTIC************************************************25.03.98
* Solution of a quartic equation
* ref.: J. E. Hacke, Amer. Math. Monthly, Vol. 48, 327-328, (1941)
* NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
* ******************************************************************
* dd(0:4) (i) vector containing the polynomial coefficients
* sol(1:4) (o) results, real part
* soli(1:4) (o) results, imaginary part
* Nsol (o) number of real solutions
* ==================================================================
* 17-Oct-2004 / Raoul Rausch
* Conversion from Fortran to C
*
*
*-----------------------------------------------------------------------------
*/
static int quartic(final double dd[], final double sol[],
final double soli[]) {
final double AA[] = new double[4], z[] = new double[3];
double a, b, c, d, f, p, q, r, zsol, xK2, xL, xK, sqp, sqm;
int ncube, i;
int Nsol = 0;
if (dd[4] == 0.0) {
return 0;
}
a = dd[4];
b = dd[3];
c = dd[2];
d = dd[1];
f = dd[0];
p = (-3.0 * Math.pow(b, 2) + 8.0 * a * c) / (8.0 * Math.pow(a, 2));
q = (Math.pow(b, 3) - 4.0 * a * b * c + 8.0 * d * Math.pow(a, 2))
/ (8.0 * Math.pow(a, 3));
r = (-3.0 * Math.pow(b, 4) + 16.0 * a * Math.pow(b, 2) * c - 64.0
* Math.pow(a, 2) * b * d + 256.0 * Math.pow(a, 3) * f)
/ (256.0 * Math.pow(a, 4));
// Solve cubic resolvent
AA[3] = 8.0;
AA[2] = -4.0 * p;
AA[1] = -8.0 * r;
AA[0] = 4.0 * p * r - Math.pow(q, 2);
ncube = cubic(AA, z);
zsol = -1.e99;
for (i = 0; i < ncube; i++)
zsol = Math.max(zsol, z[i]); // Not sure C has max fct
z[0] = zsol;
xK2 = 2.0 * z[0] - p;
xK = Math.sqrt(xK2);
xL = q / (2.0 * xK);
sqp = xK2 - 4.0 * (z[0] + xL);
sqm = xK2 - 4.0 * (z[0] - xL);
for (i = 0; i < 4; i++)
soli[i] = 0.0;
if ((sqp >= 0.0) && (sqm >= 0.0)) {
sol[0] = 0.5 * (xK + Math.sqrt(sqp));
sol[1] = 0.5 * (xK - Math.sqrt(sqp));
sol[2] = 0.5 * (-xK + Math.sqrt(sqm));
sol[3] = 0.5 * (-xK - Math.sqrt(sqm));
Nsol = 4;
} else if ((sqp >= 0.0) && (sqm < 0.0)) {
sol[0] = 0.5 * (xK + Math.sqrt(sqp));
sol[1] = 0.5 * (xK - Math.sqrt(sqp));
sol[2] = -0.5 * xK;
sol[3] = -0.5 * xK;
soli[2] = Math.sqrt(-.25 * sqm);
soli[3] = -Math.sqrt(-.25 * sqm);
Nsol = 2;
} else if ((sqp < 0.0) && (sqm >= 0.0)) {
sol[0] = 0.5 * (-xK + Math.sqrt(sqm));
sol[1] = 0.5 * (-xK - Math.sqrt(sqm));
sol[2] = 0.5 * xK;
sol[3] = 0.5 * xK;
soli[2] = Math.sqrt(-0.25 * sqp);
soli[3] = -Math.sqrt(-0.25 * sqp);
Nsol = 2;
} else if ((sqp < 0.0) && (sqm < 0.0)) {
sol[0] = -0.5 * xK;
sol[1] = -0.5 * xK;
soli[0] = Math.sqrt(-0.25 * sqm);
soli[1] = -Math.sqrt(-0.25 * sqm);
sol[2] = 0.5 * xK;
sol[3] = 0.5 * xK;
soli[2] = Math.sqrt(-0.25 * sqp);
soli[3] = -Math.sqrt(-0.25 * sqp);
Nsol = 0;
}
for (i = 0; i < 4; i++)
sol[i] -= b / (4.0 * a);
return Nsol;
}
/*-------------------- Global Function Description Block ----------------------
*
* ***CUBIC************************************************08.11.1986
* Solution of a cubic equation
* Equations of lesser degree are solved by the appropriate formulas.
* The solutions are arranged in ascending order.
* NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
* ******************************************************************
* A(0:3) (i) vector containing the polynomial coefficients
* X(1:L) (o) results
* L (o) number of valid solutions (beginning with X(1))
* ==================================================================
* 17-Oct-2004 / Raoul Rausch
* Conversion from Fortran to C
*
*
*-----------------------------------------------------------------------------
*/
static public int cubic(final double A[], final double X[]) {
final double PI = 3.1415926535897932;
final double THIRD = 1. / 3.;
final double U[] = new double[3];
double W, P, Q, DIS, PHI;
int i;
int L;
if (A[3] != 0.0) {
// cubic problem
W = A[2] / A[3] * THIRD;
P = Math.pow((A[1] / A[3] * THIRD - Math.pow(W, 2)), 3);
Q = -.5 * (2.0 * Math.pow(W, 3) - (A[1] * W - A[0]) / A[3]);
DIS = Math.pow(Q, 2) + P;
if (DIS < 0.0) {
// three real solutions!
// Confine the argument of ACOS to the interval [-1;1]!
PHI = Math.acos(Math
.min(1.0, Math.max(-1.0, Q / Math.sqrt(-P))));
P = 2.0 * Math.pow((-P), (5.e-1 * THIRD));
for (i = 0; i < 3; i++)
U[i] = P * Math.cos((PHI + 2 * ((double) i) * PI) * THIRD)
- W;
X[0] = Math.min(U[0], Math.min(U[1], U[2]));
X[1] = Math.max(Math.min(U[0], U[1]), Math.max(Math.min(U[0],
U[2]), Math.min(U[1], U[2])));
X[2] = Math.max(U[0], Math.max(U[1], U[2]));
L = 3;
} else {
// only one real solution!
DIS = Math.sqrt(DIS);
X[0] = CBRT(Q + DIS) + CBRT(Q - DIS) - W;
L = 1;
}
} else if (A[2] != 0.0) {
// quadratic problem
P = 0.5 * A[1] / A[2];
DIS = Math.pow(P, 2) - A[0] / A[2];
if (DIS > 0.0) {
// 2 real solutions
X[0] = -P - Math.sqrt(DIS);
X[1] = -P + Math.sqrt(DIS);
L = 2;
} else {
// no real solution
L = 0;
}
} else if (A[1] != 0.0) {
// linear equation
X[0] = A[0] / A[1];
L = 1;
} else {
// no equation
L = 0;
}
/*
* ==== perform one step of a newton iteration in order to minimize
* round-off errors ====
*/
for (i = 0; i < L; i++) {
X[i] = X[i] - (A[0] + X[i] * (A[1] + X[i] * (A[2] + X[i] * A[3])))
/ (A[1] + X[i] * (2.0 * A[2] + X[i] * 3.0 * A[3]));
}
return L;
}
static double CBRT(final double Z) {
final double THIRD = 1. / 3.;
if (Z > 0)
return Math.pow(Z, THIRD);
else if (Z < 0)
return -Math.pow(Math.abs(Z), THIRD);
else
return 0;
}
/**
* @param args
*/
public static void main(final String[] args) {
final double a[] = { 0, 3, -3, -6, 5 }, x[] = new double[5], y[] = new double[5];
final int n = quartic(a, x, y);
System.out.println(n + " solutions!");
for (int i = 0; i < n; i++)
System.out.println(x[i]);
}
}