To Solve a Cubic Equation

• Download Cardano.zip - 62.1 KB

Cardano's Method

Gerolamo Cardano published a method to solve a cubic equation in 1545. There is a description of this method on Wikipedia. But it is not too detailed and on the German Wikipedia. In fact, the last part is missing and without this part, one cannot implement it into an algorithm. So I thought I could try to pick up there where the Wikipedia description ends :)

The Wikipedia description starts with the qubic equation

And even though some details are missing, the Wikipedia description is OK until the part:

  and  

Now things become interesting and most descriptions get vague.

There are the three possibilities: D >= 0 , D = 0 and D < 0.

D > 0

If D > 0 the root of D is real and we get one real solution and two complex solutions for u as well as for v:

The real solutions are

  and  

The two complex solutions are of the same length as u0 and v0, each one has the complex angle 120° for one and 240° for the other solution. That means with:

  and  

For the complex angles we get:

  and  

And:

  and  

So we get three possible solutions for u and three possible solutions for v and if we combine them, we get nine possible solutions. That’s pretty much it, but not all of them are valid. Furthermore, above we had the clause u*v = - p/3 and this is not the case for all nine combinations of u and v. It’s valid only for the combinations uo * vo, u1 * v2 and u2 * v1. Let’s have look at u1 * v2. That’s:

If we separate:

And:

That only works for ε1 * ε2. If we combine ε1 * ε1, there remains an imaginary part and that can’t give - p/3. This is a part that is most often missing in the description of Cardano’s formula:

So, given that D > 0 there are the three solutions. But only the first one is real.

y1 = u0 + v0

y2 = u1 + v2

y3 = u2 + v1

And on the way to this the substitution:

x = y – a/3 = u + v – a/3

has been made. This should not be forgotten as we are looking for "x" finally :)

If we put this into code it looks like this:

if (d > 0)
{
    u = Xroot(-q / 2.0 + Math.Sqrt(d), 3.0);
    v = Xroot(-q / 2.0 - Math.Sqrt(d), 3.0);
    x1.real = u + v - a / 3.0;
    x2.real = -(u + v) / 2.0 - a / 3.0;
    x2.imag = Math.Sqrt(3.0) / 2.0 * (u - v);
    x3.real = x2.real;
    x3.imag = -x2.imag;
}

D = 0

If D = 0, we have the same three solutions, but as D = 0, u0 = v0 and

And y3 = y2. That means we get only two solutions.

In code that’s:

if (d == 0)
{
    u = Xroot(-q / 2.0, 3.0);
    v = Xroot(-q / 2.0, 3.0);
    x1.real = u + v - a / 3.0;
    x2.real = -(u + v) / 2.0 - a / 3.0;
}

D < 0

The case D < 0 is the so called "Casus irreducibilis." Now the root of D becomes complex and the cubic roots for u and v become complex too.

  and  

To solve this complex root we interpret the complex number as number of the length r and complex angle α.

With the substitution from further above:

4((q/2)² + (p/3)³) = 4D

We get the same length for u and v:

The angle becomes for u:

The angle becomes for v:

To get the cubic root out of this now we have to calculate the cubic root of r and basically divide α by 3 that’s it. At least for the one first solution. The cubic root of a complex number has not just one but three solutions. The angle α we get once by multiplying α/3 by 3, twice by multiplying (2π + α)/3 by 3 and thirdly by multiplying (4π + α)/3 by 3. That makes three solutions for u and three solutions for v. They look like:

And as α_v = 2π - α

That would give nine possible solutions again. But as we are interested in real solutions, we can reduce that to three real solutions. For the three combinations

u0 + v2 u1 + v1 u2 + v0

the imaginary part becomes 0 and the solutions are

This is the most complicate part:

if (d < 0)
{
    r = Math.Sqrt(-p * p * p / 27.0);
    alpha = Math.Atan(Math.Sqrt(-d) / -q * 2.0);
    if (q > 0)                         // if q > 0 the angle becomes 2 *PI - alpha
        alpha = 2.0 * Math.PI - alpha;

    x1.real = Xroot(r, 3.0) * (Math.Cos((6.0 * Math.PI - alpha) / 3.0)
	    + Math.Cos(alpha / 3.0)) - a / 3.0;
    x2.real = Xroot(r, 3.0) * (Math.Cos((2.0 * Math.PI + alpha) / 3.0)
		+ Math.Cos((4.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
    x3.real = Xroot(r, 3.0) * (Math.Cos((4.0 * Math.PI + alpha) / 3.0)
		+ Math.Cos((2.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
}

But that’s all it takes :)