python - 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++ 引用实现，其中复数除法与超过 15 位小数的 NumPy float 除法一致。我想以相同的精度使用 NumPy 复数除法。有办法实现吗？

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) .这是 provably true除非方 block 出现上溢或下溢，或者方 block 不规范，或者实现不符合IEEE。