double precision function d9gmit (a, x, algap1, sgngam, alx)
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c compute tricomi-s incomplete gamma function for small x.
  c
        double precision a, x, algap1, sgngam, alx, ae, aeps, algs, alg2,
       1  bot, eps, fk, s, sgng2, t, te, d1mach, dlngam, dexp, dlog
        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 (x.le.0.d0) call seteru (24hd9gmit  x should be gt 0, 24, 1, 2)
  c
        ma = a + 0.5d0
        if (a.lt.0.d0) ma = a - 0.5d0
        aeps = a - dble(float(ma))
  c
        ae = a
        if (a.lt.(-0.5d0)) ae = aeps
  c
        t = 1.d0
        te = ae
        s = t
        do 20 k=1,200
          fk = k
          te = -x*te/fk
          t = te/(ae+fk)
          s = s + t
          if (dabs(t).lt.eps*dabs(s)) go to 30
   20   continue
        call seteru (54hd9gmit  no convergence in 200 terms of taylor-s se
       1ries, 54, 2, 2)
  c
   30   if (a.ge.(-0.5d0)) algs = -algap1 + dlog(s)
        if (a.ge.(-0.5d0)) go to 60
  c
        algs = -dlngam(1.d0+aeps) + dlog(s)
        s = 1.0d0
        m = -ma - 1
        if (m.eq.0) go to 50
        t = 1.0d0
        do 40 k=1,m
          t = x*t/(aeps-dble(float(m+1-k)))
          s = s + t
          if (dabs(t).lt.eps*dabs(s)) go to 50
   40   continue
  c
   50   d9gmit = 0.0d0
        algs = -dble(float(ma))*dlog(x) + algs
        if (s.eq.0.d0 .or. aeps.eq.0.d0) go to 60
  c
        sgng2 = sgngam * dsign (1.0d0, s)
        alg2 = -x - algap1 + dlog(dabs(s))
  c
        if (alg2.gt.bot) d9gmit = sgng2 * dexp(alg2)
        if (algs.gt.bot) d9gmit = d9gmit + dexp(algs)
        return
  c
   60   d9gmit = dexp (algs)
        return
  c
        end