Blame view
fvn_fnlib/r9lgit.f
1.02 KB
38581db0c git-svn-id: https... |
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 |