!*********************************************
!* Program to demonstrate the Gamma Function *
!* ----------------------------------------- *
!* Reference:                                * 
!* "Numerical Recipes, by W.H. Press, B.P.   *
!*  Flannery, S.A. Teukolsky and T. Vetter-  *
!*  ling, Cambridge University Press, 1986"  *
!*  [BIBLI 08].                              *
!* ----------------------------------------- *
!* SAMPLE RUN:                               *
!*                                           *
!*       X        Gamma(X)                   *   
!*  -------------------------                *
!*    0.5000    1.772453851                  *
!*    1.0000    1.000000000                  *
!*    1.5000    0.886226925                  *
!*    2.0000    1.000000000                  *
!*    2.5000    1.329340388                  *
!*    3.0000    2.000000000                  *
!*    3.5000    3.323350970                  *
!*    4.0000    5.999999999                  *
!*    4.5000   11.631728395                  *
!*    5.0000   23.999999996                  *
!*                                           *
!*********************************************               
PROGRAM Test_Gamma

  real*8 x, y
  integer i  

  print *,' '
  print *,'      X        Gamma(X)   '
  print *,' -------------------------'
  x=0.d0
  do i=1, 10
    x=x+0.5d0
    y = Gamma(x)
    write(*,10)  x, y
  end do
  print *,' '

  stop

10 format(F10.4,F15.9)

end

!*******************************************
!*           FUNCTION  GAMMA(X)            *
!* --------------------------------------- *
!* Returns the value of Gamma(x) in double *
!* precision as EXP(LN(GAMMA(X))) for X>0. *
!*******************************************
real*8 Function Gamma(xx)
  real*8 xx,cof(6),stp,half,one,fpf,x,tmp,ser
  DATA cof,stp /76.18009173d0,-86.50532033d0,24.01409822d0,  &
       -1.231739516d0,0.120858003d-2,-0.536382d-5,2.50662827465d0/
  DATA half,one,fpf /0.5d0,1.0d0,5.5d0/
  x=xx-one
  tmp=x+fpf
  tmp=(x+half)*DLOG(tmp)-tmp
  ser=one
  do j=1,6
    x=x+one
    ser=ser+cof(j)/x
  end do
  Gamma = DEXP(tmp+DLOG(stp*ser))
  return
end

!end of file gamma.f90