File contents

      function dlngam(x,psi,d,iflag)
c  calculates logarithm of the gamma function, its derivative psi(x), 
c  and the derivative of psi(x) for positive argument
c  gamma(x) is obtained by exponentiation
c  derivative of gamma(x) is gamma(x)*psi(x)
c  method used is continued fraction representation for arguments 
c  greater than 12 and iteration downward for arguments less than 12
c  maximum error is approximately 1 part in 10**14
c  written by R. B. Shirts, Brigham Young University
c  last altered August 2001
c  parameter iflag must be set to 0 on first call to set up arrays
c  of constants
      implicit real*8 (a-h,o-z)
      dimension a(5), b(5), c(5)
      data zero,one,half,big,tol/0.d0,1.d0,0.5d0,1.d99,1.d-4/
      data zeta2,zeta3/1.64493406684822643647d0,1.2020569031595942854d0/
c  const is 0.5*dlog(2*pi)
      data gamma,const/5.772156649015329d-1,9.18938533204673d-1/
      if (iflag.eq.0) then
c  set continued fraction parameters for lngamma
c  see Abramowitz and Stegun, 6.1.48 (p. 258)
          a(5) = one/12.d0
          a(4) = one/30.d0
          a(3) = 53.d0/210.d0
          a(2) = 195.d0/371.d0
          a(1) = 22999.d0/22737.d0
c  set continued fraction parameters for psi
          b(5) = -one/12.d0
          b(4) = one/10.d0
          b(3) = 79.d0/210.d0
          b(2) = 1205.d0/1659.d0
          b(1) = 262445.d0/209429.d0
c  set continued factions parameters for psi derivative
          c(5) = one/6.d0
          c(4) = 0.2d0
          c(3) = 18.d0/35.d0
          c(2) = 20.d0/21.d0
          c(1) = 50.d0/33.d0
          iflag=1
      endif
      if (x.le.zero) then
c  set values for negative argument
          dlngam=big
          psi=-big
          d=big
          return
      endif
      if (x.le.tol) then
c  set values for arguments near zero
          xsq=x*x
          dlngam=-dlog(x)-gamma*x+half*xsq*zeta2-x*xsq*zeta3/3.d0
          psi=-one/x-gamma+zeta2*x-zeta3*xsq
          d=one/xsq+zeta2-2.d0*zeta3*x
          return
      endif
      n=12-idint(x)
      if (n.le.0) then
          z=x
      else
          z=x+dfloat(n)
      endif
c  form continued faction representation of asympotic expansions
      frac1=zero
      frac2=zero
      frac3=zero
      do 2 i=1,5
          frac1=a(i)/(z+frac1)
          frac2=b(i)/(z+frac2)
          frac3=c(i)/(z+frac3)
    2 continue
      dlnz=dlog(z)
      dlngam=(z-half)*dlnz-z+const+frac1
      psi=dlnz+(frac2-half)/z
      d=(one+(half+frac3)/z)/z
      if (n.le.0) then
          return
      else
c  use recurrence if x < 12
      do 1 i=1,n
          z=z-one
c  Abramowitz and Stegun 6.1.15 (p. 256)
          dlngam = dlngam-dlog(z)
c  Abramowitz and Stegun 6.3.5 (p. 258)
          psi = psi-one/z
          d=d+one/(z*z)
    1 continue
      end if
      return
      end