fix: use Loader's saddle-point form for the gamma/beta prefixes and pmfs

[agene0001]

Stacked on #417 (beta_reg robustness); review only the last commit until that merges.

Densities and incomplete-function prefixes in the gamma/beta family are written as exp(a·ln x − x − ln_gamma(a)) and friends. Every term grows like a·ln a while the sum stays O(ln a), so at a = 1e4 each term is ~1e5 and the ~1e-11 absolute error of the sum becomes the accuracy ceiling of everything downstream — this is what limited beta_reg to ~6e-11 even where its recurrence converged, and put Binomial::pmf at ~7e-11 relative.

Adds the two kernels from Loader, Fast and Accurate Computation of Binomial Probabilities (2000):

• stirling_delta(z) — the Stirling-series remainder of ln_gamma, O(1/12z). Table-free: the recurrence δ(z) = δ(z+1) + (z+½)·ln1p(1/z) − 1 lifts any argument to where a 6-term series is good to 1.4e-18, so there is no interpolation table and no recursive dependency on ln_gamma. ≤1.2e-15 absolute, which is what matters since every caller adds it to an exponent.
• bd0(x, np) = x·ln(x/np) − (x − np), with a double root at x = np. Evaluated by a series in (x−np)/(x+np) so the cancellation happens analytically: ~1e-16 relative uniformly through the root.

Rewriting the prefixes (gamma_lr/gamma_ur, beta_reg, Binomial::pmf, Poisson::pmf) in these terms keeps every piece O(1):

	before	after
beta_reg, a = b = 1e4	6e-11	8.7e-15
beta_reg, a = b = 1e6	3e-9	2.6e-13
gamma_ur, a ≥ 16	5.3 median ulp, p99 970k	0.38 median, p99 3071
Binomial::pmf	685k median ulp	2.45
Poisson::pmf	649 median ulp	3.16

Two details reviewers should look at:

1. Both prefixes fall back to the direct form for small parameters (max param < 16): there's no cancellation to remove there, and without the gate FisherSnedecor(1,1).cdf got worse (the recurrence costs 16 steps at a = 0.5).
2. bd0 is sensitive enough that the rounding of n·p itself matters (∂bd0/∂np = 1 − x/np: half an ulp of np at 6e5 moves the result 2.5e-13 ≈ 1000 ulp of the pmf), so n, 1−p, and the products are carried as double-doubles via bd0_dd. R's dbinom has this same limitation.

Also guards bd0's direct branch where (x/np).ln() under/overflows — bd0(1e-300, 5e99) gave −inf where the answer is +5e99, which surfaced as NaN (or a silent 1.0 where the truth is ~1e-30118) in beta_reg at extreme parameter ratios.

This is the same log-space reasoning as #381, applied to the rest of the family.

🤖 Generated with Claude Code

[coderabbitai]

Important

Draft PR not reviewed

Draft PRs are not automatically reviewed by default.

• Trigger a manual review

To automatically review draft PRs, update your CodeRabbit configuration:

reviews:
  auto_review:
    drafts: true

---

_{Comment @coderabbitai help to get the list of available commands.}