Next: Definition of Up: The Goodness of Fit Previous: The Euler Gamma Function:

Definition of

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

 (2.21)

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

 (2.22)

Now, if we were to define a new function , then we would have:
 = (2.23) (2.24)

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

 (2.25)

What this means is that 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 and make use of the recurrence relation in order to outline precisely the computational procedure

 = = = = (2.26)

As you see, 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.
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.

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 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 , which saves us having to call x**a, which is more involved computationally.

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