# Third degree equation resolution

On page 17 of the Handbook of Mathematical Functions With Formulas, Graphs and Mathematical Tables by Abramowitz and Stegun, it is written that given the third degree equation:

$$z^3 + a_2\,z^2 + a_1\,z + a_0 = 0$$

and calculating the following values first:

$$q = \frac{1}{3}\,a_1 - \frac{1}{9}\,a_2^2\,, \quad r = \frac{1}{6}\left(a_1\,a_2 - 3\,a_0\right) - \frac{1}{27}\,a_2^3$$

and then these other values:

$$s_1 = \left(r + \left(q^3 + r^2\right)^{\frac{1}{2}}\right)^{\frac{1}{3}} \,, \quad s_2 = \left(r - \left(q^3 + r^2\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}\,,$$

the three desired roots of the cubic equation are:

\begin{align*} & z_1 = (s_1 + s_2) - \frac{a_2}{3}\,, \\ & z_2 = - \frac{1}{2}\,(s_1 + s_2) - \frac{a_2}{3} + \frac{\text{i}\,\sqrt{3}}{2}\,(s_1 - s_2)\,, \\ & z_3 = - \frac{1}{2}\,(s_1 + s_2) - \frac{a_2}{3} - \frac{\text{i}\,\sqrt{3}}{2}\,(s_1 - s_2)\,.\end{align*}

Unfortunately, writing the following code:

a2 = 0;
a1 = 1;
a0 = 1;

NSolve[z^3 + a2 z^2 + a1 z + a0 == 0, z]

q = 1/3 a1 - 1/9 a2^2;
r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;

s1 = (r + (q^3 + r^2)^(1/2))^(1/3);
s2 = (r - (q^3 + r^2)^(1/2))^(1/3);

z1 = (s1 + s2) - a2/3 // N
z2 = -1/2 (s1 + s2) - a2/3 + I Sqrt[3]/2 (s1 - s2) // N
z3 = -1/2 (s1 + s2) - a2/3 - I Sqrt[3]/2 (s1 - s2) // N


I get:

{{z -> -0.682328}, {z -> 0.341164 - 1.16154 I}, {z -> 0.341164 + 1.16154 I}}

0.835342 + 0.876227 I

0.341164 - 0.590913 I

-1.17651 - 0.285314 I


so there is definitely something that does not work, but I do not know what! Thank you!

• Maybe check out equations 22, 23, 49-56 on this page and compare. Commented Oct 12, 2017 at 22:27
• One thing to do is plug in the two sets to see which, if either, give residuals near zero. Commented Oct 12, 2017 at 23:19

As already noted by other answers, the version of Cardano in Abramowitz and Stegun uses a branch cut for the cube root that is not the same as the branch cut used by Mathematica.

So, here's how to implement Cardano manually for this case, as adapted from Numerical Recipes:

myCubic[ac_, bc_, cc_, dc_] :=
Module[{a, b, c, da, db, ds, dt, du, q, r},
{a, b, c} = {bc, cc, dc}/ac;
q = (a^2 - 3 b)/9; r = (2 a^3 + 9 (3 c - a b))/54;
ds = Sqrt[r^2 - q^3];
If[Re[Conjugate[r] ds] < 0, ds = -ds];
da = -(r + ds)^(1/3); db = If[da == 0, 0, q/da];
dt = da + db; du = I Sqrt[3] (da - db)/2;
{-dt/2 - du, dt, -dt/2 + du} - a/3]


Test:

With[{a = 1, b = 2, c = 3, d = 4},
FromDigits[{a, b, c, d}, #] & /@ N[myCubic[a, b, c, d]] // Chop]
{0, 0, 0}

With[{a = I - 3, b = 5 + 4 I, c = 2 - 2 I, d = 6},
FromDigits[{a, b, c, d}, #] & /@ N[myCubic[a, b, c, d]] // Chop]
{0, 0, 0}


The cube root (1/3 power) in A&S is CubeRoot in Mathematica:

a2 = 0;
a1 = 1;
a0 = 1;

NSolve[z^3 + a2 z^2 + a1 z + a0 == 0, z]

q = 1/3 a1 - 1/9 a2^2;
r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;

s1 = (r + (q^3 + r^2)^(1/2)) // CubeRoot;
s2 = (r - (q^3 + r^2)^(1/2)) // CubeRoot;

z1 = (s1 + s2) - a2/3 // N
z2 = -1/2 (s1 + s2) - a2/3 + I Sqrt[3]/2 (s1 - s2) // N
z3 = -1/2 (s1 + s2) - a2/3 - I Sqrt[3]/2 (s1 - s2) // N
(*
{{z -> -0.682328}, {z -> 0.341164 - 1.16154 I}, {z -> 0.341164 + 1.16154 I}}

-0.682328
0.341164 + 1.16154 I
0.341164 - 1.16154 I
*)

• Or Surd[#, 3] &.... Commented Oct 13, 2017 at 1:11
• Oh, right. That's what was going on. Commented Oct 13, 2017 at 2:55

Because of the cube roots, there are three, in general complex, solutions for s1 and s2. You have a positive discriminant, $q^3+r^2=31/108$, corresponding to one real root and two complex conjugate roots.

Test all possible solutions for s1 and s2 with something like the following:

Block[{a2 = 0, a1 = 1, a0 = 1, q, r, s1, s2, x},
q = 1/3 a1 - 1/9 a2^2;
r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;
Print["Discriminant: ", q^3 + r^2, " > 0"];
s1 = Chop[x /. {ToRules[Reduce[x^3 == r + (q^3 + r^2)^(1/2), x]]}];
s2 = Chop[x /. {ToRules[Reduce[x^3 == r - (q^3 + r^2)^(1/2), x]]}];
Flatten[
Table[{
Chop[N[(s1[[i]] + s2[[j]]) - a2/3]],
Chop[N[-1/2 (s1[[i]] + s2[[j]]) - a2/3 + I Sqrt[3]/2 (s1[[i]]-s2[[j]])]],
Chop[N[-1/2 (s1[[i]] + s2[[j]]) - a2/3 - I Sqrt[3]/2 (s1[[i]]-s2[[j]])]]},
{i, 1, 3}, {j, 1, 3}], 1]
] // Column


Pick the solution with one real root and two complex conjugate roots. If the discriminant is less than zero, then there are three real roots to pick. If the discriminant equals zero, then there are three real roots with at least two of them equal.