lgama Subroutine

subroutine lgama ( kf, x, gl )

  !*****************************************************************************80
  !
  !! LGAMA computes the gamma function or its logarithm.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    15 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) KF, the argument code.
  !    1, for gamma(x);
  !    2, for ln(gamma(x)).
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) GL, the function value.
  !
  implicit none

  real ( kind = 8 ), save, dimension ( 10 ) :: a = (/ &
       8.333333333333333D-02, &
       -2.777777777777778D-03, &
       7.936507936507937D-04, &
       -5.952380952380952D-04, &
       8.417508417508418D-04, &
       -1.917526917526918D-03, &
       6.410256410256410D-03, &
       -2.955065359477124D-02, &
       1.796443723688307D-01, &
       -1.39243221690590D+00 /)
  real ( kind = 8 ) gl
  real ( kind = 8 ) gl0
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kf
  integer ( kind = 4 ) n
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x2
  real ( kind = 8 ) xp

  x0 = x

  if ( x == 1.0D+00 .or. x == 2.0D+00 ) then
     gl = 0.0D+00
     if ( kf == 1 ) then
        gl = 1.0D+00
     end if
     return
  else if ( x <= 7.0D+00 ) then
     n = int ( 7.0D+00 - x )
     x0 = x + n
  end if

  x2 = 1.0D+00 / ( x0 * x0 )
  xp = 6.283185307179586477D+00
  gl0 = a(10)

  do k = 9, 1, -1
     gl0 = gl0 * x2 + a(k)
  end do

  gl = gl0 / x0 + 0.5D+00 * log ( xp ) + ( x0 - 0.5D+00 ) * log ( x0 ) - x0

  if ( x <= 7.0D+00 ) then
     do k = 1, n
        gl = gl - log ( x0 - 1.0D+00 )
        x0 = x0 - 1.0D+00
     end do
  end if

  if ( kf == 1 ) then
     gl = exp ( gl )
  end if

  return
end subroutine lgama