# 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. – march Oct 12 '17 at 22:27
• One thing to do is plug in the two sets to see which, if either, give residuals near zero. – Daniel Lichtblau Oct 12 '17 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}

• Thanks to everyone, now I understand! – TeM Oct 13 '17 at 5:56

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] &.... – Michael E2 Oct 13 '17 at 1:11
• Oh, right. That's what was going on. – march Oct 13 '17 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.