如何使用scipy和超大数字

Question

我想使用scipy的特殊函数来处理非常大的数字

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

math.exp抛出溢出错误，这并不奇怪。所以我试着用十进制模块代替

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)

但出现以下错误：

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来自scipy

正确的处理方法是什么？你知道吗

Answer 1

一种选择是使用^{}。它包括^{}的实现。你知道吗

例如

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、C++或FORTRAN实现的，并且限制在64位浮点值上，最大值约为1.8E308:

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