Blame view

fvn_fnlib/d9lgit.f 1.16 KB
38581db0c   daniau   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
40
41
42
        double precision function d9lgit (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
        double precision a, x, algap1, ax, a1x, eps, fk, hstar, p, r, s,
       1  sqeps, t, d1mach, dlog, dsqrt
        external d1mach
        data eps, sqeps / 2*0.d0 /
  c
        if (eps.ne.0.d0) go to 10
        eps = 0.5d0*d1mach(3)
        sqeps = dsqrt (d1mach(4))
  c
   10   if (x.le.0.d0 .or. a.lt.x) call seteru (
       1  35hd9lgit  x should be gt 0.0 and le a, 35, 2, 2)
  c
        ax = a + x
        a1x = ax + 1.0d0
        r = 0.d0
        p = 1.d0
        s = p
        do 20 k=1,200
          fk = k
          t = (a+fk)*x*(1.d0+r)
          r = t/((ax+fk)*(a1x+fk)-t)
          p = r*p
          s = s + p
          if (dabs(p).lt.eps*s) go to 30
   20   continue
        call seteru (57hd9lgit  no convergence in 200 terms of continued f
       1raction, 57, 3, 2)
  c
   30   hstar = 1.0d0 - x*s/a1x
        if (hstar.lt.sqeps) call seteru (
       1  39hd9lgit  result less than half precision, 39, 1, 1)
  c
        d9lgit = -x - algap1 - dlog(hstar)
        return
  c
        end