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