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
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
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.
Written by Eddy.