real function gammln(xx)
  implicit none
  real :: x,y,tmp,ser
  real :: cof(6) = (/76.18009172947146, -86.50532032941677, 24.01409824083091,&
                     -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5/)
  integer :: j
!
  y = xx
  x = xx
  tmp = x+5.5
  tmp = tmp-(x+0.5)*log(tmp)
  ser = 1.000000000190015
  do j=1,5
    y = y+1
    ser = ser+cof(j)/y
  enddo
  gammln = -tmp+log(2.5066282746310005*ser/x)
end function gammln