public static int solveQuartic(double c[], double s[]) {
        /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */

        double A = c[3] / c[4];
        double B = c[2] / c[4];
        double C = c[1] / c[4];
        double D = c[0] / c[4];

        /*  substitute x = y - A/4 to eliminate cubic term:
        x^4 + px^2 + qx + r = 0 */

        double sq_A = A * A;
        double p = -3.0 / 8 * sq_A + B;
        double q = 1.0 / 8 * sq_A * A - 1.0 / 2 * A * B + C;
        double r = -3.0 / 256 * sq_A * sq_A + 1.0 / 16 * sq_A * B - 1.0 / 4 * A * C + D;

        double coeffs[] = new double[4];
        int num;
        if (isZero(r)) {
            /* no absolute term: y(y^3 + py + q) = 0 */
            coeffs[0] = q;
            coeffs[1] = p;
            coeffs[2] = 0;
            coeffs[3] = 1;

            num = solveCubic(coeffs, s);

            s[num++] = 0;
        } else {
            /* solve the resolvent cubic ... */
            coeffs[0] = 1.0 / 2 * r * p - 1.0 / 8 * q * q;
            coeffs[1] = -r;
            coeffs[2] = -1.0 / 2 * p;
            coeffs[3] = 1;

            solveCubic(coeffs, s);

		    /* ... and take the one real solution ... */

            double z = s[0];

		    /* ... to build two quadric equations */

            double u = z * z - r;
            double v = 2 * z - p;

            if (isZero(u))
                u = 0;
            else if (u > 0)
                u = sqrt(u);
            else
                return 0;

            if (isZero(v))
                v = 0;
            else if (v > 0)
                v = sqrt(v);
            else
                return 0;

            coeffs[0] = z - u;
            coeffs[1] = q < 0 ? -v : v;
            coeffs[2] = 1;

            num = solveQuadric(coeffs, s);

            coeffs[0] = z + u;
            coeffs[1] = q < 0 ? v : -v;
            coeffs[2] = 1;

            double s2[] = {s[num], s[num + 1]};
            num += solveQuadric(coeffs, s2);
        }

        /* resubstitute */

        double sub = 1.0 / 4 * A;

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

        return num;
    }