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, dcos, dexp, dfac, dgamma, dgamr, dint,
     1  dlnrel, dlog, dsin, dsqrt
      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