Gamma function

Implement one algorithm (or more) to compute the Gamma (Γ) function (in the real field only).If your language has the function as builtin or you know a library which has it, compare your implementation's results with the results of the builtin/library function. The Gamma function can be defined as

$\Gamma (x)=\int _{0}^{\infty }dt\;t^{x-1}e^{-t}$

This suggests a straightforward (but inefficient) way of computing the Γ through numerical integration. Better suggested methods:

Fortran

This code shows two methods: Numerical Integration through Simpson formula, and Lanczos approximation. The results of testing are printed altogether with the values given by the function gamma; this function is defined in the Fortran 2008 standard, and supported by GNU Fortran (and other vendors) as extension; if not present in your compiler, you can remove the last part of the print in order to get it compiled with any Fortran 95 compliant compiler.

Works with: Fortran version 2008

Works with: Fortran version 95 with extensions

<lang fortran>program ComputeGammaInt

 implicit none

 integer :: i

 write(*, "(3A15)") "Simpson", "Lanczos", "Builtin"
 do i=1, 10
    write(*, "(3F15.8)") my_gamma(i/3.0), lacz_gamma(i/3.0), gamma(i/3.0)
 end do

contains

 pure function intfuncgamma(x, y) result(z)
   real :: z
   real, intent(in) :: x, y
   
   z = x**(y-1.0) * exp(-x)
 end function intfuncgamma

 function my_gamma(a) result(g)
   real :: g
   real, intent(in) :: a

   real, parameter :: small = 1.0e-4
   integer, parameter :: points = 100000

   real :: infty, dx, p, sp(2, points), x
   integer :: i
   logical :: correction

   x = a

   correction = .false.
   ! value with x<1 gives \infty, so we use
   ! \Gamma(x+1) = x\Gamma(x)
   ! to avoid the problem
   if ( x < 1.0 ) then
      correction = .true.
      x = x + 1
   end if

   ! find a "reasonable" infinity...
   ! we compute this integral indeed
   ! \int_0^M dt t^{x-1} e^{-t}
   ! where M is such that M^{x-1} e^{-M} ≤ \epsilon
   infty = 1.0e4
   do while ( intfuncgamma(infty, x) > small )
      infty = infty * 10.0
   end do

   ! using simpson
   dx = infty/real(points)
   sp = 0.0
   forall(i=1:points/2-1) sp(1, 2*i) = intfuncgamma(2.0*(i)*dx, x)
   forall(i=1:points/2) sp(2, 2*i - 1) = intfuncgamma((2.0*(i)-1.0)*dx, x)
   g = (intfuncgamma(0.0, x) + 2.0*sum(sp(1,:)) + 4.0*sum(sp(2,:)) + &
        intfuncgamma(infty, x))*dx/3.0

   if ( correction ) g = g/a

 end function my_gamma

 function lacz_gamma(a) result(g)
   real, intent(in) :: a
   real :: g

   real, parameter :: pi = 3.14159265358979324
   integer, parameter :: cg = 7

   ! these precomputed values are taken by the sample code in Wikipedia,
   ! and the sample itself takes them from the GNU Scientific Library
   real, dimension(0:8), parameter :: p = &
        (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, &
        771.32342877765313, -176.61502916214059, 12.507343278686905, &
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)

   real :: t, w, x
   integer :: i

   x = a

   if ( x < 0.5 ) then
      g = pi / ( sin(pi*x) * gamma(1.0-x) )
   else
      x = x - 1.0
      t = p(0)
      do i=1, cg+2
         t = t + p(i)/(x+real(i))
      end do
      w = x + real(cg) + 0.5
      g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t
   end if
 end function lacz_gamma

end program ComputeGammaInt</lang>

Output:

        Simpson        Lanczos        Builtin
     2.65968132     2.67893839     2.67893839
     1.35269761     1.35411859     1.35411787
     1.00000060     1.00000024     1.00000000
     0.88656044     0.89297968     0.89297950
     0.90179849     0.90274525     0.90274531
     0.99999803     1.00000036     1.00000000
     1.19070935     1.19063985     1.19063926
     1.50460517     1.50457609     1.50457561
     2.00000286     2.00000072     2.00000000
     2.77815390     2.77816010     2.77815843