d9lgic.f
853 Bytes
double precision function d9lgic (a, x, alx)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the log complementary incomplete gamma function for large x
c and for a .le. x.
c
double precision a, x, alx, eps, fk, p, r, s, t, xma, xpa,
1 d1mach, dlog
external d1mach
data eps / 0.d0 /
c
if (eps.eq.0.d0) eps = 0.5d0*d1mach(3)
c
xpa = x + 1.0d0 - a
xma = x - 1.d0 - a
c
r = 0.d0
p = 1.d0
s = p
do 10 k=1,300
fk = k
t = fk*(a-fk)*(1.d0+r)
r = -t/((xma+2.d0*fk)*(xpa+2.d0*fk)+t)
p = r*p
s = s + p
if (dabs(p).lt.eps*s) go to 20
10 continue
call seteru (57hd9lgic no convergence in 300 terms of continued f
1raction, 57, 1, 2)
c
20 d9lgic = a*alx - x + dlog(s/xpa)
c
return
end