fix: evaluate the Lanczos sum in a well-conditioned form

[agene0001]

Stacked on #410 (needs its tightened ulps_eq! for the digamma pole test); review only the last commit until that merges.

ln_gamma/gamma evaluate Pugh's Lanczos approximation as the partial-fraction sum d₀ + Σ dₖ/(z+k−1). The residues alternate in sign, so the sum cancels badly — condition number ~3600 around z = 50, i.e. ~440 eps of relative error in the sum alone.

The same rational function written as a single fraction N(z)/D(z), derived from those residues in exact rational arithmetic, has all-positive coefficients in both polynomials (D(z) = z(z+1)…(z+9), whose expanded coefficients are unsigned Stirling numbers of the first kind, exact in f64). For z > 0 both Horner evaluations then have condition number 1. No new external constants — they are derived from the residues already in the file, and agree with the partial-fraction form to 3.2e-17 over [0.5, 3000]. A reversed form in 1/z covers z ≥ 1e29 where z¹⁰ would overflow.

Also in the same evaluation:

• lanczos_power compensates the base of ((p+g)/e)^p — powf amplifies base error by the exponent, the dominant remaining term (~190 ulp by x = 122). The correction is first-order and gated on staying small, since applied blindly it can flip the sign of an overflowing result.
• Integer arguments ≤ 171 come from the exact factorial table (Γ(2) == 1 exactly, was 1.0000000000000002).
• gamma above x ≈ 169.7 halves the exponent and squares: the power alone overflows there while the product is representable to x ≈ 171.61. gamma(171.6) was inf, now 1.5858969096673e308 (mpmath: 1.5858969096673029e308).
• sin_pi/tan_pi reduce by the period before sin/tan — x − round(x) is exact, so reflection-formula error stays relative to the fractional part instead of ulp(π)/dist_to_pole.

Measured against mpmath at 45 dps:

	before	after
ln_gamma p99	45.9 ulp	8.2 ulp
gamma (positive), median	285 ulp	19 ulp
gamma (negative), median/max	109 / 1173 ulp	3.4 / 22.8 ulp
digamma (negative), median/max	19.5 / 1.5e10 ulp	1.6 / 538 ulp

The remaining ~19–35 ulp in gamma is the approximation floor of the f64-rounded Pugh coefficient set itself (~3.7e-15 — confirmed by evaluating the formula in exact arithmetic), so evaluation error is now below the fit.

Four test expectations move as a consequence, each verified against mpmath: one Binomial::sf and two FisherSnedecor literals were fitted to the old output, and the Dirichlet/MultivariateStudent doctest values were 34 and 22 ulp from truth (new outputs are within 2).

🤖 Generated with Claude Code