gammainc NotImplementedError

[Roban]

I've been getting some NotImplementedErrors when calculating
generalized incomplete gamma function values. Here's an example:

>>> gammainc(-1, 0.3, 0.301)
mpf('0.0081998709378296071')
>>> gammainc(-1, 0.3, 0.3001)
------------------------------------------------------------
Traceback (most recent call last):
File "<ipython console>", line 1, in <module>
File "/usr/local/lib/python2.6/dist-packages/mpmath/functions/
expintegrals.py", line 166, in gammainc
return +ctx._gamma3(z, a, b, regularized)
File "/usr/local/lib/python2.6/dist-packages/mpmath/functions/
expintegrals.py", line 244, in _gamma3
raise NotImplementedError
NotImplementedError

But it works again with

>>> gammainc(-1+1e-10, 0.3, 0.3001)
mpf('0.00082281592767685652')

The problem is in the _expintegrals.gamma3 function, which chooses not
to use the difference of two upper incomplete gamma function when the
difference is much smaller than the function values themselves.
Normally this is a wise choice, but when z in a non-positive integer
it is a pole of the lower incomplete gamma function, so the
NotImplementedError is reached.

So my question is, is this justified (i.e. do the current algorithms
really produce inaccurate answers in this situation), or should the
code actually just use the difference of upper function in the case of
a pole of the lower function? If the later, I'd be glad to submit a
patch (in whatever way is easiest).

Thanks,
Roban

[Fredrik Johansson]

Hi Roban,

Basically, _gamma3 uses two different formulas: it either computes a
difference of two upper incomplete gamma functions, or a difference of
two lower incomplete gamma functions.

Right now it just tries one formula with a fixed precision and checks
if it's good, and otherwise it tries the other with the same fixed
precision. If neither is good, it raises NotImplementedError rather
than returning an inaccurate value.

A better algorithm would be to increase the precision in a loop until
either formula gives full accuracy (and maybe raise NoConvergence when
the precision gets too large). This could be done by changing a few
lines in the existing code. You're very welcome to write a patch.

To improve efficiency further, one could switch to a series expansion
for computing gamma(z,a,b) when b = a+epsilon for sufficiently small
epsilon, say |epsilon| < max(|a|,|z|) * 2^(-prec/4). The series wrt
epsilon goes something like

gamma(z,a,a+epsilon) = exp(-a) * a^z * (epsilon / a) * (1 + epsilon *
(3z-3-3a) / (6*a) + ...)

This should also be quite straightforward, in case you are interested
in writing a patch.

Fredrik