lngamma Function

FUNCTION lngamma(x) RESULT(fn_val)
  ! Logarithm to base e of the gamma function.
  !
  ! Accurate to about 1.e-14.
  ! Programmer: Alan Miller
  ! Latest revision of Fortran 77 version - 28 February 1988
  REAL (dp), INTENT(IN) :: x
  REAL (dp)             :: fn_val
  !       Local variables
  REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03,  &
       a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04,  &
       temp, arg, product, lnrt2pi = 9.189385332046727D-1,       &
       pi = 3.141592653589793D0
  LOGICAL   :: reflect
  !       lngamma is not defined if x = 0 or a negative integer.
  IF (x > 0.d0) GO TO 10
  IF (x /= INT(x)) GO TO 10
  fn_val = 0.d0
  RETURN
  !       If x < 0, use the reflection formula:
  !               gamma(x) * gamma(1-x) = pi * cosec(pi.x)
10 reflect = (x < 0.d0)
  IF (reflect) THEN
     arg = 1.d0 - x
  ELSE
     arg = x
  END IF
  !       Increase the argument, if necessary, to make it > 10.
  product = 1.d0
20 IF (arg <= 10.d0) THEN
     product = product * arg
     arg = arg + 1.d0
     GO TO 20
  END IF
  !  Use a polynomial approximation to Stirling's formula.
  !  N.B. The real Stirling's formula is used here, not the simpler, but less
  !       accurate formula given by De Moivre in a letter to Stirling, which
  !       is the one usually quoted.
  arg = arg - 0.5D0
  temp = 1.d0/arg**2
  fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &
       (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)
  IF (reflect) THEN
     temp = SIN(pi * x)
     fn_val = LOG(pi/temp) - fn_val
  END IF
  RETURN
END FUNCTION lngamma