NumPy complex128 division inconsistent with float64 division

Question

I'm doing complex division in a context where numerical precision really matters. I find that dividing two complex128 numbers with no imaginary part gives a different result beyond 15 decimal digits that dividing the same two numbers as 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

I have a C++ reference implementation where complex division agrees with NumPy float division beyond 15 decimal digits. I'd like to use NumPy complex division with the same precision. Is there a way to accomplish that?

Paul Panzer · Accepted Answer · 2019-10-31T03:21:05.770

This seems to work:

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

Sample run:

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

Explanation:

Let us write /// for complex by float division (x+iy)///r = x/r + iy/r

numpy seems to implement complex division A/B as A*(1/B) (1/B can be computed as B.conj()///(B.conj()*B)), indeed A/B appears to always equal a*(1/b)

We do instead (A///abs(B)) * (B.conj()///abs(B)) as abs(B)^2 = B*B.conj() this is mathematically, but not numerically, equivalent.

Now, if we had abs(B) == abs(b) then A///abs(B) = a/abs(b) and B///abs(B) = sign(b) and we could see that compl_div(A,B) indeed gives back exactly a/b.

As abs(x+iy) = sqrt(x^2+y^2) we need to show sqrt(b*b) = abs(b). This is provably true unless there is over or underflow in the square or the square is denormal or the implementation does not conform to IEEE.

NumPy complex128 division inconsistent with float64 division

1 Answers1