Decide and enforce the quantile contract at extreme p

[OldCrow]

Found by the #46 characterization sweep (PR #101; docs/ACCURACY_CHARACTERIZATION.md, Findings + appendix).

Current behavior at deep-tail p (down to 1e-300) is inconsistent across families:

• gamma and chi-squared quantiles return NaN for some valid deep-tail p;
• pareto returns +inf and student-t −inf where finite doubles exist (e.g. student-t ν=1.001 at p=1e-300 has quantile ≈ −1.6e+299, representable);
• cauchy loses all relative accuracy at p=1e-300 (rel = 1.0);
• discrete quantiles differ from the right-continuous-inverse convention by one integer at exact lattice probabilities F(k) == p — decide and document which convention is the contract.

Decision needed first: what does the library promise at extreme p — best-representable finite value, documented saturation bounds, or explicitly unsupported range? The gamma/chi-squared NaN class also overlaps the parked corvus adoption items (#47/#52): inversion quality there is bounded by the incomplete-gamma/beta cores, so full accuracy in the tails likely waits on adoption even after the contract is settled.

Unmilestoned pending that decision.

🤖 Generated with Claude Code

[OldCrow]

Defensive review 2026-08-21 — the Cauchy row here is understated: the quantile degrades gradually from moderate p, not only at 1e-300. Shipped library vs mpmath (γ = 1, x₀ = 0): p = 1e-6 → rel 1.5e-11; 1e-10 → rel 5.1e-7; 1e-15 → rel 8.9e-3; 1e-300 → rel 1.0. The sweep's quantilePGrid() jumps from 1e-15 to 1e-300, so the 1e-10/1e-15 rows never appear in the appendix. Mechanism: tan(π(p − ½)) with the argument approaching −π/2.

The obvious x₀ − γ/tan(πp) rewrite is wrong near p = ½ (rel 4.3e-7 at p = 0.5 − 1e-10 — it moves the cancellation to the other end). The correct form is a split (verified ≤ 2.7e-16 at every p probed from 0.5 − 1e-10 down to 1e-300):

const double d = p - detail::HALF;
if (std::fabs(d) <= 0.25)
    return x0 + g * std::tan(detail::PI * d);          // d exact (Sterbenz); tan well-conditioned
return p < detail::HALF ? x0 - g / std::tan(detail::PI * p)
                        : x0 + g / std::tan(detail::PI * (detail::ONE - p));

Regression rows (γ = 1, x₀ = 0, rel ≤ 4e-16): p = 1e-6 → −318309.88618274347; 1e-10 → −3183098861.8379064; 1e-15 → −318309886183790.69; 1e-300 → −3.183098861837907e+299; 0.5 − 1e-10 → −3.1415929135263347e-10 (unchanged); symmetric rows at 1 − p with the sign flipped. Note this also settles the "best-representable finite value" half of the contract question for Cauchy: the split form reaches it.