Module pochhammer 

The Pochhammer symbol $(z)_m = \Gamma(z+m)/\Gamma(z)$.

Three paths, and why the obvious one is not enough

Written out, this is a ratio of two Gamma functions, and the obvious spelling exp(lgamma(z+m) - lgamma(z)) is a disaster in exactly the region the function is most used. The two logarithms are large and nearly equal whenever m is small next to z, so the subtraction sheds the digits that carry the answer: measured against mpmath at z = 1e8, m = 1e-4, that form has no correct digits at all (2.8e-7 relative, where the true value is within 1e-3 of 1). Everything below exists to avoid forming that difference.

Integer m, small (|m| <= PRODUCT_CAP). The definition collapses to a plain product z(z+1)...(z+m-1), which forms no logarithm at all and is therefore exact to within its own multiplications: measured worst 4.1 ulp across the sweep, and 0 to 0.2 ulp on most of it. This is the dominant case, not a fast path bolted on. Hypergeometric series, binomial-style coefficients and Taylor coefficients of special functions all advance m by whole numbers. A negative integer m is the reciprocal of the same product started at z + m, which is why the sign of m only chooses a starting point and a final reciprocal. The product is also indifferent to the sign of z, so this path covers the negative half of the domain for free, poles included: (-2)_3 contains a zero factor and correctly returns 0.

It runs at Average and above, and within that is gated at runtime on any lane wanting it, returning early when every lane does. That is the shape gamma’s exact-integer branch uses, one tier lower. The tier differs because the trade does: for gamma the Lanczos path is fast and runs regardless, so its integer branch is pure accuracy spend, while here the product is also the cheaper route for integer-heavy data (a handful of multiplies and an early return, against the full Stirling evaluation). With the runtime guard skipping it outright when no lane wants it, the only shape that pays for having it is a genuinely mixed vector.

Below Average it is compiled out and integer m goes through the Stirling difference like anything else. Measured over 231 points with m in 0..20 and z across eleven magnitudes, that is 4.2 ulp median and 172 worst, against 0.00 median and 4.2 worst for the product. The visible difference is the exactness rather than the ulp count: (3)_1 is 3.0 on the product path and 3.0000000000000018 without it.

Everything else with both arguments positive. Take the Stirling difference instead of the logarithm difference. With $\ln\Gamma(x) = (x - \tfrac12)\ln x - x + \tfrac12\ln 2\pi + \mathrm{stirlerr}(x)$, the $\tfrac12 \ln 2\pi$ cancels exactly and the rest regroups so that nothing large is ever subtracted from anything large:

\ln\frac{\Gamma(x+m)}{\Gamma(x)}
  = \left(x - \tfrac12\right)\ln\!\left(1 + \frac{m}{x}\right)
  + m\left(\ln(x+m) - 1\right)
  + \mathrm{stirlerr}(x+m) - \mathrm{stirlerr}(x)

Every term is O(m) as m -> 0, which is what makes the small-m region well behaved. The log1p is doing the work the naive subtraction failed at. stirlerr is only valid at or above STIRLERR_MIN, so an argument below it is first walked up by a whole number of steps and the exact product of those steps divided back out, the same shifted-product trick pmf_parts uses, and for the same reason.

Each argument is shifted independently, which is what keeps every intermediate in range: a product is built only for an argument that is below 9, so its factors are under 18 and it can never exceed 18^9. Shifting both by a shared amount instead (the obvious spelling) walks an argument that was already fine and overflows it, and the resulting inf * 0 is a NaN sitting exactly where the answer is a perfectly good infinity. The independent shift is also the more accurate of the two, because the product it skips is a string of roundings that never happens. See the comments on the shift for the numbers.

Its accuracy is the floor of anything that exponentiates a logarithm: the relative error of the result is the absolute error of the exponent, so it tracks $|\ln (z)_m| \cdot \epsilon$ and is bounded below by nothing else. Measured against mpmath over 6924 points with z in [0.1, 8.9] and non-integer m, the median is 2.6 ulp, the 99th percentile 25 ulp and the worst 51 ulp. Individual points scale with the result’s own logarithm, reaching 259 ulp at z = 3.7, m = 100 where the value is near 1e163 and |ln| = 375. It falls to zero error where the answer approaches 1, which is precisely where the naive form was worst.

The residue. A non-integer m (or one past the cap) with z or z + m non-positive reaches neither path above, and falls back to the logarithmic form with the sign taken from lgamma_r. It inherits that form’s cancellation. This is the region where (z)_m is a ratio across Gamma’s poles and no cheap rearrangement is available. It is documented rather than fixed.

Constants

PRODUCT_CAP

Largest |m| taken by the exact product path.

Functions

pochhammer

$(z)_m = \Gamma(z+m)/\Gamma(z)$, the Pochhammer symbol, for real z and real m.