function gamit (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
      external aint, alngam, alog, exp, gamr, r1mach, r9gmit, r9lgic,
     1  r9lgit, sqrt
      data alneps, sqeps, bot / 3*0.0 /
c
      if (alneps.ne.0.0) go to 10
      alneps = -alog(r1mach(3))
      sqeps = sqrt(r1mach(4))
      bot = alog(r1mach(1))
c
 10   if (x.lt.0.0) call seteru (21hgamit   x is negative, 21, 2, 2)
c
      if (x.ne.0.0) alx = alog(x)
      sga = 1.0
      if (a.ne.0.0) sga = sign (1.0, a)
      ainta = aint (a+0.5*sga)
      aeps = a - ainta
c
      if (x.gt.0.0) go to 20
      gamit = 0.0
      if (ainta.gt.0.0 .or. aeps.ne.0.0) gamit = gamr(a+1.0)
      return
c
 20   if (x.gt.1.0) go to 40
      if (a.ge.(-0.5) .or. aeps.ne.0.0) call algams (a+1.0, algap1,
     1  sgngam)
      gamit = r9gmit (a, x, algap1, sgngam, alx)
      return
c
 40   if (a.lt.x) go to 50
      t = r9lgit (a, x, alngam(a+1.0))
      if (t.lt.bot) call erroff
      gamit = exp(t)
      return
c
 50   alng = r9lgic (a, x, alx)
c
c evaluate gamit in terms of alog(gamic(a,x))
c
      h = 1.0
      if (aeps.eq.0.0 .and. ainta.le.0.0) go to 60
      call algams (a+1.0, algap1, sgngam)
      t = alog(abs(a)) + alng - algap1
      if (t.gt.alneps) go to 70
      if (t.gt.(-alneps)) h = 1.0 - sga*sgngam*exp(t)
      if (abs(h).gt.sqeps) go to 60
      call erroff
      call seteru (32hgamit   result lt half precision, 32, 1, 1)
c
 60   t = -a*alx + alog(abs(h))
      if (t.lt.bot) call erroff
      gamit = sign (exp(t), h)
      return
c
 70   t = t - a*alx
      if (t.lt.bot) call erroff
      gamit = -sga*sgngam*exp(t)
      return
c
      end