double precision function dgamit (a, x)
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c evaluate tricomi-s incomplete gamma function defined by
  c
  c gamit = x**(-a)/gamma(a) * integral t = 0 to x of exp(-t) * t**(a-1.)
  c
  c and analytic continuation for a .le. 0.0.  gamma(x) is the complete
  c gamma function of x.  gamit is evaluated for arbitrary real values of
  c a and for non-negative values of x (even though gamit is defined for
  c x .lt. 0.0).
  c
  c      a slight deterioration of 2 or 3 digits accuracy will occur when
  c gamit is very large or very small in absolute value, because log-
  c arithmic variables are used.  also, if the parameter a is very close
  c to a negative integer (but not a negative integer), there is a loss
  c of accuracy, which is reported if the result is less than half
  c machine precision.
  c
  c ref. -- w. gautschi, an evaluation procedure for incomplete gamma
  c functions, acm trans. math. software.
  c
        double precision a, x, aeps, ainta, algap1, alneps, alng, alx,
       1  bot, h, sga, sgngam, sqeps, t, d1mach, dgamr, d9gmit, d9lgit,
       2  dlngam, d9lgic, dint, dexp, dlog, dsqrt
        external d1mach, d9gmit, d9lgic, d9lgit, dgamr,
       1  dlngam
  c
        data alneps, sqeps, bot / 3*0.d0 /
  c
        if (alneps.ne.0.d0) go to 10
        alneps = -dlog (d1mach(3))
        sqeps = dsqrt (d1mach(4))
        bot = dlog (d1mach(1))
  c
   10   if (x.lt.0.d0) call seteru (21hdgamit  x is negative, 21, 2, 2)
  c
        if (x.ne.0.d0) alx = dlog (x)
        sga = 1.0d0
        if (a.ne.0.d0) sga = dsign (1.0d0, a)
        ainta = dint (a + 0.5d0*sga)
        aeps = a - ainta
  c
        if (x.gt.0.d0) go to 20
        dgamit = 0.0d0
        if (ainta.gt.0.d0 .or. aeps.ne.0.d0) dgamit = dgamr(a+1.0d0)
        return
  c
   20   if (x.gt.1.d0) go to 30
        if (a.ge.(-0.5d0) .or. aeps.ne.0.d0) call dlgams (a+1.0d0, algap1,
       1  sgngam)
        dgamit = d9gmit (a, x, algap1, sgngam, alx)
        return
  c
   30   if (a.lt.x) go to 40
        t = d9lgit (a, x, dlngam(a+1.0d0))
        if (t.lt.bot) call erroff
        dgamit = dexp (t)
        return
  c
   40   alng = d9lgic (a, x, alx)
  c
  c evaluate dgamit in terms of dlog (dgamic (a, x))
  c
        h = 1.0d0
        if (aeps.eq.0.d0 .and. ainta.le.0.d0) go to 50
  c
        call dlgams (a+1.0d0, algap1, sgngam)
        t = dlog (dabs(a)) + alng - algap1
        if (t.gt.alneps) go to 60
  c
        if (t.gt.(-alneps)) h = 1.0d0 - sga * sgngam * dexp(t)
        if (dabs(h).gt.sqeps) go to 50
  c
        call erroff
        call seteru (32hdgamit  result lt half precision, 32, 1, 1)
  c
   50   t = -a*alx + dlog(dabs(h))
        if (t.lt.bot) call erroff
        dgamit = dsign (dexp(t), h)
        return
  c
   60   t = t - a*alx
        if (t.lt.bot) call erroff
        dgamit = -sga * sgngam * dexp(t)
        return
  c
        end