complex(kind(1.d0)) function z8lgmc (zin)

        implicit none
  c may 1978 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c compute the log gamma correction term for almost all z so that
  c   clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c8lgmc(z)
  c this routine is valid even when real(z) is negative or when cabs(z)
  c is small.
  c when real(z) is negative, c8lgmc merely returns a correction which
  c may be wrong by a multiple of 2*pi*i.
  c

        complex(kind(1.d0)) zin, z, corr,  z9lgmc, z0lgmc, zexp, zlnrel
        real(kind(1.d0)) d1mach,pi,bound,sqeps,eps,x,y,absz,test

        integer nterm,n,i,irold,ir

  c      complex(kind(1.d0)) tmp_arg

        external  z0lgmc, z9lgmc, zlnrel,
       1  d1mach
        data pi / 3.1415926535 8979324d0 /
        data bound, sqeps, eps / 3*0.0 /
  c
        if (bound.ne.0.0) go to 10
        nterm = -0.30*log(d1mach(3))
        bound = 0.1170*dble(nterm)*
       1  (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0))
        eps = 10.0*d1mach(4)
        sqeps = sqrt (d1mach(4))
  c
   10   z = zin
        x = real(z)
        y = aimag(z)
        if (x.lt.0.0) z = -z
  c
        corr = (0.0, 0.0)
        if (x.gt.bound .or. abs(y).ge.bound) go to 30
        absz = abs(z)
        if (absz.ge.bound) go to 30
  c
        if (x.eq.0.0 .and. y.eq.0.0) call seteru (
       1  57hz8lgmc  log gamma correction term is infinite for z = 0.0,
       2  57, 3, 2)
  c
        n = sqrt (bound**2-y**2) - abs(x) + 1.0
        do 20 i=1,n
          corr = corr + z0lgmc(z)
          z = z + 1.0
   20   continue
  c
   30   z8lgmc = z9lgmc(z) + corr
        if (x.gt.0.0) return
  c
        call entsrc (irold, 1)
        test = 1.0
        if (x.lt.(-0.5)) test = abs ((z-aint(x-0.5))/z)
        if (test.lt.eps) call seteru (
       1  31hz8lgmc  z is a negative integer, 31, 2, 2)
  c
        z = zin
        if (y.lt.0.0) z = conjg(z)

  c      tmp_arg = -exp(cmplx(0.,2.*pi,kind(1.d0))*z)
        corr = -cmplx(0.0,pi,kind(1.d0)) + 
       1 zlnrel(-exp(cmplx(0.,2.*pi,kind(1.d0))*z))

        if (y.lt.0.0) corr = conjg(corr)
  c
        z8lgmc = corr - z8lgmc
  c
        call erroff
        call entsrc (ir, irold)
        if (test.lt.sqeps) call seteru ( 63hc8lgmc  answer lt half precisi
       1on, z too near a negative integer, 63, 1, 1)
  c
        return
        end