double precision function dgamic (a, x)
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c evaluate the complementary incomplete gamma function
  c
  c gamic = integral from t = x to infinity of exp(-t) * t**(a-1.)  .
  c
  c gamic is evaluated for arbitrary real values of a and for non-negative
  c values of x (even though gamic is defined for x .lt. 0.0), except that
  c for x = 0 and a .le. 0.0, gamic is undefined.
  c
  c      a slight deterioration of 2 or 3 digits accuracy will occur when
  c gamic 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, alngs, alx,
       1  bot, e, eps, gstar, h, sga, sgng, sgngam, sgngs, sqeps, t,
       2  d1mach, dlngam, d9gmic, d9gmit, d9lgic, d9lgit, dint,
       3  dexp, dlog, dsqrt
        external d1mach, d9gmic, d9gmit, d9lgic, d9lgit
       1  dlngam
  c
        data eps, sqeps, alneps, bot / 4*0.0d0 /
  c
        if (eps.ne.0.d0) go to 10
        eps = 0.5d0*d1mach(3)
        sqeps = dsqrt (d1mach(4))
        alneps = -dlog (d1mach(3))
        bot = dlog (d1mach(1))
  c
   10   if (x.lt.0.d0) call seteru (21hdgamic  x is negative, 21, 2, 2)
  c
        if (x.gt.0.d0) go to 20
        if (a.le.0.d0) call seteru (
       1  47hdgamic  x = 0 and a le 0 so dgamic is undefined, 47, 3, 2)
  c
        dgamic = dexp (dlngam(a+1.d0) - dlog(a))
        return
  c
   20   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
        izero = 0
        if (x.ge.1.0d0) go to 40
  c
        if (a.gt.0.5d0 .or. dabs(aeps).gt.0.001d0) go to 30
        e = 2.0d0
        if (-ainta.gt.1.d0) e = 2.d0*(-ainta+2.d0)/(ainta*ainta-1.0d0)
        e = e - alx * x**(-0.001d0)
        if (e*dabs(aeps).gt.eps) go to 30
  c
        dgamic = d9gmic (a, x, alx)
        return
  c
   30   call dlgams (a+1.0d0, algap1, sgngam)
        gstar = d9gmit (a, x, algap1, sgngam, alx)
        if (gstar.eq.0.d0) izero = 1
        if (gstar.ne.0.d0) alngs = dlog (dabs(gstar))
        if (gstar.ne.0.d0) sgngs = dsign (1.0d0, gstar)
        go to 50
  c
   40   if (a.lt.x) dgamic = dexp (d9lgic(a, x, alx))
        if (a.lt.x) return
  c
        sgngam = 1.0d0
        algap1 = dlngam (a+1.0d0)
        sgngs = 1.0d0
        alngs = d9lgit (a, x, algap1)
  c
  c evaluation of dgamic(a,x) in terms of tricomi-s incomplete gamma fn.
  c
   50   h = 1.d0
        if (izero.eq.1) go to 60
  c
        t = a*alx + alngs
        if (t.gt.alneps) go to 70
        if (t.gt.(-alneps)) h = 1.0d0 - sgngs*dexp(t)
  c
        if (dabs(h).lt.sqeps) call erroff
        if (dabs(h).lt.sqeps) call seteru (
       1  32hdgamic  result lt half precision, 32, 1, 1)
  c
   60   sgng = dsign (1.0d0, h) * sga * sgngam
        t = dlog(dabs(h)) + algap1 - dlog(dabs(a))
        if (t.lt.bot) call erroff
        dgamic = sgng * dexp(t)
        return
  c
   70   sgng = -sgngs * sga * sgngam
        t = t + algap1 - dlog(dabs(a))
        if (t.lt.bot) call erroff
        dgamic = sgng * dexp(t)
        return
  c
        end