NumPy中complex128除法与float64除法不一致

Question

我正在进行复数除法运算，而数字精度确实很重要。我发现将两个没有虚数部分的complex128数相除，与将这两个数作为float64数相除，在小数点后15位后会得到不同的结果。

a = np.float64(1.501)
b = np.float64(1.337)

print('{:.20f}'.format(a / b))
# 1.12266267763649962852

a_com = np.complex128(1.501)
b_com = np.complex128(1.337)

print('{:.20f}'.format((a_com / b_com).real))
# 1.12266267763649940647

我有一个C++参考实现，其中复数除法与NumPy浮点数除法在15个小数位以上一致。我想使用具有相同精度的NumPy复数除法。是否有一种方法可以实现这一点？

- Jaroslaw M. Blusewicz

Answer 1

这似乎有效：

import numpy as np

def compl_div(A,B):
    A,B = np.asarray(A),np.asarray(B)
    Ba = np.abs(B)[...,None]
    A = (A[...,None].view(float)/Ba).view(complex)[...,0]
    B = (B.conj()[...,None].view(float)/Ba).view(complex)[...,0]
    return A*B

a = np.random.randn(10000)
b = np.random.randn(10000)
A = a.astype(complex)
B = b.astype(complex)

print((compl_div(A,B)==a/b).all())

print((np.sqrt(b*b)==np.abs(b)).all())

ac = a.view(complex)
bc = b.view(complex)

print(np.allclose(compl_div(ac,bc),ac/bc))

示例运行：

True    # complex without imag exactly equal float
True    # reason it works
True    # for nonzeron imag part do we actually get complex division

解释：

让我们用浮点除法的复数写成/// (x+iy)///r = x/r + iy/r

numpy似乎将复数除法A/B实现为A*(1/B)（可以将1/B计算为B.conj()///(B.conj()*B)），事实上，A/B似乎总是等于a*(1/b)

相反，我们使用(A///abs(B)) * (B.conj()///abs(B))作为abs(B)^2 = B*B.conj()这在数学上是等价的，但在数字上不等价。

现在，如果我们有abs(B) == abs(b)，那么A///abs(B) = a/abs(b)和B///abs(B) = sign(b)，我们可以看到compl_div(A,B)确实精确地返回a/b。

由于 abs(x+iy) = sqrt(x^2+y^2)，我们需要证明 sqrt(b*b) = abs(b)。除非在平方过程中出现上溢或下溢、平方为非规格化数或实现不符合IEEE标准，否则这是可证明的真实。

- Paul Panzer