Approximating the
Incomplete Gamma Function

Daniel D. Warner
July 25, 1999

The Incomplete Gamma Function and its complement are defined by

[Maple Math] , and [Maple Math] .

From the definition of the Gamma Function we have the identity

[Maple Math] .

The Incomplete Gamma Function is of particular interest in approximating the Cummulative Poisson Probability Function and the Cummualtive
[Maple Math] Probability Function.

[Maple Math] and [Maple Math] .

The Incomplete Gamma Function is also useful in approximating the Error Function and the Exponential Integrals

[Maple Math] and [Maple Math] .

Except for the Exponential Integrals, these applications are restricted to
[Maple Math] and [Maple Math] . Although the Incomplete Gamma Function can be defined for nonpositive [Maple Math] and negative [Maple Math] , we will not consider that situation in this study.

There is a rich collection of material on approximations for the Incomplete Gamma Function, starting with the material in Abramowitz and Stegun and ending with the sophisticated algorithm published by Walter Gautschi in 1979. Both Gautschi and the authors of the Numerical Recipes book employ a series expansion and a continued fraction expansion. In the Numerical Recipes approach the
[Maple Math] is calculated using a series when [Maple Math] , and [Maple Math] is calculated using a continued fraction when [Maple Math] . Gautschi's algorithm takes a similar approach but there are more regions and the shape of the regions is more subtle. In part, the added complexity is caused by the fact that he is developing an approximation that willl work for nonpositive values of [Maple Math] . However, some of the complexity comes from avoiding unnecessary loss of precision through cancillation. That will be a concern in the approximation developed here.

> IG := (a,x) -> (1/GAMMA(a)) * int(exp(-t)*t^(a-1),t = x .. infinity);
Ig := (a,x) -> (1/GAMMA(a)) * int(exp(-t)*t^(a-1),t = 0..x);

[Maple Math]

[Maple Math]

> plot({IG(0.5,x),IG(1.0,x),IG(3.0,x),IG(10.0,x)},x=0..15);

[Maple Plot]

> plot({Ig(0.5,x),Ig(1.0,x),Ig(3.0,x),Ig(10.0,x)},x=0..15);

[Maple Plot]

The following continued fraction, due to Legendre, will converge for positive [Maple Math] and for arbitrary real [Maple Math] .

[Maple Math] [Maple Math]

The approach in the Numerical Recipes is to compute this continued fraction using a general continued fraction algorithm. However, Gautschi introduces two transformations which reduce the calculation to the calculation of a series, whose convergence is easier to verify.
This yields the approximation
[Maple Math] [Maple Math] ,
where
[Maple Math] , and [Maple Math] ,
where
[Maple Math] , and [Maple Math] ,
and where
[Maple Math] = [Maple Math] .

The following calculations serve to verify the transformations.

> t[0] := 1:
s[0] := 1:
z[0] := t[0]/(x+1-a);
a[1] := 1*(a-1)/((x+2*1-a)^2-1):
s[1] := 1/(1+a[1]*s[0]):
t[1] := t[0]*(s[1]-1):
z[1] := normal((t[0]+t[1])/(x+1-a));
a[2] := 2*(a-2)/((x+2*2-a)^2-1):
s[2] := 1/(1+a[2]*s[1]):
t[2] := t[1]*(s[2]-1):
z[2] := normal((t[0] + t[1] + t[2])/(x+1-a));

[Maple Math]

[Maple Math]

[Maple Math]

> z[0] := normal(1/(x+(1-a)/(1)));
z[1] := normal(1/(x+(1-a)/(1+1/(x+(2-a)/(1)))));
z[2] := normal(1/(x+(1-a)/(1+1/(x+(2-a)/(1+2/(x+(3-a)/(1)))))));

[Maple Math]

[Maple Math]

[Maple Math]

And the following calculations verify the actual algorithm.

> p0 := 0:
q0 := (x-1-a)*(x+1-a):
r0 := 4*(x+1-a):
s0 := -a+1:
u0 := 0:
t0 := 1:
z0 := t0/(x+1-a);
p1 := p0 + s0:
q1 := simplify(q0 + r0):
r1 := r0 + 8:
s1 := s0 + 2:
v := p1*(1+u0):
u1 := simplify(v/(q1-v)):
t1 := u1*t0:
z1 := normal((t0+t1)/(x+1-a));
p2 := p1 + s1:
q2 := q1 + r1:
r2 := r1 + 8:
s2 := s1 + 2:
v := p2*(1+u1):
u2 := v/(q2-v):
t2 := u2*t1:
z2 := normal((t0+t1+t2)/(x+1-a));

[Maple Math]

[Maple Math]

[Maple Math]

Now consider [Maple Math] , [Maple Math] and [Maple Math] , then [Maple Math] will be positive and [Maple Math] . This means that the series will start off as an alternating series whose terms have monotone decreasing magnitudes. On the other hand, once [Maple Math] is larger than [Maple Math] then [Maple Math] can grow, but since [Maple Math] , we can show that [Maple Math] remains between 1 and 2.

> IGCF := proc( aa,xx )
local x, a, tol, p, q, r, s, u, v, t, w, z, k;
x := evalf( xx );
a := evalf( aa );
tol := evalf(10^(-15));
p := 0;
q := (x-1.0-a)*(x+1.0-a);
r := 4.0*(x+1.0-a);
s := 1.0-a;
u := 0;
t := 1;
z := t;
for k from 1 to 100 do
w := t;
p := p + s;
q := q + r;
r := r + 8.0;
s := s + 2.0;
v := p*(1.0+u);
u := v/(q-v);
t := u*t;
z := z+t;
if ((0<w) and (0<u) and (w < tol*(1-u))) then
break
fi;
od;
RETURN( evalf( z/(x+1-a) ) );
end:

> IgS := proc( aa, xx )
local a, x, tol, t, z, s, r, k;
x := eval( xx );
a := eval( aa );
tol := evalf(10^(-15));
t := 1;
z := t;
for k from 1 to 100 do
s := t;
r := x/(a+k);
t := t*r;
z := z+t;
if (s < tol*(1-r)) then
break
fi;
od;
RETURN( evalf( z ) );
end:

> IG1 := (a,x) -> `if`(x<(a+1), evalf(1-exp(-x+a*log(x)-lnGAMMA(a+1))*IgS(a,x)), evalf(exp(-x+a*log(x)-lnGAMMA(a))*IGCF(a,x))):
Ig1 := (a,x) -> `if`(x<(a+1), evalf(exp(-x+a*log(x)-lnGAMMA(a+1))*IgS(a,x)), evalf(1-exp(-x+a*log(x)-lnGAMMA(a))*IGCF(a,x))):

> plot({IG1(0.5,x),IG1(1.0,x),IG1(3.0,x),IG1(10.0,x)},x=0..15);

[Maple Plot]

> plot({Ig1(0.5,x),Ig1(1.0,x),Ig1(3.0,x),Ig1(10.0,x)},x=0..15);

[Maple Plot]

> Digits := 20;

[Maple Math]

> plot3d((IG - IG1), 0..10, 0..15, title="Plot of the Absolute Errors",axes=BOX);

[Maple Plot]

> plot3d((IG - IG1)/IG, 0..10, 0..15, title="Plot of the Relative Errors",axes=BOX);

[Maple Plot]

> readlib(C);
C(IGCF,ansi);

[Maple Math]

#include <math.h>

double IGCF(double aa,double xx)

{

  double x;

  double a;

  double tol;

  double p;

  double q;

  double r;

  double s;

  double u;

  double v;

  double t;

  int w;

  double z;

  int k;

  {

    x = 0.1E1*xx;

    a = 0.1E1*aa;

    tol = 0.1E-14;

    p = 0.0;

    q = (x-0.1E1-a)*(x+0.1E1-a);

    r = 0.4E1*x+0.4E1-0.4E1*a;

    s = 0.1E1-a;

    u = 0.0;

    t = 1.0;

    z = t;

    for(k = 1;k <= 100.0;k++)

    {

      w = t;

      p += s;

      q += r;

      r += 0.8E1;

      s += 0.2E1;

      v = p*(0.1E1+u);

      u = v/(q-v);

      t = u*t;

      z += t;

      if( 0.0 < w && 0.0 < u && w < tol*(1.0-u) )

        break;

    }

    return(0.1E1*z/(x+1.0-a));

  }

}

> C(IgS,ansi);

#include <math.h>

double IgS(double aa,double xx)

{

  double a;

  double x;

  double tol;

  double t;

  double z;

  int s;

  double r;

  int k;

  {

    x = (* xx);

    a = (* aa);

    tol = 0.1E-14;

    t = 1.0;

    z = t;

    for(k = 1;k <= 100.0;k++)

    {

      s = t;

      r = x/(a+k);

      t = t*r;

      z += t;

      if( s < tol*(1.0-r) )

        break;

    }

    return(0.1E1*z);

  }

}

>

References:

Abramowitz, M. and Stegun, I.A.,
Handbook of Mathematical Functions ,
Applied Mathematics Series, Vol. 55, National Bureau of Standards, 1964.
reprinted by Dover Publications, 1968.

Walter Gautschi,
"A Computational Procedure for Incomplete Gamma Functions",
ACM Transactions on Mathematical Software ,
Vol. 5, No. 4, December 1979, pp. 466-481.

Luke, Y. L.,
The Special Functions and Their Approximations ,
Volumes I and II, Academic Press, 1969.

William H. Press, Brian P. Flannery, Saul A. Teukolsky, and Wiliam T. Vetterling,
Numericl Recipes, The Art of Scientific Computing ,
Cambridge University Press. 1986.