CaRMtl/rene/zirkel/expression/Quartic.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]);
	}

}