r9lgit.f
1.02 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
function r9lgit (a, x, algap1)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the log of tricomi-s incomplete gamma function with perron-s
c continued fraction for large x and for a .ge. x.
c
external r1mach
data eps, sqeps / 2*0.0 /
c
if (eps.eq.0.0) eps = 0.5*r1mach(3)
if (sqeps.eq.0.0) sqeps = sqrt(r1mach(4))
c
if (x.le.0.0 .or. a.lt.x) call seteru (
1 35hr9lgit x should be gt 0.0 and le a, 35, 2, 2)
c
ax = a + x
a1x = ax + 1.0
r = 0.0
p = 1.0
s = p
do 20 k=1,200
fk = k
t = (a+fk)*x*(1.0+r)
r = t/((ax+fk)*(a1x+fk)-t)
p = r*p
s = s + p
if (abs(p).lt.eps*s) go to 30
20 continue
call seteru (57hr9lgit no convergence in 200 terms of continued f
1raction, 57, 3, 2)
c
30 hstar = 1.0 - x*s/a1x
if (hstar.lt.sqeps) call seteru (
1 39hr9lgit result less than half precision, 39, 1, 1)
c
r9lgit = -x - algap1 - alog(hstar)
c
return
end