Check if root of cubic equation is complex or not?

Question

I use this Cubic root implementation.

I have equation #1:

x³ -2 x² -5 x + 6 = 0

It gives me 3 complex roots ({real, imaginary}):

{-2, 7.4014868308343765E-17}
{1 , -2.9605947323337506E-16}
{3 , 2.9605947323337506E-16}

But in fact, the right result should be 3 non-complex roots: -2, 1, 3.

With this case, I can test by: apply 3 complex roots to the equation, it returns non-zero result (failed); apply 3 non-complex roots to the equation, it returns zero result (passed).

But there is the case where I apply both 3-complex roots and 3-non-complex roots to the equation (e.g. 47 x³ +7 x² -52 x + 0 = 0), it return non-zero (failed).

I think what causes this issue is because of this code:

/// <summary>
/// Evaluate all cubic roots of this <c>Complex</c>.
/// </summary>
public static (Complex, Complex, Complex) CubicRoots(this Complex complex)
{
    var r = Math.Pow(complex.Magnitude, 1d/3d);
    var theta = complex.Phase/3;
    const double shift = Constants.Pi2/3;
    return (Complex.FromPolarCoordinates(r, theta),
        Complex.FromPolarCoordinates(r, theta + shift),
        Complex.FromPolarCoordinates(r, theta - shift));
}

I know that floating point value can lose precision when calculating (~1E-15), but the problem is the imaginary part needs to decide weather it's zero or non-zero to tell if it's complex number or not.

I can't tell the user of my app: "hey user, if you see the imaginary part is close enough to 0, you can decide for yourself that the root's not a complex number".

Currently, I use this method to check:

const int TOLERATE = 15;
bool isRemoveImaginary = System.Math.Round(root.Imaginary, TOLERATE) == 0; //Remove imaginary if it's too close to zero

But I don't know if this method is appropriate, what if the TOLERATE = 15 is not enough. Or is it the right method to solve this problem?

So I want to ask, is there any better way to tell the root is complex or not?

Answer 1

Thank you Mark Dickinson.

So according to Wikipedia:

delta > 0: the cubic has three distinct real roots

delta < 0: the cubic has one real root and two non-real complex conjugate roots.

The delta D = (B*B - 4*A*A*A)/(-27*a*a)

My ideal is:

• delta > 0: remove all imaginary numbers of 3 roots.

• delta < 0: find the real root then remove its imaginary part if any (to make sure it's real). Leave the other 2 roots untouched. Now I have 2 ideas to find the real root:

Ideal #1

In theory, the real root should have imaginary = 0, but due to floating point precision, imaginary can deviate from 0 a little (e.g. imaginary = 1E-15 instead of 0). So the idea is: the 1 real root among 3 roots should have the imaginary whose value is closest to 0.

Code:

NumComplex[] arrRoot = { x1, x2, x3 };

if (delta > 0)
{
    for (var idxRoot = 0; idxRoot < arrRoot.Length; ++idxRoot)
        arrRoot[idxRoot] = arrRoot[idxRoot].RemoveImaginary();
}
else
{
    //The root with imaginary closest to 0 should be the real root,
    //the other two should be non-real.
    var realRootIdx = 0;
    var absClosest = double.MaxValue;
    double abs;

    for (var idxRoot = 0; idxRoot < arrRoot.Length; ++idxRoot)
    {
        abs = System.Math.Abs(arrRoot[idxRoot].GetImaginary());
        if (abs < absClosest)
        {
            absClosest = abs;
            realRootIdx = idxRoot;
        }
    }

    arrRoot[realRootIdx] = arrRoot[realRootIdx].RemoveImaginary();
}

The code above can be wrong if there are 3 roots ({real, imaginary}) like this:

{7, -1E-99}
{3, 1E-15}//1E-15 caused by floating point precision, 1E-15 should be 0
{7, 1E-99}//My code will mistake this because this is closer to 0 than 1E-15.

Maybe if that case does happen in real life, I will come up with a better way to pick the real root.

Idea #2

Take a look at how the 3 roots calculated:

x1 = FromPolarCoordinates(r, theta);
x2 = FromPolarCoordinates(r, theta + shift);
x3 = FromPolarCoordinates(r, theta - shift);

3 roots have the form (know this by tests, not proven by math):

x1 = { A }
x2 = { B, C }
x3 = { B, -C }

Use math knowledge to prove which one among the 3 roots is the real one.

Trial #1: Maybe the root x1 = FromPolarCoordinates(r, theta) is always real? (failed) untrue because the following case proved that guess is wrong: -53 x³ + 6 x² + 14 x - 54 = 0 (Thank Mark Dickinson again)

I don't know if math can prove something like: while delta < 0: if B < 0 then x3 is real, else x1 is real?

So until I get better idea, I'll just use idea #1.