16

I need to solve a cubic equation (ax^3 + bx^2 + c*x + d = 0) analytically and in real numbers, preferably in pure javascript (no libs). As there could be 1 to 3 roots, I think an array of numbers is a reasonable result type.

P.S. Provided my own solution below, hope it'll be useful.

Alexander Shtuchkin
  • 1,797
  • 15
  • 20

1 Answers1

30

Here you go. Includes handling degenerate cases. Main algorithm is mostly from wikipedia article.

function cuberoot(x) {
    var y = Math.pow(Math.abs(x), 1/3);
    return x < 0 ? -y : y;
}

function solveCubic(a, b, c, d) {
    if (Math.abs(a) < 1e-8) { // Quadratic case, ax^2+bx+c=0
        a = b; b = c; c = d;
        if (Math.abs(a) < 1e-8) { // Linear case, ax+b=0
            a = b; b = c;
            if (Math.abs(a) < 1e-8) // Degenerate case
                return [];
            return [-b/a];
        }

        var D = b*b - 4*a*c;
        if (Math.abs(D) < 1e-8)
            return [-b/(2*a)];
        else if (D > 0)
            return [(-b+Math.sqrt(D))/(2*a), (-b-Math.sqrt(D))/(2*a)];
        return [];
    }

    // Convert to depressed cubic t^3+pt+q = 0 (subst x = t - b/3a)
    var p = (3*a*c - b*b)/(3*a*a);
    var q = (2*b*b*b - 9*a*b*c + 27*a*a*d)/(27*a*a*a);
    var roots;

    if (Math.abs(p) < 1e-8) { // p = 0 -> t^3 = -q -> t = -q^1/3
        roots = [cuberoot(-q)];
    } else if (Math.abs(q) < 1e-8) { // q = 0 -> t^3 + pt = 0 -> t(t^2+p)=0
        roots = [0].concat(p < 0 ? [Math.sqrt(-p), -Math.sqrt(-p)] : []);
    } else {
        var D = q*q/4 + p*p*p/27;
        if (Math.abs(D) < 1e-8) {       // D = 0 -> two roots
            roots = [-1.5*q/p, 3*q/p];
        } else if (D > 0) {             // Only one real root
            var u = cuberoot(-q/2 - Math.sqrt(D));
            roots = [u - p/(3*u)];
        } else {                        // D < 0, three roots, but needs to use complex numbers/trigonometric solution
            var u = 2*Math.sqrt(-p/3);
            var t = Math.acos(3*q/p/u)/3;  // D < 0 implies p < 0 and acos argument in [-1..1]
            var k = 2*Math.PI/3;
            roots = [u*Math.cos(t), u*Math.cos(t-k), u*Math.cos(t-2*k)];
        }
    }

    // Convert back from depressed cubic
    for (var i = 0; i < roots.length; i++)
        roots[i] -= b/(3*a);

    return roots;
}
Thomas Gales
  • 107
  • 2
  • 8
Alexander Shtuchkin
  • 1,797
  • 15
  • 20
  • 1
    Props for writing that, assuming it works! I might suggest you edit the question; I think it is attracting downvotes because it seems to be demanding an answer with no work, when in fact you're providing the answer here and you've done the work! – Phil H Nov 28 '14 at 09:38
  • Hardcoding epsilons is a bad practice. – Stepan Yakovenko Nov 02 '15 at 06:10
  • Thank you very much. I ported it to php and it is working perfectly. The best solution I found. – miile7 Mar 28 '19 at 09:57
  • why are you using <1e-8 instead of === 0? Is it in order to avoid floating point error? – martian17 Jun 06 '20 at 17:05
  • yes, to avoid rounding errors. Please feel free to adjust epsilons as needed for your application. – Alexander Shtuchkin Jun 07 '20 at 02:59
  • 1
    @AlexanderShtuchkin thanks, it was super useful! The only remark - JS has a built-in cube root method (Math.cbrt), no need to create a custom one: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/cbrt – Andrew Simontsev Jul 02 '20 at 10:32