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