4

I am playing with the python library mpmath, in particular evaluating the incomplete gamma function. This is part of a root finding routine but its evaluation is extremely slow for some combinations of complex-valued arguments.

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

Here the evaluation of the mpmath.gammainc function gets stuck:

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

On the other hand, Mathematica returns me the result almost instantly:

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

In other cases, for different arguments mpmath and Mathematica return the same output:

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

Do you have some idea of the reason of this behaviour? Can this be considered a problem of mpmath or is this a mathematical problem of quadrature? Unfortunately scipy does not offer an implementation of gamma functions for complex arguments, so it not an option.

linello
  • 8,451
  • 18
  • 63
  • 109

1 Answers1

3

Apparently, a bug in mpmath causes it to go into an infinite loop when evaluating gammainc in some cases. Worth reporting on the mpmath tracker. But at least for the case you mentioned, a workaround is to write the three-argument incomplete gamma function as the difference of two upper incomplete gamma functions (reference). The upper incomplete gamma function is computed by passing two arguments to gammainc (i.e., z1 is implicitly taken to be positive infinity). 

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

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

prints (2.3801203496987e-7 + 5.54827238374855e-7j) in agreement with WolframAlpha.

  • Can you say why you think it's infinite recursion? Perhaps a source link, if you've identified where the code goes wrong? I would expect such infinite recursion to overflow the stack and crash relatively quickly. – user2357112 Aug 26 '18 at 00:26
  • If I try to evaluate `gammainc(12.5+17.5j, 0.0, -12.5)` and interrupt the computation after a while, the stack trace consists mostly of three lines: "expintegrals.py, line 163, expintegrals.py, line 166, expintegrals.py, line 241", repeated over and over again. Maybe an infinite loop rather than recursion. –  Aug 26 '18 at 00:33
  • After trying it out, it's definitely recursion (though I don't know whether it's theoretically infinite), and it hits the recursion limit pretty quickly. – user2357112 Aug 26 '18 at 00:47