Fail to implement Cardano method. Cube root of a c

Question

In order to improve np.roots performance on cubic equation, I try to implement Cardan(o) Method :

def cardan(a,b,c,d):
    #"resolve P=ax^3+bx^2+cx+d=0" 
    #"x=z-b/a/3=z-z0 => P=z^3+pz+q"
    z0=b/3/a
    a2,b2 = a*a,b*b    
    p=-b2/3/a2 +c/a
    q=(b/27*(2*b2/a2-9*c/a)+d)/a
    D=-4*p*p*p-27*q*q+0j
    r=sqrt(-D/27)
    J=-0.5+0.86602540378443871j # exp(2i*pi/3)
    u=((-q+r)/2)**(1/3)
    v=((-q-r)/2)**(1/3)
    return u+v-z0,u*J+v/J-z0,u/J+v*J-z0

It works well when roots are real:

In [2]: P=poly1d([1,2,3],True)
In [3]: roots(P)
Out[3]: array([ 3.,  2.,  1.])
In [4]: cardan(*P)
Out[4]: ((3+0j), (1+0j), (2+1.110e-16j))

But fails in the complex case :

In [8]: P=poly1d([1,-1j,1j],True)
In [9]: P
Out[9]: poly1d([ 1., -1.,  1., -1.])
In [10]: roots(P)
Out[10]: array([  1.0000e+00+0.j,   7.771e-16+1.j,   7.771e-16-1.j])
In [11]: cardan(*P)
Out[11]: ((1.366+0.211j),(5.551e-17+0.577j),(-0.366-0.788j))

I guess that the problem is the evaluation of u and v by cube roots . Theory say uv=-p/3, but here uv=pJ/3: (u,v) is not a good pair of roots.

What is the best way to obtain a correct pair in all cases ?

EDIT

After @Sally post, I can precise the problem. The good pair is not always (u,v), it can be (u,vJ) or (uJ,v). So the question could be :

• is there a straightforward method to catch the good pair ?

And ultimate : for the moment, by compiling this code with Numba, it's 20x faster than np.roots .

• is there a better algorithm to compute the three roots of a cubic equation ?

Answer 1

Since the sign of r as one of the square roots is free (resp. the role of u and v is interchangeable in u+v and of u^3,v^3 as roots of a quadratic polynomial) set

u3 = (abs(q+r)>abs(q-r))? -(q+r) : -(q-r)

u = u3**(1/3)
v = -p/(3*u)

This makes sure that the divisor is always as large as possible, reducing the error in the quotient and minimizing the number of cases where division by (near) zero might become an issue.

Answer 2

You correctly identified the problem: there are three possible values of cubic root in a complex plane, resulting in 9 possible pairs of ((-q+r)/2)**(1/3) and ((-q-r)/2)**(1/3). Of these 9, only 3 pairs lead to the correct roots: namely, those for which u*v = -p/3. An easy fix is to replace the formula for v with v=-p/(3*u). This is probably also a speed-up: division should be faster than taking cubic root.

However u might be equal or close to zero, in which case the division becomes suspect. Indeed, already in your first example it makes the precision slightly worse. Here is a numerically robust approach: insert these two lines before the return statement.

choices = [abs(u*v*J**k+p/3) for k in range(3)]
v = v*J**choices.index(min(choices))

This loops over the three candidates for v, picking whichever minimizes the absolute value of u*v+p/3. Perhaps one can slightly improve the performance here by storing the three candidates so that the winner does not have to be recomputed.