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