program tgam
c
c -- the gamma function
c -- figure 11.11
c
      integer out
      real x, ans
c
      out = 6
10    write(out, 101)
      read(*, 102) x
      if (x .eq. 0.0) goto 10
      if (x .lt. -22.0) goto 99
      ans = gamma(x)
      write(out, 104) ans
      goto 10
99    stop
101   format(/' arg? ' )
102   format(e10.0)
104   format(' gamma =', f11.4)
      end
      function gamma(x)
c
c -- the gamma functionction by infinite series
c
      integer i, j
      real x, y, z, pi, gam
      data pi/ 3.141593/
c
c -- function within a function
c
      g(y) = sqrt(2 * pi/y) *
     *  exp(y * alog(y) + (1 - 1/(30*y*y))/(12*y)-y)
     *  / ((y - 2) * (y - 1))
c
      if (x .lt. 0.0) goto 10
      gamma = g(x + 2.0)
      goto 99
c
c -- increment argument until positive
c
10    z = x
      j = -1
      y = x
20     j = j + 1
       y = y + 1.0
       if (y .lt. 0.0) goto 20
      gam = g(y + 2.0)
      do 30 i = 0, j
       gam = gam / (x + i)
30    continue
      gamma = gam
99    return
      end