如何处理 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

一种选择是使用 mpmath. It includes an implementation of lambertw。

例如，

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

如何处理 scipy 和非常大的数字

How to work with scipy and extremely large numbers

python

scientific-computing

scipy