Python логарифм неправильно вычисляется очень малые комплексные числа

Question

мне нужно использовать Python для вычисления логарифмов объектов вида

log(1 - n0 - n1*1j) 

где n0 и n1 являются очень малые числа ~ 1.0e-27 и 1j является мнимое число.

Использование cmath.log дает неправильный ответ

print cmath.log(1.0 - 1e-27 - 1.0e-27j) 
(5e-55-1e-27j) 

Использование mpmath я могу получить правильный результат, но только если я фраза аргумент правильно

import mpmath as mp 
mp.prec = 100 
print mp.log( mp.mpc(1.0) - mp.mpc(0.00000000000000000000000001) - mp.mpc(real='0.0', imag='1.0e-27')) 

дает

(-1.0000389695486766657204483072e-26 - 1.00000000000000000000000001e-27j) 

(что является правильным ответом) тогда как

print mp.log( mp.mpc(0.99999999999999999999999999) - mp.mpc(real='0.0', imag='1.0e-27')) 

дает

(5.0e-55 - 1.0e-27j) 

Что здесь происходит? Могу ли я получить правильный ответ, используя cmath.log()?

Answer 1

Python использует стандарт IEEE binary64, обычно называемый двойной точностью, для представления чисел с плавающей запятой. binary64 имеет только около 15 десятичных цифр; цифры, которые занимают более 15 цифр, не могут быть точно представлены. Как вы можете увидеть здесь:

>>> 1 - 1e-27 
1.0 

потери точности вызывает 1 - 1e-27 округлить только 1.0. Более сложные математические библиотеки предлагают функции для вычисления логарифма на 1 больше, чем введенный вами номер. Например, в numpy:

>>> numpy.log1p(-1e-27-1e-27j) 
-1e-27j 

... который, к сожалению, также неправильно. Это меня удивило. Я не знаю, почему это произошло. Следующим лучшим вариантом, вероятно, будет использование математической библиотеки произвольной точности, такой как mpmath. Обязательно инициализируйте номера произвольной точности строками, а не плавает. Если вы используете литералы с плавающей запятой, вы потеряете тонну точности, прежде чем библиотека произвольной точности вступит в игру.

Answer 2

Двойной формат имеет точность примерно 15 десятичных цифр. cmath.log() теряет точность в вашем случае, так как он суммирует слишком разные значения величин (1. и ваше небольшое число).

Если вам не нужна высокая скорость в вашем алгоритме, вы можете использовать ряд Тейлора для log(1+x). Т.е .:

журнал (1 + х) = х - х ** 2/2 + х ** 3/3 + ...

Для еще большей точности можно реализовать Kahan summation algorithm.

Кроме того, можно использовать формулу:

журнал (1 + х) = math.log1p (х) + 1j * cmath.phase (1 + х),

, где math.log1p вычисляет логарифм именно для маленький x.