未能实施 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

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

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

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

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

我想问题出在 u 和 v 的立方根求值上。 Theory 说 uv=-p/3，但这里 uv=pJ/3: (u,v) 不是一对好的根。

在所有情况下获得正确配对的最佳方法是什么？

编辑

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

有没有简单的方法来捕捉好对？

最终：目前，通过使用 Numba 编译此代码，它比 np.roots 快 20 倍。

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

Answer 1

您正确地识别了问题：复平面上的立方根有 3 个可能的值，因此有 9 对可能的 ((-q+r)/2)**(1/3) 和 ((-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 的绝对值最小的那个。也许可以通过存储三个候选者来略微提高性能，这样就不必重新计算获胜者。

Answer 2

因为作为平方根之一的 r 的符号是自由的（分别 u 和 v 的作用在 u+v 和 u^3,v^3作为二次多项式的根）集合

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

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

这可确保除数始终尽可能大，从而减少商的误差并最大限度地减少被（接近）零除可能成为问题的情况。

未能实施 Cardano 方法。复数的立方根

Fail to implement Cardano method. Cube root of a complex number

python

algorithm

math

numpy

numba