2

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 ?
B. M.
  • 18,243
  • 2
  • 35
  • 54

2 Answers2

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.

  • Thank's for this confirmation and suggestions. I guess there is a more direct and symmetric approach, but I don't see it at the moment. I edit my post in this goal. – B. M. Mar 06 '16 at 07:59
  • Final version here. http://stackoverflow.com/questions/35795663/fastest-way-to-find-the-smallest-positive-real-root-of-quartic-polynomial-4-degr/35829660#35829660 . I prefer avoid the division by `u` to limit rounding issues, but the precision is not very good. – B. M. Mar 06 '16 at 20:35
2

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thanks. It will nevertheless fails on `q=p=0`, which is not an rare case. – B. M. Mar 06 '16 at 17:45
  • Naturally, you have to check,catch and treat these special cases beforehand. People have tried robust implementations of Cardano, as I read somewhere, it gets rather lengthy. --`p==0 => z=(-q)**(1/3)*J^k`, `q==0 => z = 0, +-(-p)**(1/2)` etc. – Lutz Lehmann Mar 06 '16 at 18:39