fix: safeguard Gamma::inverse_cdf Newton step to avoid NaN for shape < 1

[gaoflow]

Summary

Gamma::inverse_cdf(p) returns NaN for any shape < 1 once p is small enough. ChiSquared (df < 2) and Erlang delegate to the same solver and inherit it — this is the NaN reported in #361 and #373.

Gamma::new(0.5, 1.0).unwrap().inverse_cdf(0.01);    // NaN — scipy: 7.8544e-5
Gamma::new(0.18, 0.18).unwrap().inverse_cdf(0.25);  // NaN — scipy: 1.6167e-3   (#373)
ChiSquared::new(1.0).unwrap().inverse_cdf(0.01);    // NaN — scipy: 1.5709e-4   (#361)

Root cause

inverse_cdf builds a valid bracket [low, high] by bisection and then refines it with an unguarded Newton–Raphson loop. For shape < 1 the true quantile sits extremely close to 0 while the post-bisection iterate is still far from it, so the first Newton step overshoots to x < 0, outside the support. There cdf(x) == 0 and pdf(x) == 0, so the update is (0 - p) / 0 = -inf; the next step is inf - inf = NaN, which propagates because the loop never re-brackets. The bracket invariant from the bisection phase is discarded by the Newton loop.

Fix

Merge the two phases into one safeguarded loop (Numerical Recipes' rtsafe): keep [low, high] as an invariant, tighten it by the sign of cdf(x) - p each step, and accept a Newton step only when it is finite and strictly inside the bracket, otherwise take a bisection step. The iterate can no longer leave the support, so the result is always finite; shape ≥ 1 results are unchanged.

I diff-tested the whole inverse_cdf surface against scipy.stats / mpmath (dps=60) to scope this: the Gamma family (Gamma/ChiSquared/Erlang) was the only one producing NaN. Beta/StudentsT/FisherSnedecor go through the separately-safeguarded AS 64 inv_beta_reg, and Chi/InverseGamma use the bracketed default search — none share this failure — so the fix stays in the Gamma solver.

Tests

Added regression tests in gamma.rs (all NaN/panic before the fix): shape < 1 quantiles are finite; values checked against scipy/mpmath constants; round-trip cdf(inverse_cdf(p)) ≈ p and monotonicity across the shape/p grid; and ChiSquared/Erlang delegation, including the large-df upper-tail case from #361. Full test suite, cargo fmt --check and cargo clippy are clean.

[YeungOnion]

I like the solution, we might benefit from this kind of search for other inverse cdf steps for float domains. Thanks for this!