class QuarticRootFinder {

    public static ImmutableList<Complex> find(double a, double b, double c, double d, double e) {
        if (a == 0) {
            return CubicRootFinder.find(b, c, d, e);
        }

        //See https://en.wikipedia.org/wiki/Quartic_function#Solution_methods - subsection: Ferrari's Solution

        //Simplify equation so its of the form x^4 + bx^3 + cx^2 + dx + e = 0
        b = b / a;
        c = c / a;
        d = d / a;
        e = e / a;

        //Define y^4 + py^2 + qy + r = 0 (where x = y - b/4)
        Complex yToXAdjustment = new Complex(-b / 4, 0);
        double p = c - 3 * Math.pow(b, 2) / 8;
        double q = d + (Math.pow(b, 3) / 8) - (b * c / 2);
        double r = e - (3 * Math.pow(b, 4) / 256) + (Math.pow(b, 2) * c / 16) - (b * d / 4);

        //Special case of y^4 = 0
        if (p == 0 && q == 0 && r == 0) {
            return ImmutableList.of(yToXAdjustment, yToXAdjustment, yToXAdjustment, yToXAdjustment);
        }

        //Solve the following cubic for m
        ImmutableList<Complex> rootsOfCubic = CubicRootFinder.find(1,  p, (Math.pow(p, 2) / 4) - r, -Math.pow(q, 2) / 8);
        Complex m = rootsOfCubic.stream()
                .filter(complex -> Math.abs(complex.getImaginary()) >  EPSILON || Math.abs(complex.getReal()) > EPSILON)
                .findFirst()
                .orElseThrow(() -> new IllegalStateException("No non-zero solution for m found"));

        Complex yRoot1 = calculateRoot(p, q, m, true, true);
        Complex yRoot2 = calculateRoot(p, q, m, true, false);
        Complex yRoot3 = calculateRoot(p, q, m, false, true);
        Complex yRoot4 = calculateRoot(p, q, m, false, false);

        return ImmutableList.of(
                yRoot1.add(yToXAdjustment),
                yRoot2.add(yToXAdjustment),
                yRoot3.add(yToXAdjustment),
                yRoot4.add(yToXAdjustment));
    }

    private  QuarticRootFinder();

}