# Numerically stable replacement for generalised incomplete gamma function [closed]

I am looking to replace the generalised incomplete gamma function (which appears in a solution to a problem I've posted about here) with a numerically stable formula involving other functions. This is because I want to move over to other languages which do not have this function implemented.

To illustrate my point, the generalised incomplete gamma function when expanded is,

Gamma[a, b, c] = Gamma[a, b] - Gamma[a, c]


So naively you might think that the LHS and RHS should return the answer. The RHS however is numerically unclever and so has issues when $$a$$ is large and $$b$$ and $$c$$ are small. For example,

Gamma[100, 0.1, 0.01]


returns,

-9.05734*10^-103

whereas,

Gamma[100, 0.1] - Gamma[100, 0.01]


returns,

0.

My question is what does the generalised incomplete gamma function actually do in the background? And so, how can I replicate this behaviour using simpler functions?

Note 1: I do not want answers explicitly involving calls to NIntegrate or Integrate since these are inefficient to implement. Rather I would like a solution involving other functions.

Note 2: I know the definition of the incomplete Gamma in Mathematica is the upper incomplete gamma. If there were a way to write a numerically stable form using the lower incomplete gamma, that may avoid the issue since the Gamma[a] terms would cancel out.

• You are aware that the lower incomplete gamma function is expressible in terms of three-argument Gamma[]? Gamma[p, a, b] == Gamma[p, 0, b] - Gamma[p, 0, a] // FunctionExpand Mar 13 '19 at 0:09
• @J.M.isslightlypensive thanks -- was being slow! Mar 13 '19 at 0:27
• Maybe one of these formulas will help: functions.wolfram.com/GammaBetaErf/Gamma3/26/01 Mar 13 '19 at 1:05
• The expressions in that Wolfram functions link are all hypergeometric expressions, so they may not be that useful if you cannot evaluate a hypergeometric function. Mar 13 '19 at 1:37
• @J.M.isslightlypensive the GSL can do hypergeometric functions, and can be added to any other language: gnu.org/software/gsl/doc/html/… Mar 13 '19 at 2:31

## 2 Answers

As J. M. comments, Gamma[a, z1, z2] == Gamma[a, 0, z2] - Gamma[a, 0, z1]. The lower incomplete gamma function can be expressed as a confluent hypergeometric function

z^a/a Hypergeometric1F1[a, a + 1, -z]


Gamma[a, 0, z]

This form can be implemented in any language that supports the $$_1F_1$$ confluent hypergeometric function. For languages that don't, this function can be added through the GNU Scientific Library's gsl_sf_hyperg_1F1_int or gsl_sf_hyperg_1F1 function, or through its series definition that converges very quickly: (this is equal to @MichaelE2's solution)

Hypergeometric1F1[a, a + 1, -z] == Sum[(-z)^k/((1 + k/a) k!), {k, 0, ∞}]


True

The GNU Scientific Library can also calculate the regularized lower incomplete gamma function directly via gsl_sf_gamma_inc_P; but this method is underflow-problematic for very small values of $$z$$ when the parameter $$a$$ is large.

Here's a C code that calculates the lower incomplete gamma function and the generalized incomplete gamma function, for real and integer first argument (different algorithms). This can be incorporated into any other programming language.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_hyperg.h>

// lower incomplete gamma function
double ligamma(const double a, const double z) {
return pow(z,a)*gsl_sf_hyperg_1F1(a,a+1,-z)/a; }
double ligamma_int(const int a, const double z) {
return pow(z,a)*gsl_sf_hyperg_1F1_int(a,a+1,-z)/a; }

// generalized incomplete gamma function
double gigamma(const double a, const double z1, const double z2) {
return ligamma(a,z2)-ligamma(a,z1); }
double gigamma_int(const int a, const double z1, const double z2) {
return ligamma_int(a,z2)-ligamma_int(a,z1); }

int main() {
printf("%.16g\n", gigamma_int(100, 0.1, 0.01));
printf("%.16g\n", gigamma(100.0, 0.1, 0.01));
}


-9.057341758336185e-103

-9.057341758336185e-103

• The scaled incomplete gamma is called GammaRegularized[] in Mathematica. It also supports a three-argument form, so scaled lower incomplete gamma would be GammaRegularized[a, 0, z]. Mar 13 '19 at 13:20
• Yes, the reason why one would have both $P$ and $Q$ around in a computing environment is precisely to avoid having to do $1-P$ or $1-Q$ (similar remarks hold for the $\operatorname{erf}$/$\operatorname{erfc}$ pair). Mar 13 '19 at 13:36
• Unfortunately the regularization of $P$ and $Q$ leads to problems at large $a$ because $\Gamma(a)$ becomes enormous. It would be great to have unregularized functions available directly in the GSL; but as I wrote they can be calculated from confluent hypergeometrics. I digress from Mathematica. Mar 13 '19 at 13:39
• Thanks @J.M.isslightlypensive for the formatting help. Mar 13 '19 at 14:11

For a large and b and c small, the series expansion works pretty well:

Block[{a = 100., b = 0.1, c = 0.01},
{Sum[((-1)^n (c^(a + n) - b^(a + n)))/((a + n) n!), {n, 0, 9}],
N[Gamma @@ Rationalize[{a, b, c}], 16]}
]
Subtract @@ %/First[%]
(*
{-9.05734*10^-103, -9.057341758336136*10^-103}
5.69173*10^-15
*)


Just go out until the term is less than the sum times \$MachineEpsilon. If you need help with that, I can do it, but first I was wondering whether the relative error is ok. It's about 25 times worse than the theoretically best possible approximation.