next up previous index
Next: Definition of Up: The Goodness of Fit Previous: The Euler Gamma Function:

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:

\begin{displaymath}\Gamma(a + 1) = a \Gamma(a) \Rightarrow \frac{\Gamma(a + 1)}{\Gamma(a)} = a
\end{displaymath} (2.21)

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

\begin{displaymath}(n+1)! = n! (n+1) \Rightarrow \frac{(n+1)!}{n!} = n+1
\end{displaymath} (2.22)

Now, if we were to define a new function $\Gamma'(n) = \Gamma(n+1)$, then we would have:
$\displaystyle \Gamma'(n+1)$ = $\displaystyle \Gamma(n+2) = (n+1)\Gamma(n+1) = (n+1)\Gamma'(n)$ (2.23)
  $\textstyle \Rightarrow$ $\displaystyle \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

\begin{displaymath}n! = \Gamma(n+1)
\end{displaymath} (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

$\displaystyle \gamma(a, x)$ = $\displaystyle e^{-x}x^a\sum_{n=0}^\infty
\frac{\Gamma(a)x^n}{\Gamma(a + 1 + n)}$  
  = $\displaystyle 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)$  
  = $\displaystyle 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)$  
  = $\displaystyle 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-xxa 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 xa. Figure 2.12 lists our full implementation of small_gamma.
  
Figure 2.12: Function small_gamma from module euler
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}real...
...og(x))
end if\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

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 xa we multiply by $\exp(-x + a \log{x})$, which saves us having to call x**a, which is more involved computationally.


next up previous index
Next: Definition of Up: The Goodness of Fit Previous: The Euler Gamma Function:
Zdzislaw Meglicki
2001-02-26