2

I need to use python to calculate the logarithms of objects of the form

log( 1 - n0 - n1*1j)

where n0 and n1 are very small numbers ~ 1.0e-27 and 1j is the the imaginary number.

Using cmath.log gives the wrong answer

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

Using mpmath I can obtain the correct result but only if I phrase the argument correctly

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

gives

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

(which is the correct answer) whereas

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

gives

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

What's going on here? Can I get the right answer just using cmath.log()?

egoburnswell
  • 77
  • 7
  • 1
    The "correct answer" seems to me to be off by an order of magnitude in the real part. [Wolfram Alpha](http://www.wolframalpha.com/input/?i=log1p%28-1e-27-1e-27i%29) and my intuition both suggest that the real part should be closer to `-1e-27`. – user2357112 Jan 16 '14 at 11:37
  • It is just larger than 1 if it was just smaller than one it would be something like 9.99e-27, I think. – egoburnswell Jan 16 '14 at 11:43
  • 2
    `0.00000000000000000000000001` is `1e-26`, not `1e-27`. Why didn't you initialize it with a string, like `mp.mpc(real='0.0', imag='1.0e-27')`? – user2357112 Jan 16 '14 at 11:46
  • You're right, if initialized as 'mp.log(mp.mpc(real='0.99999999999999999999999999', imag='-1.0e-27') )' it does give the correct answer. – egoburnswell Jan 16 '14 at 11:55

2 Answers2

3

Python uses the IEEE binary64 standard, commonly referred to as double precision, to represent floating-point numbers. binary64 only has about 15 decimal digits worth of precision; numbers that take more than 15 digits can't be accurately represented. As you can see here:

>>> 1 - 1e-27
1.0

Precision loss causes 1 - 1e-27 to round to just 1.0. More sophisticated math libraries offer functions to compute the logarithm of 1 more than the number you input. For example, in numpy:

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

...which unfortunately is also wrong. That surprised me. I'm not sure why that happened. The next best bet would probably be to use an arbitrary-precision math library like mpmath. Make sure to initialize your arbitrary-precision numbers with strings, rather than floats. If you use floating-point literals, you'll lose a ton of precision before the arbitrary-precision library comes into play.

user2357112
  • 260,549
  • 28
  • 431
  • 505
1

Double format has precision roughly 15 decimal digits. cmath.log() loses precision in your case since it sums up too different in magnitude numbers (1. and your small number).

If you do not need high speed in your algorithm, you can use Taylor series for log(1+x). I.e.:

log(1+x) = x - x**2/2 + x**3/3 + ...

For even more precision you can implement Kahan summation algorithm.

You can also use the formula:

log(1+x) = math.log1p(x) + 1j*cmath.phase(1+x),

where math.log1p calculates logarithm precisely for small x.

unixmin
  • 163
  • 7
  • Here, since `x ~ 1e-26`, first order approximation `log(1+x) = x` gives precision up to `~1e-52`, which is probably enough. – aland Jan 16 '14 at 12:04