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