Solving (real) cubic equations

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 +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.

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, 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.