import java.awt.geom.QuadCurve2D;

public class CubicSolver {
    public static int solveCubicOld(double eqn[], double res[]) {
        // From Numerical Recipes, 5.6, Quadratic and Cubic Equations
        double d = eqn[3];
        if (d == 0.0) {
            // The cubic has degenerated to quadratic (or line or ...).
            return QuadCurve2D.solveQuadratic(eqn, res);
        }
        double a = eqn[2] / d;
        double b = eqn[1] / d;
        double c = eqn[0] / d;
        int roots = 0;
        double Q = (a * a - 3.0 * b) / 9.0;
        double R = (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) / 54.0;
        double R2 = R * R;
        double Q3 = Q * Q * Q;
        a = a / 3.0;
        if (R2 < Q3) {
            double theta = Math.acos(R / Math.sqrt(Q3));
            Q = -2.0 * Math.sqrt(Q);
            if (res == eqn) {
                // Copy the eqn so that we don't clobber it with the
                // roots.  This is needed so that fixRoots can do its
                // work with the original equation.
                eqn = new double[4];
                System.arraycopy(res, 0, eqn, 0, 4);
            }
            res[roots++] = Q * Math.cos(theta / 3.0) - a;
            res[roots++] = Q * Math.cos((theta + Math.PI * 2.0)/ 3.0) - a;
            res[roots++] = Q * Math.cos((theta - Math.PI * 2.0)/ 3.0) - a;
            fixRoots(res, eqn);
        } else {
            boolean neg = (R < 0.0);
            double S = Math.sqrt(R2 - Q3);
            if (neg) {
                R = -R;
            }
            double A = Math.pow(R + S, 1.0 / 3.0);
            if (!neg) {
                A = -A;
            }
            double B = (A == 0.0) ? 0.0 : (Q / A);
            res[roots++] = (A + B) - a;
        }
        return roots;
    }

    public static int solveCubicNew(double eqn[], double res[]) {
        // From Graphics Gems:
        // http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
        final double d = eqn[3];
        if (d == 0) {
            return QuadCurve2D.solveQuadratic(eqn, res);
        }
        /* normal form: x^3 + Ax^2 + Bx + C = 0 */
        final double A = eqn[2] / d;
        final double B = eqn[1] / d;
        final double C = eqn[0] / d;


        //  substitute x = y - A/3 to eliminate quadratic term:
        //     x^3 +px + q = 0

        double sq_A = A * A;
        double p = 1.0/3 * (-1.0/3 * sq_A + B);
        double q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);

        /* use Cardano's formula */

        double cb_p = p * p * p;
        double D = q * q + cb_p;

        int num;
        // XXX: we consider 0 to be anything within 1e-9 of 0.
        // Is this really right? Maybe we should use a bound that changes
        // with the number being tested (math.ulp would do that, but
        // what input do we give it? And do we scale it, or will
        // within(D, 0, math.ulp(somenumber)) be good enough).
        if (iszero(D)) {
            // XXX: do we really need iszero for q? All we do with it is
            // take it's cube root, which works fine even for Double.MIN_VALUE
            // Then again, if we remove it, we will get two extremely close
            // roots for equations where there is only one zero root.
            // We probably should use something like iszero, but much more
            // scrict - so, within(q, 0, 2*Math.ulp(0));
            if (iszero(q)) { /* one triple solution */
                res[ 0 ] = 0;
                num = 1;
            } else { /* one single and one double solution */
                final double u = Math.cbrt(-q);
                res[ 0 ] = 2*u;
                res[ 1 ] = -u;
                num = 2;
            }
        } else if (D < 0) { /* Casus irreducibilis: three real solutions */
            final double phi = 1.0/3 * Math.acos(-q / Math.sqrt(-cb_p));
            final double t = 2 * Math.sqrt(-p);

            if (res == eqn) {
                // Copy the eqn so that we don't clobber it with the
                // roots.  This is needed so that fixRoots can do its
                // work with the original equation.
                eqn = new double[4];
                System.arraycopy(res, 0, eqn, 0, 4);
            }
            res[ 0 ] =  ( t * Math.cos(phi));
            res[ 1 ] =  (-t * Math.cos(phi + Math.PI / 3));
            res[ 2 ] =  (-t * Math.cos(phi - Math.PI / 3));
            num = 3;
        } else { /* one real solution */
            final double sqrt_D = Math.sqrt(D);
            final double u = Math.cbrt(sqrt_D - q);
            final double v = - Math.cbrt(sqrt_D + q);

            res[ 0 ] =  u + v;
            num = 1;
        }

        /* resubstitute */

        final double sub = 1.0/3 * A;

        for (int i = 0; i < num; ++i)
            res[ i ] -= sub;

        if (num == 3) {
            fixRoots(res, eqn);
        }

        return num;
    }

    private static boolean iszero(double x) {
        return within(x, 0, 1e-9);
    }
    
    private static boolean within(final double x, final double y, final double err) {
        final double d = y - x;
        return (d <= err && d >= -err);
    }

    private static void fixRoots(double res[], double eqn[]) {
        final double EPSILON = 1E-5;
        for (int i = 0; i < 3; i++) {
            double t = res[i];
            if (Math.abs(t) < EPSILON) {
                res[i] = findZero(t, 0, eqn);
            } else if (Math.abs(t - 1) < EPSILON) {
                res[i] = findZero(t, 1, eqn);
            }
        }
    }

    private static double solveEqn(double eqn[], int order, double t) {
        double v = eqn[order];
        while (--order >= 0) {
            v = v * t + eqn[order];
        }
        return v;
    }

    private static double findZero(double t, double target, double eqn[]) {
        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) {
                // At a local minima - must return
                return t;
            }
            double y = solveEqn(eqn, 3, t);
            if (y == 0) {
                // Found it! - return it
                return t;
            }
            // assert(slope != 0 && y != 0);
            double delta = - (y / slope);
            // assert(delta != 0);
            if (origdelta == 0) {
                origdelta = delta;
            }
            if (t < target) {
                if (delta < 0) return t;
            } else if (t > target) {
                if (delta > 0) return t;
            } else { /* t == target */
                return (delta > 0
                        ? (target + java.lang.Double.MIN_VALUE)
                        : (target - java.lang.Double.MIN_VALUE));
            }
            double newt = t + delta;
            if (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;

    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;
    }

    public static void trySolve3(double r1, double r2, double r3) {
        // (x-r1)*(x-r2)*(x-r3)
        // = (xx - r1x - r2x + r1r2)*(x-r3)
        // = xxx - r1xx - r2xx + r1r2x - r3xx + r1r3x + r2r3x - r1r2r3
        // = xxx - (r1 + r2 + r3)xx + (r1r2 + r2r3 + r3r1)x - r1r2r3
        System.out.println("solving: (x - "+r1+") * (x - "+r2+") * (x - "+r3+") = 0");
        double eqn[] = new double[4];
        eqn[3] = 1.0;
        eqn[2] = - (r1 + r2 + r3);
        eqn[1] = r1*r2 + r2*r3 + r3*r1;
        eqn[0] = - (r1 * r2 * r3);
        double res[] = new double[4];
        int n = solveCubicOld(eqn, res);
        System.out.print("old method returned "+n+" roots: ");
        for (int i = 0; i < n; i++) {
            System.out.print(res[i]+", ");
        }
        System.out.println();
        n = solveCubicNew(eqn, res);
        System.out.print("new method returned "+n+" roots: ");
        for (int i = 0; i < n; i++) {
            System.out.print(res[i]+", ");
        }
        System.out.println();
        System.out.println();
    }

    public static void trySolve2(double r1, double r2) {
        trySolve3(r1, r2, r2);
        trySolve3(r1, r1, r2);
    }

    public static void trySolve1(double r) {
        trySolve3(r, r, r);
    }

    public static void trySolve(double r1, double r2, double r3) {
        trySolve3(r1, r2, r3);
        trySolve2(r1, r2);
        trySolve2(r2, r3);
        trySolve2(r3, r1);
        trySolve1(r1);
        trySolve1(r2);
        trySolve1(r3);
    }

    public static void main(String argv[]) {
        trySolve(1, 2, 3);
        trySolve(15, 103, 27);
    }
}