13.3.4.1 The Gamma Function

13.3.4.1  The Gamma Function

The gamma function was first introduced by Euler to generalize the factorial function to non-integer values. It has been studied by eminent mathematicians such as Gauss, Legendre, Liouvlle, Hermite and Weirstarss. Some famous mathematical constants such as the Euler constant g are related to the gamma function.

The gamma function is denoted by G(n) and is defined in terms of a special and useful integral. The definition is given below.

                                        (11)

The integral can be shown to converge for all values n>0. If we integrate by parts and take the limiting values, we can obtain the following recurrence relation for describing the gamma function.

        G(n+1)  =  nG(n)       (12)

By evaluating the integral 11 for n=1, we can show that G(1)=1. From Equation 12, G(n) can be determined for all n>0 when the values for the gamma function have been determined for a certain interval of unit length, e.g., 1£n<2. In particular, it can be shown that if n is a positive integer, the following holds for n=1,2,3,×××.

             G(n+1)                      =                     n!                (13)

For this reason, G(n) is sometimes called the factorial function. It also can be shown that

                                        (14)

The recurrence relation 12 is a difference equation that has Equation 11 as a solution.

The following program simply prints the value of the gamma function for a sequence of arguments.

 Program 13.5

#!/usr/bin/perl
#gamma.pl
use Math::Cephes qw (:gammas);
use strict;
for (my $i = -1; $i <= 100; $i = $i + 3.5 ){
    my $gammaVal = Math::Cephes::gamma ($i);
    print "gamma($i) = $gammaVal\n";
  }

 The output of the program is given below.

gamma domain error

gamma(-1) = nan

gamma(2.5) = 1.32934038817914

gamma(6) = 120

gamma(9.5) = 119292.461994609

gamma(13) = 479001600

gamma(16.5) = 5189998453040.12

gamma(20) = 1.21645100408832e+17

gamma(23.5) = 5.36130358754441e+21

gamma(27) = 4.03291461126606e+26

gamma(30.5) = 4.82269693349091e+31

gamma(34) = 8.68331761881189e+36

gamma(37.5) = 2.25511578410651e+42

gamma(41) = 8.15915283247898e+47

gamma(44.5) = 3.99612665510252e+53

gamma(48) = 2.58623241511168e+59

gamma(51.5) = 2.16668377073774e+65

gamma(55) = 2.30843697339241e+71

gamma(58.5) = 3.07979367469713e+77

gamma(62) = 5.07580213877225e+83

gamma(65.5) = 1.02102976358767e+90

gamma(69) = 2.48003554243683e+96

gamma(72.5) = 7.20403236376959e+102

gamma(76) = 2.48091408113954e+109

gamma(79.5) = 1.00493279593549e+116

gamma(83) = 4.75364333701284e+122

gamma(86.5) = 2.60868045964056e+129

gamma(90) = 1.65079551609085e+136

gamma(93.5) = 1.19790604809448e+143

gamma(97) = 9.9167793487095e+149

The error statements on top show that the gamma value cannot be computed by the package for negative values. nan seen on the second line stands for not a number. Although the gamma function can be defined for negative integers using a process called analytical continuation, it has not been implemented in the Math::Cephes module. There are some computations, particularly in number theory, where the logarithm of the gamma function often appears. Math::Cephes provides an implementation of such a function also although we do not discuss it here. Several variations of the gamma function are implemented as well. Of these, the so-called incomplete gamma function and its derivative are commonly used in statistics.