无法实现Cardano方法。复杂麻木的立方根

Question

为了提高np.roots三次方程的性能，我尝试实现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

当根是真实的时，效果很好：

^{pr2}$

但在复杂的情况下却失败了：

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))

{{cd2}的求值是{cd2}的问题。 Theory说uv=-p/3，但是这里uv=pJ/3：(u,v)不是一对好的根。在

在任何情况下，获得正确配对的最佳方法是什么？在

编辑

在@Sally post之后，我可以精确地解决这个问题。好的一对并不总是(u,v)，它可以是(u,vJ)或{}。所以问题可能是：

• 有没有一种直接的方法来抓住好的一对？在

终极：目前，通过使用Numba编译这段代码，它比np.roots快20倍。在

• 有没有更好的算法来计算三次方程的三次根？在

Answer 1

因为作为平方根之一的r的符号是自由的（resp。u和{}在u+v和{}作为二次多项式的根）集中是可互换的

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

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

这可以确保除数总是尽可能大，从而减少商的误差，并将被（接近）零除法可能成为问题的情况的数量减到最小。在

Answer 2

您正确地识别了这个问题：在一个复平面中有三个可能的立方根值，导致了9个可能的((-q+r)/2)**(1/3)和{}对。在这9对中，只有3对可以找到正确的根：即u*v=-p/3的根。一个简单的解决方法是将v的公式替换为v=-p/(3*u)。这可能也是一种加速：除法应该比求立方根快。在

但是u可能等于或接近于零，在这种情况下，除法变得可疑。实际上，在您的第一个例子中，它会使精度稍差一点。这里有一个健壮的数值方法：在return语句之前插入这两行。在

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

这个循环遍历v的三个候选者，选择使u*v+p/3的绝对值最小的一个。也许可以通过存储三个候选项来稍微提高性能，这样就不必重新计算胜出者。在