Skip to content

Instantly share code, notes, and snippets.

@weepy
Last active March 16, 2022 18:18
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save weepy/6009631 to your computer and use it in GitHub Desktop.
Save weepy/6009631 to your computer and use it in GitHub Desktop.
Cubic Roots Solver Function in JS
function CubicSolve(a, b, c, d){
b /= a;
c /= a;
d /= a;
var discrim, q, r, dum1, s, t, term1, r13;
q = (3.0*c - (b*b))/9.0;
r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
r /= 54.0;
discrim = q*q*q + r*r;
var roots = [ {real: 0, i: 0}, {real: 0, i: 0}, {real: 0, i: 0} ]
term1 = (b/3.0);
if (discrim > 0) { // one root real, two are complex
s = r + Math.sqrt(discrim);
s = ((s < 0) ? -Math.pow(-s, (1.0/3.0)) : Math.pow(s, (1.0/3.0)));
t = r - Math.sqrt(discrim);
t = ((t < 0) ? -Math.pow(-t, (1.0/3.0)) : Math.pow(t, (1.0/3.0)));
roots[0].real = -term1 + s + t;
term1 += (s + t)/2.0;
roots[1].real = roots[2].real = -term1;
term1 = Math.sqrt(3.0)*(-t + s)/2;
roots[1].i = term1;
roots[2].i = -term1;
return roots;
} // End if (discrim > 0)
// The remaining options are all real
if (discrim == 0){ // All roots real, at least two are equal.
r13 = ((r < 0) ? -Math.pow(-r,(1.0/3.0)) : Math.pow(r,(1.0/3.0)));
roots[0].real = -term1 + 2.0*r13;
roots[2].real = roots[1].real = -(r13 + term1);
return roots;
} // End if (discrim == 0)
// Only option left is that all roots are real and unequal (to get here, q < 0)
q = -q;
dum1 = q*q*q;
dum1 = Math.acos(r/Math.sqrt(dum1));
r13 = 2.0*Math.sqrt(q);
roots[0].real = -term1 + r13*Math.cos(dum1/3.0);
roots[1].real = -term1 + r13*Math.cos((dum1 + 2.0*Math.PI)/3.0);
roots[2].real = -term1 + r13*Math.cos((dum1 + 4.0*Math.PI)/3.0);
return roots;
}
@martian17
Copy link

Great work. One problem is that this does not work when a===0. This can be avoided by using quadratic and linear solver when a===0 or b===0, I think.

@weepy
Copy link
Author

weepy commented Jun 19, 2020

good point!

@hkhoont
Copy link

hkhoont commented Feb 18, 2021

Line 27 should be

roots[2].real = roots[1].real = -term1;

@weepy
Copy link
Author

weepy commented Feb 18, 2021

Good catch

@joshi1983
Copy link

joshi1983 commented Mar 16, 2022

Math.cbrt(x) will work more efficiently and correctly than Math.pow(x, 1/3) in some cases so this code would be better with that change.

Math.cbrt(-8) returns -2 but Math.pow(-8, 1/3) returns NaN for example. The -2 is a more desirable result.

The browser support for Math.cbrt is very good now. All major browsers support Math.cbrt including Chrome, Firefox, Edge, Safari. Internet Explorer 11 doesn't but that browser isn't at all popular now.

@weepy
Copy link
Author

weepy commented Mar 16, 2022

Yes that's right Math.pow doesn't work well for 1/3 for negative numbers.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment