NumPy complex128除法与float64除法不一致

Question

我是在一个数字精度非常重要的背景下做复杂除法的。我发现将两个没有虚数部分的复数128除以15个十进制数字，得到的结果与将两个数除以float64得到的结果不同。你知道吗

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位数以上的Nuffy浮点除法一致。我想用同样精度的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，除非平方中有过流或下流，或者平方是非规范的，或者实现不符合IEEE。你知道吗