complex(kind(1.d0)) function zlngam (zin) implicit none c august 1980 edition. w. fullerton c3, los alamos scientific lab. c eventually clngam should make use of c8lgmc for all z except for c z in the vicinity of 1 and 2. complex(kind(1.d0)) zin, z, corr, zlnrel, z9lgmc real(kind(1.d0)) d1mach,pi,sq2pil,bound,dxrel,rmax,cabsz real(kind(1.d0)) x,y,argsum,zarg integer irold,ir,i,n external z9lgmc, zarg, zlnrel,d1mach data pi / 3.1415926535 8979324d0 / data sq2pil / 0.9189385332 0467274d0 / c data bound, dxrel, rmax / 3*0.0d0 / c if (bound.ne.0.) go to 10 n = -0.30*log(d1mach(3)) c bound = n*(0.1*eps)**(-1/(2*n-1))/(pi*exp(1)) bound = 0.1171*dble(n)*(0.1*d1mach(3))**(-1./(2.*dble(n)-1.)) dxrel = sqrt (d1mach(4)) rmax = d1mach(2)/log(d1mach(2)) c 10 z = zin x = real(zin) y = aimag(zin) c corr = (0.0, 0.0) cabsz = abs(z) if (cabsz.gt.rmax) call seteru ( 1 44hzlngam abs(z) so large result may overflow, 44, 2, 2) if (x.ge.0.0 .and. cabsz.gt.bound) go to 50 if (x.lt.0.0 .and. abs(y).gt.bound) go to 50 c if (cabsz.lt.bound) go to 20 c c use the reflection formula for real(z) negative, cabs(z) large, and c abs(aimag(y)) small. c call entsrc (irold, 1) if (y.gt.0.0) z = conjg (z) corr = exp (-cmplx(0.0,2.0*pi,kind(1.d0))*z) if (real(corr).eq.1.0 .and. aimag(corr).eq.0.0) call seteru ( 1 31hzlngam z is a negative integer, 31, 3, 2) c zlngam = sq2pil + 1.0 - cmplx(0.0,pi,kind(1.d0))*(z-0.5) - 1 zlnrel(-corr) + (z-0.5)*log(1.0-z) - z - z9lgmc(1.0-z) if (y.gt.0.0) zlngam = conjg (zlngam) c call erroff call entsrc (ir, irold) if (abs(y).gt.dxrel) return c X in sequence field on next line is to preserve trailing blank if (0.5* abs((1.-corr)*zlngam/z).lt.dxrel) call seteru (68hzlngam X 1 answer lt half precision because z too near negative integer, 68, 2 1, 1) return c c use the recursion relation for cabs(z) small. c 20 if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30 if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru ( 68hzlngam ans 1wer lt half precision because z too near negative integer, 68,1,1) c 30 n = sqrt (bound**2 - y**2) - x + 1.0 argsum = 0.0 corr = (1.0, 0.0) do 40 i=1,n argsum = argsum + zarg(z) corr = z*corr z = 1.0 + z 40 continue c if (real(corr).eq.0.0 .and. aimag(corr).eq.0.0) call seteru ( 1 31hzlngam z is a negative integer, 31, 3, 2) corr = -cmplx (log(abs(corr)), argsum, kind(1.d0)) c c use stirling-s approximation for large z. c 50 zlngam = sq2pil + (z-0.5)*log(z) - z + corr + z9lgmc(z) return c end