z8lgmc.f 2.17 KB
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