function alngam (x)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
c
      external  gamma, r1mach, r9lgmc
      data sq2pil / 0.9189385332 0467274e0/
c sq2pil = alog(sqrt(2.*pi)),  sqpi2l = alog (sqrt(pi/2.))
      data sqpi2l / 0.2257913526 4472743e0/
      data pi     / 3.1415926535 8979324e0/
c
      data xmax, dxrel / 0., 0. /
c
      if (xmax.ne.0.) go to 10
      xmax = r1mach(2)/alog(r1mach(2))
      dxrel = sqrt (r1mach(4))
c
 10   y = abs(x)
      if (y.gt.10.0) go to 20
c
c alog (abs (gamma(x))) for  abs(x) .le. 10.0
c
      alngam = alog (abs (gamma(x)))
      return
c
c alog (abs (gamma(x))) for abs(x) .gt. 10.0
c
 20   if (y.gt.xmax) call seteru (
     1  38halngam  abs(x) so big alngam overflows, 38, 2, 2)
c
      if (x.gt.0.) alngam = sq2pil + (x-0.5)*alog(x) - x + r9lgmc(y)
      if (x.gt.0.) return
c
      sinpiy = abs (sin(pi*y))
      if (sinpiy.eq.0.) call seteru (31halngam  x is a negative integer,
     1  31, 3, 2)
c
      alngam = sqpi2l + (x-0.5)*alog(y) - x - alog(sinpiy) - r9lgmc(y)
c
      if (abs((x-aint(x-0.5))*alngam/x).lt.dxrel) call seteru (
     1'alngam  answer lt half precision because x too near negative ',
     2   68, 1, 1)
      return
c
      end