评估 mpmath gammainc 函数的困难

Question

我正在使用 python 库 mpmath，特别是评估不完整的 gamma 函数。 这是求根例程的一部分，但对于复值参数的某些组合，它的计算速度非常慢。

import mpmath as mp
from mpmath import gammainc
def Gamma(a,z0,z1):
    return gammainc(a,a=z0,b=z1,regularized=False)

mpmath.gammainc 函数的求值在这里卡住了：

>> Gamma(mp.mpc(12.5+17.5j), mp.mpf(0.0), mp.mpf(-12.5))

另一方面，Mathematica 几乎立即返回结果：

In[1]:= Gamma[12.5 + 17.5 I, 0, -12.5]
Out[1]:= 2.38012*10^-7 + 5.54827*10^-7 I

在其他情况下，对于不同的参数 mpmath 和 Mathematica 返回相同的输出：

数学

In[2]: Gamma[3.5 I, 0, 10]
Out[2]:= 0.0054741 + 0.000409846 I

Python mpmath

>> Gamma(3.5j,0,10)
mpc(real='0.0054741038497352953', imag='0.00040984640431357779')

你知道这种行为的原因吗？这可以被认为是mpmath 的问题还是正交的数学问题？ 不幸的是，scipy 没有为复杂参数提供gamma 函数的实现，所以它不是一个选项。

Answer 1

显然，在某些情况下，mpmath 中的一个错误会导致它在评估 gammainc 时进入无限循环。值得报告mpmath tracker。但至少对于您提到的情况，一种解决方法是将三参数不完全 gamma 函数编写为两个 upper 不完全 gamma 函数的差 (reference)。通过将 两个 参数传递给 gammainc 来计算上不完全 gamma 函数（即，z1 被隐含地视为正无穷大）。

def Gamma(a, z0, z1):
    return gammainc(a, z0) - gammainc(a, z1)

print(Gamma(12.5+17.5j, 0.0, -12.5)) 

打印(2.3801203496987e-7 + 5.54827238374855e-7j) 与 WolframAlpha 一致。