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.

closed as off-topic by MarcoB, m_goldberg, Alex Trounev, LCarvalho, José Antonio Díaz NavasMar 21 at 19:50

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – MarcoB, m_goldberg, Alex Trounev, LCarvalho, José Antonio Díaz Navas
If this question can be reworded to fit the rules in the help center, please edit the question.

• 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 – J. M. is away Mar 13 at 0:09
• @J.M.isslightlypensive thanks -- was being slow! – ben18785 Mar 13 at 0:27
• Maybe one of these formulas will help: functions.wolfram.com/GammaBetaErf/Gamma3/26/01 – Roman Mar 13 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. – J. M. is away Mar 13 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/… – Roman Mar 13 at 2:31

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]. – J. M. is away Mar 13 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). – J. M. is away Mar 13 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. – Roman Mar 13 at 13:39
• Thanks @J.M.isslightlypensive for the formatting help. – Roman Mar 13 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.