Solving cubic polynomials

There are many sources for how to solve cubic polynomials, but a summary is helpful when implementing a new code. This page is an attempt to record that.

We will consider the polynomial

<math> x^3 + a x^2 + b x + c = 0 </math>

First, calculate the coefficients of the 'depressed cubic' <math>y^3 + y x + y = 0</math> obtained by the substitution <math>y = x - a/3</math>:

<math>p = \frac{3 b - a^2}{3}</math>

<math>q = \frac{9 a b - 27 c - 2 a^3}{27}</math>

This can then be converted into a quadratic in <math>z^3</math> via the 'magic' Vieta substitution <math>y = z - \frac{p}{3 z}</math>,

<math>\left(z^3\right)^2 - q z^3 - \left(\frac{p}{3}\right)^3 = 0</math>

for which the solution is

<math>z^3 = \frac{q}{2} \pm \sqrt{\left(\frac{q}[2}\right)^2 + \left(\frac{p}{3})\right)}