Although a closed form solution exists for the roots of polynomials of degree ≤ 4, the general formulae for cubics (and quartics) is ugly. Various simplifications can be made; commonly, the cubic \(a_3\,x^3+a_2\,x^2+a_1\,x+a_0\) is transformed by substituting \(x = t-a_2/3a_3\), giving

\(\displaystyle t^3+p\,t+q=0\),

where

\(\begin{align*} p&=\frac{1}{a_{3}}\,\left(a_{1}-\frac{a_{2}^2}{3\,a_{3}}\right)\\ q&=\frac{1}{a_{3}}\,\left(a_{0}-\frac{a_{2}\,a_{1}}{3\,a_{3}}+\frac{2\,a_{2}^3}{27\,a_{3}^2}\right). \end{align*}\)

Several strategies exist for solving this depressed cubic, but most involve dealing with complex numbers, even when we only want to find real roots. Fortunately, there are a couple ways to find solutions with only real arithmetic. First, we need to know the value of the discriminant \(\Delta=4\,p^3+27\,q^2\).

If \(\Delta<0\), then there are three real roots; this is the tough case. Luckily, we can use the uncommon trigonometric identity \(4\,\cos^3(\theta)-3\,\cos(\theta)-cos(3\,\theta)=0\) to calculate the roots with trig functions. Making another substitution, \(t=2\,\sqrt{-p/3}\,\cos(\theta)\), gives us

\(\displaystyle 4\,\cos^3(\theta)-3\,\cos(\theta)-\frac{3\,q}{2\,p}\,\sqrt{\frac{-3}{p}}=0\),

so

\(\displaystyle \cos(3\,\theta)=\frac{3\,q}{2\,p}\,\sqrt{\frac{-3}{p}}\).

Thus the three roots are:

\(\displaystyle t_n=2\,\sqrt{\frac{-p}{3}}\,\cos\left(\frac{1}{3}\,\cos^{-1}\left(\frac{3\,q}{2\,p}\,\sqrt{\frac{-3}{p}}\right)-\frac{2\,\pi\,n}{3}\right)\), for \(n\in\lbrace0,1,2\rbrace\).

We can save a cosine calculation by noting that \(t_2=-t_0|_{q=-q}\) and \(t_1=-(t_0+t_2)\). When generated in this way, the roots are ordered from smallest to largest (\(t_0\) ≤ \(t_1\) ≤ \(t_2\)).

For the \(\Delta>0\) case, there is a single real root. We can apply the general cubic formula and avoid dealing with complex numbers, provided we choose the signs right. In this case,

\(\displaystyle t=-\frac{|q|}{q}\,\left(C-\frac{p}{3\,C}\right)\), where \(\displaystyle C=\sqrt[3]{\sqrt{\frac{\Delta}{108}}+\frac{|q|}{2}}\).

When \(\Delta=0\), there is a root of at least multiplicity 2, and we can avoid calculating any radicals. If \(p=0\), the only root is \(t_0=0\); otherwise, the duplicate root is \(-3\,q/2\,p\) and the simple root is \(3\,q/p\).

My C implementation of this solution method can be seen in the dmnsn_solve_cubic() function here.