double precision function dpoch (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
        double precision a, x, absa, absax, alnga, alngax, ax, b, pi,
       1  eps, sqeps, cospix, sinpix, cospia, sinpia, errpch, den, err,
       2  sgnga, sgngax, dfac, dlnrel, d9lgmc, dgamma, dgamr, dint,
       3  d1mach, dcos, dexp, dlog, dsin, dsqrt
        external d1mach, d9lgmc, dfac, dgamma, dgamr,
       1  dlnrel
        data pi / 3.1415926535 8979323846 2643383279 503 d0 /
        data eps, sqeps / 2*0.0d0 /
  c
        if (eps.ne.0.0d0) go to 10
        eps = d1mach(4)
        sqeps = dsqrt(eps)
  c
   10   ax = a + x
        if (ax.gt.0.0d0) go to 30
        if (dint(ax).ne.ax) go to 30
  c
        if (a.gt.0.0d0 .or. dint(a).ne.a) call seteru (
       1  48hdpoch   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
        dpoch = 1.0d0
        if (x.eq.0.d0) return
  c
        n = x
        if (dmin1(a+x,a).lt.(-20.0d0)) go to 20
  c
        ia = a
        dpoch = (-1.0d0)**n * dfac(-ia)/dfac(-ia-n)
        return
  c
   20   dpoch = (-1.0d0)**n * dexp ((a-0.5d0)*dlnrel(x/(a-1.0d0))
       1  + x*dlog(-a+1.0d0-x) - x + d9lgmc(-a+1.0d0) - d9lgmc(-a-x+1.d0))
        return
  c
  c a+x is not zero or a negative integer.
  c
   30   dpoch = 0.0d0
        if (a.le.0.0d0 .and. dint(a).eq.a) return
  c
        n = dabs(x)
        if (dble(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
        dpoch = 1.0d0
        if (n.eq.0) return
        do 40 i=1,n
          dpoch = dpoch * (a+dble(float(i-1)))
   40   continue
        return
  c
   50   absax = dabs(a+x)
        absa = dabs(a)
        if (dmax1(absax,absa).gt.20.0d0) go to 60
  c
        if (dabs(x).gt.absax/sqeps) call seteru (67hdpoch   answer lt half
       1 precision, a+x cannot be computed accurately, 67, 1, 1)
        dpoch = dgamma(a+x) * dgamr(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 (dabs(x).gt.0.5d0*absa) go to 70
  c
  c 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.0d0) b = -a - x + 1.0d0
        dpoch = dexp ((b-0.5d0)*dlnrel(x/b) + x*dlog(b+x) - x
       1  + d9lgmc(b+x) - d9lgmc(b) )
  c
        if (a.ge.0.0d0 .or. dpoch.eq.0.0d0) return
        call entsrc (irold, 1)
        cospix = dcos (pi*x)
        call erroff
        sinpix = dsin (pi*x)
        call erroff
        cospia = dcos (pi*a)
        call erroff
        sinpia = dsin (pi*a)
        call erroff
        call entsrc (irold2, irold)
  c
        errpch = dabs(x)*(1.0d0+dlog(b))
        den = cospix + cospia*sinpix/sinpia
        err = (dabs(x)*(dabs(sinpix) + dabs(cospia*cospix/sinpia))
       1  + dabs(a*sinpix)/sinpia**2)*pi
        err = errpch + err/dabs(den)
        if (err.gt.1.0d0/eps) call seteru (73hdpoch   answer has no precis
       1ion, a or a+x too close to a negative integer, 73, 3, 2)
        if (err.gt.1.0d0/sqeps) call seteru (74hdpoch   answer lt half pre
       2cision, a or a+x too close to a negative integer, 74, 2, 1)
  c
        dpoch = dpoch/den
        return
  c
   70   call dlgams (a+x, alngax, sgngax)
        call dlgams (a, alnga, sgnga)
        dpoch = sgngax * sgnga * dexp(alngax-alnga)
  c
        return
        end