NumPy complex128 разделение несовместимо с разделением float64

Question

Я делаю сложное деление в контексте, где численная точность действительно имеет значение. Я обнаружил, что деление двух комплексных чисел 128 без мнимой части дает другой результат, кроме 15 десятичных цифр, который разделяет те же два числа, что и float64. деление после 15 десятичных цифр. Я бы хотел использовать комплексное деление 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). Это доказуемо верно , если в квадрате нет переполнения или недостатка, или квадрат не является нормальным или реализация не соответствует IEEE.