NumPy complex128除法与float64除法不一致答案

【问题标题】：NumPy complex128 division inconsistent with float64 divisionNumPy complex128除法与float64除法不一致
【发布时间】：2020-02-26 05:31:45
【问题描述】：

我在数值精度非常重要的上下文中进行复杂除法。我发现将两个没有虚部的 complex128 数字相除会得到超过 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 个十进制数字的 NumPy 浮点除法一致。我想以相同的精度使用 NumPy 复数除法。有没有办法做到这一点？

【问题讨论】：

标签： python numpy

【解决方案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

解释：

让我们写/// for complex by float 除法(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。

【讨论】：