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

Question

I need to compute the integral of the following function within ranges that start as low as -150:

import numpy as np
from scipy.special import ndtr

def my_func(x):
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))

The problem is that this part of the function

np.exp(x ** 2)

tends toward infinity -- I get inf for values of x less than approximately -26.

And this part of the function

2 * ndtr(x * np.sqrt(2))

which is equivalent to

from scipy.special import erf

1 + erf(x)

tends toward 0.

So, a very, very large number times a very, very small number should give me a reasonably sized number -- but, instead of that, python is giving me nan.

What can I do to circumvent this problem?

Answer 1

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))
# nan

print(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 mpmath

def 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))
# inf

print(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

Answer 2

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

Answer 3

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 , which you correctly identified would be . Opening the brackets you can integrate both parts of the summation.

Scipy has this imaginary error function implemented

The second part is harder:

This is a generalized hypergeometric function. Sadly it looks like scipy does not have an implementation of it, but this package claims it does.

Here I used indefinite integrals without constants, knowing the from to values it is clear how to use definite ones.