Definition of $\gamma (a, x)$

We'll write a subprogram for calculating $\gamma (a, x)$, because its mathematical definition is less scary. This definition will look simpler still when we recall the following recurrence relation:

$$\Gamma(a + 1) = a \Gamma(a) \Rightarrow \frac{\Gamma(a + 1)}{\Gamma(a)} = a$$ (2.21)

Recall that the factorial function, n!, has a very similar recurrence:

$$(n+1)! = n! (n+1) \Rightarrow \frac{(n+1)!}{n!} = n+1$$ (2.22)

Now, if we were to define a new function $\Gamma'(n) = \Gamma(n+1)$, then we would have:

$$\Gamma'(n+1) \Gamma(n+2) = (n+1)\Gamma(n+1) = (n+1)\Gamma'(n)$$ = (2.23)

$$\Rightarrow \frac{\Gamma'(n+1)}{\Gamma'(n)} = n+1$$ (2.24)

This relationship is the same as the recurrence relationship for a Factorial, and indeed we have that

$$n! = \Gamma(n+1)$$ (2.25)

What this means is that $\Gamma$ is an extremely rapidly rising function, which is why UNIX defines only its logarithm in libm.a. It is safer that way.

Let us again have a look at the series expansion for $\gamma (a, x)$and make use of the recurrence relation in order to outline precisely the computational procedure

$$\gamma(a, x) e^{-x}x^a\sum_{n=0}^\infty \frac{\Gamma(a)x^n}{\Gamma(a + 1 + n)}$$ =

$$e^{-x}x^a \left( \frac{\Gamma(a)x^0}{\Gamma(a + 1 + 0)} + \frac{\... ...}{\Gamma(a + 1 + 2)} + \frac{\Gamma(a) x^3}{\Gamma(a + 1 + 3)} + \cdots \right)$$ =

$$e^{-x}x^a \left( \frac{\Gamma(a)}{a \Gamma(a)} + \frac{\Gamma(a) ... ...)} + \frac{\Gamma(a) x^3}{(a + 3) (a + 2) (a + 1) a \Gamma(a)} + \cdots \right)$$ =

$$e^{-x}x^a \left( \frac{1}{a} + \frac{x}{a (a + 1)} + \frac{x x}{a (a + 1) (a + 2)} + \frac{x x x}{a (a + 1) (a + 2) (a + 3)} + \cdots \right)$$ = (2.26)

As you see, $\Gamma (a)$ falls out altogether. And the algorithm, neglecting the e^-xx^a factor for the time being, becomes as follows:

1.

evaluate 1/a first and save in on, say, g, and at the same time on g_sum

2.

multiply g by x/(a+1), save the result back on g, and add it to whatever is in g_sum

3.

multiply g by x/(a+2), save the result back on g, and add it to g_sum

4.

multiply g by x/(a+3), save the result back on g, and add it to g_sum

5.

$\ldots$ and so on ad infinitum

In Fortran this algorithm would look as follows:

g = 1/a
g_sum = g
do i = 1, infinity
   g = g * x / (a + i)
   g_sum = g_sum + g
end do

On many processors incrementing a number by one is faster than adding an integer to it. So we can optimise the algorithm presented above as follows:

g = 1/a; a_incr = a; g_sum = g
do i = 1, infinity
   a_incr = a_incr + 1
   g = g * x / a_incr
   g_sum = g_sum + g
end do

We must now stop this loop. We stop it when adding a yet another term doesn't visibly change the result, i.e., when g is less than the number of significant digits we want in the final result. We also restructure the do loop, because we don't make any use of its index i.

real(kind=long), parameter :: epsilon = 1.0E-9
...
g = 1/a; a_incr = a; g_sum = g
do while(abs(g) .gt. abs(g_sum) * epsilon)
   a_incr = a_incr + 1
   g = g * x / a_incr
   g_sum = g_sum + g
end do

Before exiting the subprogram all we need to do is to multiply g_sum by e^-x x^a. Figure 2.12 lists our full implementation of small_gamma.

  

Figure 2.12: Function small_gamma from module euler

There are some small additions to small_gamma. First we check for a possibility that the function may have been called with x < 0. The incomplete gamma function is not defined for x < 0, so if this is the case, we stop the execution of the program and flag an error message. For x=0 $\gamma(a, x) = 0$ too, so in this case we simply return 0. Only for x > 0 the real computation starts. In the last statement of the function instead of multiplying by e^-x x^a we multiply by $\exp(-x + a \log{x})$, which saves us having to call x**a, which is more involved computationally.