double precision function d9gmic (a, x, alx)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c compute the complementary incomplete gamma function for a near
c a negative integer and for small x.
c
      double precision a, x, alx, alng, bot, eps, euler, fk, fkp1, fm,
     1  s, sgng, t, te, d1mach, dlngam, dexp, dlog
      external d1mach, dexp, dlngam, dlog
      data euler / 0.5772156649 0153286060 6512090082 40 d0 /
      data eps, bot / 2*0.d0 /
c
      if (eps.ne.0.d0) go to 10
      eps = 0.5d0*d1mach(3)
      bot = dlog (d1mach(1))
c
 10   if (a.gt.0.d0) call seteru (
     1  41hd9gmic  a must be near a negative integer, 41, 2, 2)
      if (x.le.0.d0) call seteru (25hd9gmic  x must be gt zero, 25, 3,2)
c
      m = -(a - 0.5d0)
      fm = m
c
      te = 1.0d0
      t = 1.0d0
      s = t
      do 20 k=1,200
        fkp1 = k + 1
        te = -x*te/(fm+fkp1)
        t = te/fkp1
        s = s + t
        if (dabs(t).lt.eps*s) go to 30
 20   continue
      call seteru (57hd9gmic  no convergence in 200 terms of continued f
     1raction, 57, 4, 2)
c
 30   d9gmic = -alx - euler + x*s/(fm+1.0d0)
      if (m.eq.0) return
c
      if (m.eq.1) d9gmic = -d9gmic - 1.d0 + 1.d0/x
      if (m.eq.1) return
c
      te = fm
      t = 1.d0
      s = t
      mm1 = m - 1
      do 40 k=1,mm1
        fk = k
        te = -x*te/fk
        t = te/(fm-fk)
        s = s + t
        if (dabs(t).lt.eps*dabs(s)) go to 50
 40   continue
c
 50   do 60 k=1,m
        d9gmic = d9gmic + 1.0d0/dble(float(k))
 60   continue
c
      sgng = 1.0d0
      if (mod(m,2).eq.1) sgng = -1.0d0
      alng = dlog(d9gmic) - dlngam(fm+1.d0)
c
      d9gmic = 0.d0
      if (alng.gt.bot) d9gmic = sgng * dexp(alng)
      if (s.ne.0.d0) d9gmic = d9gmic +
     1  dsign (dexp(-fm*alx+dlog(dabs(s)/fm)), s)
c
      if (d9gmic.eq.0.d0 .and. s.eq.0.d0) call seteru (
     1  25hd9gmic  result underflows, 25, 1, 0)
      return
c
      end