Failed to implement Cardano method. The cubic root of the complex number

Question

To improve the performance of np.roots according to the cubic equation, I am trying to implement the Cardan (o) method :

 def cardan(a,b,c,d): #"resolve P=ax^3+bx^2+cx+d=0" #"x=zb/a/3=z-z0 => P=z^3+pz+q" z0=b/3/a a2,b2 = a*a,b*bp=-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=((-qr)/2)**(1/3) return u+v-z0,u*J+v/J-z0,u/J+v*J-z0 

Works well when the 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 a crash in a difficult 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 assume that the problem is estimating u and v roots of the cube. The theory says uv=-p/3 , but here uv=pJ/3 : (u,v) not a good pair of roots.

What is the best way to get the right pair in all cases?

EDIT

After posting @Sally, I can pinpoint the problem. A good couple is not always (u,v) , it can be (u,vJ) or (uJ,v) . So the question may be as follows:

• Is there an easy way to catch a good pair?

And final: at the moment, compiling this code with Numba , it is 20 times faster than np.roots .

• Is there a better algorithm for calculating the three roots of a cubic equation?

Answer 1

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

However, u may be equal to or close to zero, in which case division becomes suspicious. Indeed, already in your first example, this makes accuracy a little worse. Here's a numerically reliable 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 crosses three candidates for v, choosing depending on what minimizes the absolute value of u*v+p/3 . Perhaps you can improve performance a bit here by retaining the three candidates so that the winner should not be recounted.

Answer 2

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

 u3 = (abs(q+r)>abs(qr))? -(q+r) : -(qr) u = u3**(1/3) v = -p/(3*u) 

This ensures that the divisor is always as large as possible, reducing the error in the factor and minimizing the number of cases when dividing by (almost) zero can be a problem.