Tricking numpy/python into representing very large and very small numbers

python numpy floating-point numbers scipy

I think a combination of @askewchan's solution and scipy.special.log_ndtr will do the trick:

from scipy.special import log_ndtr_log2 = np.log(2)_sqrt2 = np.sqrt(2)def my_func(x):    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))def my_func2(x):    return np.exp(x * x + _log2 + log_ndtr(x * _sqrt2))print(my_func(-150))# nanprint(my_func2(-150)# 0.0037611803122451198

For x <= -20, log_ndtr(x) uses a Taylor series expansion of the error function to iteratively compute the log CDF directly, which is much more numerically stable than simply taking log(ndtr(x)).

Update

As you mentioned in the comments, the exp can also overflow if x is sufficiently large. Whilst you could work around this using mpmath.exp, a simpler and faster method is to cast up to a np.longdouble which, on my machine, can represent values up to 1.189731495357231765e+4932:

import mpmathdef my_func3(x):    return mpmath.exp(x * x + _log2 + log_ndtr(x * _sqrt2))def my_func4(x):    return np.exp(np.float128(x * x + _log2 + log_ndtr(x * _sqrt2)))print(my_func2(50))# infprint(my_func3(50))# mpf('1.0895188633566085e+1086')print(my_func4(50))# 1.0895188633566084842e+1086%timeit my_func3(50)# The slowest run took 8.01 times longer than the fastest. This could mean that# an intermediate result is being cached  100000 loops, best of 3: 15.5 µs per# loop%timeit my_func4(50)# The slowest run took 11.11 times longer than the fastest. This could mean# that an intermediate result is being cached  100000 loops, best of 3: 2.9 µs# per loop

python numpy floating-point numbers scipy

There already is such a function: erfcx. I think erfcx(-x) should give you the integrand you want (note that 1+erf(x)=erfc(-x)).

python numpy floating-point numbers scipy

Not sure how helpful will this be, but here are a couple of thoughts that are too long for a comment.

You need to calculate the integral of $2 \cdot e^{x^2} \cdot f(\sqrt{2}x)$ , which you correctly identified would be $e^{x^2}*(1 + erf(x))$ . Opening the brackets you can integrate both parts of the summation.