c234567
      Function gammln(xx)
      double precision gammln,xx
! Returns the value ln[Γ(xx)] for xx>0
      integer j
      double precision ser,stp,tmp,x,y,cof(6)
!integer arithmetic will be done in double precision, a nicety that you can
!omit if five-figure accuracy is good enough.
      SAVE cof,stp
      DATA cof,stp/76.18009172947146,-86.50532032941677,24.01409824083091,
     $       -1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5,
     $        2.5066282746310005/
      x=xx
      y=x
      tmp=x+5.5
      tmp=(x+0.5)*log(tmp)-tmp
      ser=1.000000000190015
      do   j=1,6
        y=y+1.
        ser=ser+cof(j)/y
      enddo
      gammln=tmp+log(stp*ser/x)
      !print *, "gammln=",gammln
      return
      END