package de.jtdev.jcurve.ui;

public class Solver {

    /**
     * Solves the linear equation 0 = c1*x + c0 for x. Result will be stored in
     * the given arrays. The number of found solutions will be returned.
     * @param c0   First parameter of the equation
     * @param c1   Second parameter of the equation
     * @param real Array of doubles with minimum length of 1, which the real
     *             part of the result is stored in
     * @param imag Array of doubles with minimum length of 1, which the
     *             imaginary part of the solution is stored in (will be set to
     *             0.0 if a solution is found).
     * @return If the equation has no solution 0 is returned. In case of
     *         infinit sollutions -1 is returned to distiguish this two cases.
     *         Otherwise 1 to indicate a solution was found.
     */
    public static int solveLinear(double c0, double c1, double[] real, double[] imag) {
        if (c1 == 0.0d) {
            if (c0 == 0.0d) {
                return -1;
            } else {
                return 0;
            }
        }

        real[0] = -c0 / c1;
        imag[0] = 0.0d;
        return 1;
    }

    /**
     * Solves the quadric equation 0 = c2*x^2 + c1*x + c0 for x. Result will be
     * stored in the given arrays, which might be complex. The number of found
     * solutions will be returned.
     *
     * Note: Same solution might be returned twice.
     *
     * @param c0   First parameter of the equation
     * @param c1   Second parameter of the equation
     * @param c2   Third parameter of the equation
     * @param real Array of doubles with minimum length of 2, which the real
     *             part of the result is stored in
     * @param imag Array of doubles with minimum length of 2, which the
     *             imaginary part of the solution is stored in (will be set to
     *             0.0 if a solution is found).
     * @return If the equation has no solution 0 is returned. In case of
     *         infinit sollutions -1 is returned to distiguish this two cases.
     *         Otherwise the number of found sollutions is returned
     */
    public static int solveQuadric(double c0, double c1, double c2, double[] real, double[] imag) {
        if (c2 == 0.0d)
            return solveLinear(c0, c1, real, imag);

        double q = c0 / c2;
        double p = c1 / c2;
        double d = p*p / 4.0d - q;

        if (d >= 0.0d) {
            real[0] = -p / 2.0d + Math.sqrt(d);
            real[1] = -p / 2.0d - Math.sqrt(d);
            imag[0] = 0.0d;
            imag[1] = 0.0d;
        } else {
            real[0] = -p / 2.0d;
            real[1] = real[0];
            imag[0] = Math.sqrt(-d);
            imag[1] = -imag[0];
        }
        return 2;
    }

    /**
     * Solves the cubic equation 0 = c3*x^3 + c2*x^2 + c1*x + c0 for x. Result
     * will be stored in the given arrays, which might be complex. The number
     * of found solutions will be returned.
     *
     * Note: Same solution might be returned multiple times.
     *
     * @param c0   First parameter of the equation
     * @param c1   Second parameter of the equation
     * @param c2   Third parameter of the equation
     * @param c3   Fourth parameter of the equation
     * @param real Array of doubles with minimum length of 3, which the real
     *             part of the result is stored in
     * @param imag Array of doubles with minimum length of 3, which the
     *             imaginary part of the solution is stored in (will be set to
     *             0.0 if a solution is found).
     * @return If the equation has no solution 0 is returned. In case of
     *         infinit sollutions -1 is returned to distiguish this two cases.
     *         Otherwise the number of found sollutions is returned
     */
    public static int solveCubic(double c0, double c1, double c2, double c3, double[] real, double[] imag) {
        if (c3 == 0.0d)
            return solveQuadric(c0, c1, c2, real, imag);

        double t = c0 / c3;
        double s = c1 / c3;
        double r = c2 / c3 / 3.0d;

        double p = s - 3.0d*r*r;
        double q = 2.0d*r*r*r - r*s + t;
        double d = q*q/4.0d + p*p*p/27.0d;

        if (d < 0.0d) { // only real sollution, case d == 0.0d is handled in else
            double R = 2.0d * Math.sqrt(Math.abs(p/3.0d));
            double phi = Math.acos(-4.0d*q/Math.pow(R,3.0d));
            double y1 =  R * Math.cos(phi/3.0d);
            double y2 = -R * Math.cos((phi-Math.PI)/3.0d);
            double y3 = -R * Math.cos((phi+Math.PI)/3.0d);
            real[0] = y1-r;
            real[1] = y2-r;
            real[2] = y3-r;
            imag[0] = 0;
            imag[1] = 0;
            imag[2] = 0;
        } else {
            double h, u, v;

            h = -0.5d * q + Math.sqrt(d);
            if (h < 0.0d) {
                u = -Math.pow(-h, 1.0d/3.0d);
            } else {
                u =  Math.pow( h, 1.0d/3.0d);
            }

            h = -0.5d * q - Math.sqrt(d);
            if (h < 0.0d) {
                v = -Math.pow(-h, 1.0d/3.0d);
            } else {
                v =  Math.pow( h, 1.0d/3.0d);
            }

            double y1 = u + v;
            double y2 = -y1 / 2.0d;
            double y3 = (v-u) * Math.sqrt(0.75d);
            real[0] = y1-r;
            real[1] = y2-r;
            real[2] = real[1];
            imag[0] = 0;
            imag[1] = y3;
            imag[2] = -imag[1];
        }
        return 3;
    }

    /**
     * Solves the cubic equation 0 = c4*x^4 + c3*x^3 + c2*x^2 + c1*x + c0 for
     * x. Result will be stored in the given arrays, which might be complex.
     * The number of found solutions will be returned.
     *
     * Note: Same solution might be returned multiple times.
     *
     * @param c0   First parameter of the equation
     * @param c1   Second parameter of the equation
     * @param c2   Third parameter of the equation
     * @param c3   Fourth parameter of the equation
     * @param c4   Fifth parameter of the equation
     * @param real Array of doubles with minimum length of 4, which the real
     *             part of the result is stored in
     * @param imag Array of doubles with minimum length of 4, which the
     *             imaginary part of the solution is stored in (will be set to
     *             0.0 if a solution is found).
     * @return If the equation has no solution 0 is returned. In case of
     *         infinit sollutions -1 is returned to distiguish this two cases.
     *         Otherwise the number of found sollutions is returned
     */
    public static int solveQuartic(double c0, double c1, double c2, double c3, double c4, double[] real, double[] imag) {
        if (c4 == 0.0d)
            return solveCubic(c0, c1, c2, c3, real, imag);

        double d = c0 / c4;
        double c = c1 / c4;
        double b = c2 / c4;
        double a = c3 / c4;

        double h1 = 0.25d*a*c - b*b/12.0d - d;
        double h2 = a*b*c/24.0d - 0.125d*c*c - b*b*b/108.0d - 0.125d*a*a*d + b*d/3.0d;

        solveCubic(h2, h1, 0.0d, 1.0d, real, imag);

        double rl = real[0];
        h1 = 0.25d*a*a - 2.0d*b/3.0d + 2.0d*rl + 1e-15;
        h2 = b*b/36.0d + rl*rl + b*rl/3.0d - d + 1e-15;

        if (Math.abs(h1) < 1e-10)
            h1 = 0.0d;
        if (Math.abs(h2) < 1e-10)
            h2 = 0.0d;

        double alpha = Math.sqrt(h1);
        double beta = Math.sqrt(h2);

        double probe = 2.0d*alpha*beta - (a*b/6.0d + a*rl - c);
        if (Math.abs(probe) > 0.000001) {
            alpha = -alpha;
            probe = 2.0d*alpha*beta - (a*b/6.0d + a*rl - c);
        }

        h1 = 0.5d*a + alpha;
        h2 = b/6.0d + rl + beta;
        solveQuadric(h2, h1, 1.0d, real, imag);

        real[2] = real[0];
        imag[2] = imag[0];
        real[3] = real[1];
        imag[3] = imag[1];

        h1 = 0.5d*a - alpha;
        h2 = b/6.0d + rl - beta;
        solveQuadric(h2, h1, 1.0d, real, imag);
        return 4;
    }

}