How to work with scipy and extremely large numbers

Question

I would like to use special scipy functions to work with extremely large numbers

alpha = 9999
def y(t):
    return 1 / (special.lambertw(alpha * math.exp(alpha-t)) + 1)

math.expcauses an overflow error, which is not surprising. So I tried using a decimal module instead

alpha = 9999
def y(t):
    exp = decimal.Decimal(math.exp(1))
    exp = exp ** alpha
    exp = exp * decimal.Decimal(math.exp(-t))
    return 1 / (special.lambertw(alpha * math.exp(alpha-t)) + 1)

But get the following error:

TypeError: ufunc '_lambertw' not supported for the input types, and the 
inputs could not be safely coerced to any supported types according to 
the casting rule ''safe''

special.lambertw comes from scipy

What would be the right way to handle this?

Answer 1

One option is to use mpmath. It includes implementation lambertw.

For instance,

In [20]: import mpmath

In [21]: mpmath.mp.dps = 30

In [22]: alpha = 9999

In [23]: def y(t):
    ...:     return 1 / (mpmath.lambertw(alpha * mpmath.exp(alpha-t)) + 1)
    ...: 

In [24]: y(1.5)
Out[24]: mpf('0.000100015000749774938119797735952206')

, scipy . scipy- C, ++ Fortran 64- 1,8e308:

In [11]: np.finfo(np.float64).max
Out[11]: 1.7976931348623157e+308