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
  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