double precision function dlngam (x)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, dxrel, pi, sinpiy, sqpi2l, sq2pil,
     1  y, xmax, dint, dgamma, d9lgmc, d1mach, dlog, dsin, dsqrt
      external d1mach, d9lgmc, dgamma, dint, dlog, dsin, dsqrt
c
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
c sq2pil = alog (sqrt(2*pi)),  sqpi2l = alog(sqrt(pi/2))
      data sqpi2l / +.2257913526 4472743236 3097614947 441 d+0    /
      data pi / 3.1415926535 8979323846 2643383279 50 d0 /
c
      data xmax, dxrel / 2*0.d0 /
c
      if (xmax.ne.0.d0) go to 10
      xmax = d1mach(2)/dlog(d1mach(2))
      dxrel = dsqrt (d1mach(4))
c
 10   y = dabs (x)
      if (y.gt.10.d0) go to 20
c
c dlog (dabs (dgamma(x)) ) for dabs(x) .le. 10.0
c
      dlngam = dlog (dabs (dgamma(x)) )
      return
c
c dlog ( dabs (dgamma(x)) ) for dabs(x) .gt. 10.0
c
 20   if (y.gt.xmax) call seteru (
     1  39hdlngam  dabs(x) so big dlngam overflows, 39, 2, 2)
c
      if (x.gt.0.d0) dlngam = sq2pil + (x-0.5d0)*dlog(x) - x + d9lgmc(y)
      if (x.gt.0.d0) return
c
      sinpiy = dabs (dsin(pi*y))
      if (sinpiy.eq.0.d0) call seteru (
     1  31hdlngam  x is a negative integer, 31, 3, 2)
c
      dlngam = sqpi2l + (x-0.5d0)*dlog(y) - x - dlog(sinpiy) - d9lgmc(y)
c
      if (dabs((x-dint(x-0.5d0))*dlngam/x).lt.dxrel) call seteru (
     1  68hdlngam  answer lt half precision because x too near negative
     2integer, 68, 1, 1)
      return
c
      end