5

I'm trying to compute the upper incomplete gamma function defined like in this post. If I use

from scipy.special import gamma,gammainc
from numpy import linspace

a = 0
z = (2+3j)*np.linspace(0,10)
gamma(a)*(1-gammainc(a,z))

where z is a complex vector I get an error

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

Is there an alternative function to do the calculation? There doesn't seem to be a problem when I try to do this with WolframAlpha's Gamma function.

Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • The incomplete gamma function is defined through an integral with `x` (or `z`) being an integration boundary. How do you expect this integral to be performed in the complex plane? Integration in the complex plane usually depends on the integration path and not only on the boundaries. I am not aware how WolframAlpha deals with that situation but that would be the first thing to find out. Also note that for the specific case in your example (`a = 0`) there is an analytic solution for the integral and you can use the `scipy.special.expi` function to compute the result (for real value boundaries). – a_guest Jul 31 '18 at 10:38
  • @a_guest At least for arguments with positive real part it's unambiguous: there are no poles of the integrand there, and the factor `exp(-t)` provides nice decay at infinity. On the left halfplane, analytic continuation works with a branch cut along the negative real axis. –  Jul 31 '18 at 15:54

1 Answers1

4

When SciPy is not enough for tricky special functions, mpmath usually comes to the rescue. Compare

>>> mpmath.gammainc(0, 2+3j)
mpc(real='-0.024826207944199364', imag='0.020316674911044622')

and the same from Wolfram Alpha.

Being written in Python, it is slower than SciPy; it is also not vectorized. So with your data it would work like

import mpmath
result = np.array([mpmath.gammainc(w) for w in z[1:]], dtype=np.complex)

Note that I avoid passing 0 as the argument (it's a pole). The return type of mpmath.gammainc is its own mpc object type, but it can be converted back to NumPy as above.

  • Thanks, I also found a sympy function [sympy.functions.special.gamma_functions.uppergamma](http://docs.sympy.org/0.7.2/modules/functions/special.html?highlight=uppergamma#sympy.functions.special.gamma_functions.uppergamma) that addresses my problem. In terms of speed and accuracy, would you generally recommend the `mpmath` or `sympy` function(s)? I will be inputting an array containing 10^6 to 10^7 points in my real program (not the MWE in the post). – Medulla Oblongata Aug 01 '18 at 12:21
  • 1
    That SymPy function calls mpmath for numerical computations, it does not do them itself. There is no reason to go through SymPy to access an mpmath function, it adds an overhead of function calls. A large array like that will take time. Look into decreasing the precision in mpmath in case it helps, http://mpmath.org/doc/1.0.0/basics.html#setting-the-precision –  Aug 01 '18 at 12:45