Fast Newton-Raphson inverse_cdf for Chi and InverseGamma

[gaoflow]

Follows #390 / #382, per @YeungOnion's suggestion to extend the brent-like + Newton-Raphson inverse_cdf from Gamma (#382) to the distributions that still fall back to the generic default.

Chi and InverseGamma were the last two continuous distributions using the generic bisection default for inverse_cdf (the two named in #390). This gives each a custom solver in the same vein as Gamma.

Approach

A shared internal::newton_raphson_quantile(p, cdf, sf, pdf) brackets the quantile to a factor of two, then refines it with safeguarded Newton steps (Numerical Recipes' rtsafe): [low, high] is only ever tightened, and a Newton step that is non-finite or leaves the bracket falls back to bisection. Chi and InverseGamma call it with their own cdf/sf/pdf.

Two things were needed to keep it accurate and fast across the whole range, not just the middle:

• Relative convergence. fix: safeguard Gamma::inverse_cdf Newton step to avoid NaN for shape < 1 #382's prec::convergence folds in an absolute 1e-9 tolerance, which is meaningless when the quantile itself is ~1e-12, so the solver tests the relative Newton step instead. It also checks the raw Newton step for convergence before tightening the bracket — otherwise a converged step that rounds onto the endpoint just moved to x is rejected and the search bisects the whole span back (this showed up as ~30 wasted iterations near the median).
• sf in the upper half. As cdf saturates to one it can no longer resolve a deep upper-tail quantile (InverseGamma's reaches ~6e23 at p = 1 − 1e-12); inverting sf there keeps it well conditioned, as in Converge the default continuous inverse_cdf instead of a fixed 16 iterations #390.

Validation

Checked against scipy.stats, cross-checked with mpmath at 60 digits, over p ∈ [1e-12, 1 − 1e-12] for Chi(k), k ∈ {1,2,3,5,10} and InverseGamma(a,b), (a,b) ∈ {(1,1),(0.5,0.5),(2,1),(3,2),(5,3)}:

• worst-case relative error 3.4e-15; round-trip cdf(inverse_cdf(p)) ≈ p to ~3e-15.
• converges in a handful of Newton steps (medians ~5, vs ~48 fixed bisection steps) — roughly 60% fewer cdf/pdf evaluations across the grid, with the bisection fallback preserving the default's worst-case robustness.

The test_inverse_cdf_reference / test_inverse_cdf_round_trip cases added in #390 (scipy/mpmath values pinned in both tails) now exercise the new solver and still pass; a new internal::test_newton_raphson_quantile covers the shared routine directly against the Exponential closed form. cargo test, cargo fmt --check and cargo clippy are clean.

Gamma keeps its own #382 solver here; happy to fold it into this helper too as a follow-up if you'd prefer a single routine.

[codecov]

Codecov Report

❌ Patch coverage is 92.39130% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.75%. Comparing base (1029458) to head (c3fe60f).

Files with missing lines	Patch %	Lines
src/distribution/chi.rs	81.25%	3 Missing ⚠️
src/distribution/inverse_gamma.rs	81.25%	3 Missing ⚠️
src/distribution/internal.rs	98.33%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #404      +/-   ##
==========================================
- Coverage   94.77%   94.75%   -0.03%     
==========================================
  Files          61       61              
  Lines       13402    13494      +92     
==========================================
+ Hits        12702    12786      +84     
- Misses        700      708       +8     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.