function poch (a, x)
c august 1980 edition.  w. fullerton, c3, los alamos scientific lab.
c error handling when poch (a, x) is less than half precision is
c probably incorrect.  grossly wrong arguments are not handled right.
c
c evaluate a generalization of pochhammer-s symbol
c (a)-sub-x = gamma(a+x)/gamma(a).  for x a non-negative integer,
c poch(a,x) is just pochhammer-s symbol.
c
      external  alnrel, fac, gamma, gamr, r1mach,
     1  r9lgmc
      data pi / 3.1415926535 89793238 e0 /
      data eps, sqeps / 2*0.0 /
c
      if (eps.ne.0.0) go to 10
      eps = r1mach(4)
      sqeps = sqrt(eps)
c
 10   ax = a + x
      if (ax.gt.0.0) go to 30
      if (aint(ax).ne.ax) go to 30
c
      if (a.gt.0.0 .or. aint(a).ne.a) call seteru (
     1  48hpoch    a+x is non-positive integer but a is not, 48, 2, 2)
c
c we know here that both a+x and a are non-positive integers.
c
      poch = 1.0
      if (x.eq.0.0) return
c
      n = x
      if (amin1(a+x,a).lt.(-20.0)) go to 20
c
      poch = (-1.0)**n * fac(-int(a))/fac(-int(a)-n)
      return
c
 20   poch = (-1.0)**n * exp ((a-0.5)*alnrel(x/(a-1.0))
     1  + x*alog(-a+1.0-x) - x + r9lgmc(-a+1.) - r9lgmc(-a-x+1.) )
      return
c
c here we know a+x is not zero or a negative integer.
c
 30   poch = 0.0
      if (a.le.0.0 .and. aint(a).eq.a) return
c
      n = abs(x)
      if (float(n).ne.x .or. n.gt.20) go to 50
c
c x is a small non-positive integer, presummably a common case.
c
      poch = 1.0
      if (n.eq.0) return
      do 40 i=1,n
        poch = poch * (a+float(i-1))
 40   continue
      return
c
 50   absax = abs(a+x)
      absa = abs(a)
      if (amax1(absax,absa).gt.20.0) go to 60
c
      if (abs(x).gt.absax/sqeps) call seteru (67hpoch    answer lt half
     1precision, a+x cannot be computed accurately, 67, 1, 1)
      poch = gamma(a+x)*gamr(a)
c error handling above is probably bad when a almost = -n and when x is
c small.  maybe use reflection formula below in modified form.
      return
c
 60   if (abs(x).gt.0.5*absa) go to 70
c
c here abs(x) is small and both abs(a+x) and abs(a) are large.  thus,
c a+x and a must have the same sign.  for negative a, we use
c gamma(a+x)/gamma(a) = gamma(-a+1)/gamma(-a-x+1) *
c sin(pi*a)/sin(pi*(a+x))
c
      b = a
      if (b.lt.0.0) b = -a - x + 1.0
      poch = exp ((b-0.5)*alnrel(x/b) + x*alog(b+x) - x +
     1  r9lgmc(b+x) - r9lgmc(b) )
c
      if (a.ge.0.0 .or. poch.eq.0.0) return
      call entsrc (irold, 1)
      cospix = cos (pi*x)
      call erroff
      sinpix = sin (pi*x)
      call erroff
      cospia = cos (pi*a)
      call erroff
      sinpia = sin (pi*a)
      call erroff
      call entsrc (irold2, irold)
c
      errpch = abs(x)*(1.0+alog(b))
      den = cospix + cospia*sinpix/sinpia
      err = (abs(x)*(abs(sinpix)+abs(cospia*cospix/sinpia))
     1  + abs(a*sinpix)/sinpia**2)*pi
      err = errpch + err/abs(den)
      if (err.gt.1.0/eps) call seteru (73hpoch    answer has no precisio
     1n, a or a+x too close to a negative integer, 73, 3, 2)
      if (err.gt.1.0/sqeps) call seteru (72hpoch    answer lt half preci
     1on, a or a+x too close to a negative integer, 72, 2, 1)
c
      poch = poch/den
      return
c
 70   call algams (a+x, alngax, sgngax)
      call algams (a, alnga, sgnga)
      poch = sgngax * sgnga * exp(alngax-alnga)
c
      return
      end