Next: Definition of
Up: The Goodness of Fit
Previous: The Euler Gamma Function:
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:
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
As you see,
falls out altogether. And the algorithm,
neglecting the e^{x}x^{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.

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.0E9
...
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
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
,
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
20010226