Module poisson 

Loader’s saddle-point pieces for densities of the shape $x^k e^{-x}/\Gamma(k+1)$.

The Poisson mass $e^{-\lambda}\lambda^k/k!$, the Gamma density, and the seed of the orthonormal Laguerre functions are all this shape, and the obvious spelling exp(k ln lambda - lgamma(k+1) - lambda) computes an O(1) answer as the exponential of a difference of large terms: half an ulp of lgamma(k+1) = O(k ln k) becomes that many ulp of the result. Loader (2000, “Fast and accurate computation of binomial probabilities”, the form R’s dpois uses) rewrites it as

\frac{\lambda^k e^{-\lambda}}{k!} = \frac{e^{-\mathrm{stirlerr}(k) - \mathrm{bd0}(k, \lambda)}}{\sqrt{2\pi k}}

with the two pieces below, both small where the density is not negligible, so the exponential amplifies nothing.

Constants

STIRLERR_MIN

Below this the Stirling series in stirlerr is not accurate to binary64 at any depth (it is asymptotic, and the smallest term at n = 9 is under 1e-18, at n = 6 it is 1e-14). Callers handle n < STIRLERR_MIN some other way: a table for integers, or Gamma(n+1) directly, which is cheap and well conditioned at small argument.

Functions

bd0_series

The binomial/Poisson deviance bd0(k, lambda) = k ln(k/lambda) + lambda - k >= 0 in its peak form: (k - lambda) v + 2k sum_{j>=1} v^{2j+1}/(2j+1) for v = (k - lambda)/(k + lambda) (from ln((1+v)/(1-v)) = 2 atanh v), bd0_terms terms, full precision inside |v| < 1/5. Away from the peak the direct form is fine and pmf_parts uses it (folded with the rest of the exponent); this is only the part that needs care. diff and v are passed in because callers have them.

bd0_terms

Terms of the odd series in [bd0] by tier, inside |v| < 1/5.

exp_sum

exp(base + rest_hi + rest_lo) where base is a large exact-ish number (-lambda, x/4) and the rest pair is small: TwoSum recovers the rounding of the sum, and e^{s + lo} = e^s (1 + lo) to first order. Without it the sum rounds to half an ulp of base, which the exponential turns into hundreds of ulp.

pmf_parts

The shared core of every $x^k e^{-x}/\Gamma(k+1)$ shape here: the Poisson mass, its log, and the Laguerre-function seed. Returns (rest_hi, rest_lo, large, prod, n) such that

poisson_pmf

The Poisson mass $e^{-\lambda}\lambda^k/k!$ at real k >= 0, lambda >= 0 (LOG = false), or its log (LOG = true), through pmf_parts. Real k because the Gamma density is the same function (dgamma(x; a) = pmf(a - 1; x) for a >= 1).

stirlerr

Stirling’s error $\mathrm{stirlerr}(n) = \ln n! - \left[(n + \tfrac12)\ln n - n + \tfrac12 \ln 2\pi\right]$, for n >= STIRLERR_MIN, by the Bernoulli series

stirlerr_terms

Terms of the 1/n^2 series in stirlerr by tier, at n >= STIRLERR_MIN. The dropped term bounds the absolute error, which is the relative error of whatever density it feeds: 9 terms is 1e-18, 6 is 2.5e-15 (a measured 5 ulp at n = 9, so Average keeps all 9, the three FMAs being nothing), 5 is 1e-13, and the single 1/(12n) term is 3.8e-6. Worst forgoes the series and keeps just that.