Any cubic equation can be reduced to the form: find some x for
which x.x.x +A.x.x +B.x +C is zero

using no more than division by a
suitable scalar. By substituting y = x +A/3, so that y.y.y = x.x.x
+3.(A/3).x.x +3.(A.A/9).x +A.A.A/27, this can be reduced to the form y.y.y
−3.E.y −2.F = 0

for some E and F. (We can go further,
substituting u = y/√E or u = y/√−E to force E to be +1, 0 or
-1. This splits out three cases and gains us little in this general analysis,
but usually simplifies particular cases.)

The derivative of Q = (: y.y.y −3.E.y −2.F ←y :) gives
Q'(y) = 3.(y.y −E), which is zero when y is a square root of E; if E is
a negative real, this doesn't happen for any real y. Q's second derivative is
(: 6.y ←y :), so E's negative square root (if any) yields a local
maximum, while its positive square root yields a local minimum. Q's
derivative is large and positive for all large enough inputs; so we can infer
that Q is large and positive for all large enough positive inputs and large
and negative for all large enough negative inputs. In particular, this
implies at least one y for which Q(y) is zero (since Q has to get from
positive to negative *somewhere*).

I'm only considering the case where E and F are real. If E is zero, the equation is just y.y.y = 2.F and this has exactly one solution, with y as the cube root of 2.F. Hereafter, I'll assume E is non-zero.

At Q's local maximum, we have Q(−√E) = 2.(E.(√E) −F); if this is negative, i.e. if F>0 and E.E.E < F.F, then Q can only have one root (which must be >√E). At the local minimum, we have Q(√E) = −2.(E.(√E) +F); if this is positive, i.e. if F≤0 and E.E.E < F.F, then again Q can only have one root (which must be <−√E). If F = 0, we have Q(y) = y.(y.y −3.E), which has three zeros if E > 0, otherwise only one. Thus E.E.E < F.F implies Q has only one zero; we'll return to this case later.

If E.E.E = F.F (which, in particular, implies E>0), then F is ±E.√E and Q(−F/E) is zero; we then (because this is either a local maximum or a local minimum) have −F/E as a repeated root of Q, hence (y+F/E).(y+F/E) = y.y +2.y.F/E +E (as the square of F/E is E) as a factor of Q(y). Polynomial long division then implies Q's other factor is y −2.F/E; indeed

- (y.y +2.y.F/E +E).(y −2.F/E) = y.y.y −y.(4.F.F/E/E −E) −2.F

and F.F/E/E is E, so this is indeed Q(y). Thus, when E.E.E = F.F, Q has two zeros; a repeated one at −F/E and a simple one at 2.F/E.

If E.E.E > F.F (which, in particular, implies E>0), consider the possibility y = 2.(√E).Cos(a) for some angle a. Q(y) is thus

- Q(y) = (y.y −3.E).y −2.F = 2.E.(4.Cos(a).Cos(a) −3).(√E).Cos(a) −2.F

with F smaller than E.√E. Now, trigonometry tells us that Cos(3.a) = Cos(a).(4.Cos(a).Cos(a) −3), whence we infer that Q(y) = 2.E.(√E).Cos(3.a) −2.F is zero precisely when Cos(3.a) is F/E/√E, which we know lies between −1 and 1. From Q we know F and E, whence we can infer the possible values of 3.a, whence we can infer the possible values for Cos(a), each of which yields a value for 2.(√E).Cos(a) which must be a zero of Q.

Between 0.turn and turn/2, Cos takes all its values, so there must be exactly one a < turn/6 with Cos(3.a) = F/E/√E; the other candidates are then ±(a + n.turn/3) for integers n; taking Cos of this ignores the ± (since Cos(−t) = Cos(t) for any angle t) and gives us Cos(a +n.turn/3); since Cos is periodic with Cos(t+m.turn) = Cos(t) for any integer m and angle t, it suffices to consider n = 0, 1 and 2. Thus we get exactly three possible values for Cos(a +n.turn/3) and three possible values for y with Q(y) = 0. This is exactly how many we need.

For the remaining case, substitute y = p+q and note that y.y.y = p.p.p +q.q.q +3.p.p.q +3.p.q.q = p.p.p +q.q.q +3.p.q.(p+q), making y = p+q a solution of y.y.y = p.p.p +q.q.q +3.p.q.y; combining this with Q(y) = 0, i.e. y.y.y = 2.F +3.E.y, we obtain p.q = E and p.p.p +q.q.q = 2.F, from the first of which we infer p.p.p.q.q.q = E.E.E. This gives p.p.p and q.q.q as the inputs at which the quadratic (: t.t −2.F.t +E.E.E ←t :) yields zero output; such inputs are real iff F.F −E.E.E ≥ 0. The case we're left to consider is E.E.E < F.F, so yields real roots for this quadratic, namely F ±√(F.F −E.E.E). Taking the cube roots of these two, we get p and q; adding these we obtain the single value of y for which Q(y) = 0.

Some or all of the above make up a standard technique, attributed to Cardan
(because he published it, after wheedling it out of someone else to whom he'd
promised he'd keep it secret), for solving cubics. The
module `study.math.cardan`
in my `python`

package provides a
function, `cardan(u,s,i,c)`, which implements the above, returning ({0}:
((u.x +s).x +i).x +c ←x |) as a sequence.