Walking Randomly

I’ll trust the 709.xxxx answers. I tried an alternative calculation using the Decimal module in python:

>>> from decimal import *
>>> getcontext().prec = 400
>>> (Decimal(1.7976931348622732e+308)+Decimal(1)).ln()

The answer was: Decimal(‘709.7827128933839730844729653945406010028483308408967416882624245425182373438217354226633292092126001554896893194841985716513198090270719058674256471413044766448449620789865561439964561173965004395982567690334278537771153659333034417666750701631982343807434305084955327863646409139369174559366942984521874728821971739604744878331189497294407403384339452268952762722965577809095518932592274010897959168’)

Using base 10 logarithms and ignoring the “1p” allowed me to simplify the equation to something that doesn’t challenge the floating point limits:

>>> (308+np.log10(1.7976931348622732)) / np.log10(np.e)

Giving an answer of 709.78271289338397, which is exactly what you got in Linux and Mac OS X.

An alternative implementation, scipy.special.log1p, gives 709.78271289338397 on both Windows and Linux.

On Windows MSVC compiler maps long double to double, while *nix supports 128bit float:

http://stackoverflow.com/a/5022356/2230844

I don’t have a Windows machine, but my guess is that numpy on Windows falls back to using the naive implementation in https://github.com/numpy/numpy/blob/master/numpy/core/src/npymath/npy_math.c.src#L98, which computes log1p(x) as log(1 + x) times some compensated summation factors. The first multiplication pushes the result above the floating point cutoff, and it’s all infs from then on. On *nix, npy_log1p probably just uses log1p from the C standard library, which works fine.

To verify this, I created a small C program that links in both implementations (see https://gist.github.com/jvkersch/686051b6f7331c44add9). The native implementation in numpy (on my Macbook Pro) indeed prints out inf, while the other code path returns the 709… result

As discussed with you this morning, the problem is that on Windows numpy is using its own fallback. Presumably Windows doesn’t provide a system log1p, but Mac and Linux typically do.

numpy’s log1p is here:

https://github.com/numpy/numpy/blob/467d4e16d77a2e7c131aac53c639e82b754578c7/numpy/core/src/npymath/npy_math.c.src#L98

(for arguments “close” to the floating point maximum, `log(1+x) * x` overflows.)

As x gets large log(x) and log1p(x) get very close; the difference between them grows very small. Eventually both log(x) and log1p(x) are approximated by the same floating point number. This seems to happen at about 2**49.

observation: Since Python 2.6 python itself has a math.log1p function. So numpy could fall back to using that.

We could change numpy’s fallback so that it has an additional case for large x:

if(x > 0x1p49) {
return npy_log(x);
}

Oh yukh, I suppose the actual threshold for when log(x) and log1p(x) are the same depends on the floating point precision. But that’s numpy’s problem.

Hi,

Any reason why you would want to use log1p for such large arguments? Typically, you use log1p for problems where you want to calculate the logarithm for a number which is very close (IEEE double precision close) to 1.0.

I understand the compatibility argument, though. But you’re really way beyond the applicability range of log1p. As David Jones says, log(x) and log1p(x) for very large numbers will be the same (read: floating point-wise the same).