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