c-------------------------------------------------------------------
c     function  gamaln  computes ln(gamma(x))
c
c     reference:
c        Pike, M.C. and Hill, I.D. (1966)
c        Logarithm of gamma function
c        Algorithm 291, CACM vol. 9, no. 9, p 684
c
c     input variable:
c       x = any positive real number
c     output variable:
c       gamaln = ln(gamma(x))
c
c     author: jason c. hsu
c     last revision: 05-28-85 by jch
c
c     auxiliary function required: none
c-------------------------------------------------------------------
      doubleprecision function gamaln(xx)
      implicit double precision(a-h,o-z)
c6    implicit real*8(a-h,o-z)
      x=xx
      if (x.ge.7.d0) goto 10
c       ------------------------------------------
c       if x.lt.7, make x=f*x greater than 7,
c       use stirling's approximation to get ln(gamma(x)),
c       then substract f=-ln(f) from ln(gamma(x))
c       ------------------------------------------
        f=1.0d0
        z=x-1.d0
   30   z=z+1.d0
        if (z.ge.7.d0) goto 20
          x=z
          f=f*z
        goto 30
   20   x=x+1.d0
        f=-dlog(f)
        goto 40
   10 f=0.d0
   40 y2=1.d0/(x*x)
      gamaln=f+(x-.5d0)*dlog(x)-x+.918938533204673d0+
     *       (((-.000595238095238*y2+.000793650793651d0)*y2
     *       -.002777777777778d0)*y2+.083333333333333d0)/x
c       ------------------------------------------------------------
c       except for  f, the first line is stirling's approximation
c       for ln(gamma(x)), not ln(gamma(x+1))
c       (note that (x-.5)*ln(x) = -ln(x)+(x+.5)*ln(x), and
c         ln(2*pi)/2 = .91893 ...)
c       the second and third lines form the correction term
c       ------------------------------------------------------------
      return
      end