Fix Gamma/ChiSquared inverse-CDF collapse in the small-shape lower tail

[gaoflow]

Problem

Gamma.InvCDF and ChiSquared.InvCDF collapse in the small-shape lower tail, returning a spurious ~1.11e-15 floor and disagreeing with the library's own CDF:

Gamma.InvCDF(0.1, 1, 1e-6);   // returns 1.1102230246251567E-15
Gamma.InvCDF(0.1, 1, 1e-3);   // returns the SAME floored value
Gamma.CDF(0.1, 1, Gamma.InvCDF(0.1, 1, 1e-6));  // = 0.03359, not 1e-6

ChiSquared.InvCDF(0.2, 1e-6); // 2.22e-15 ; ChiSquared.CDF(0.2, that) = 0.0336

The true quantile of Gamma(0.1) at p = 1e-6 is ~6.07e-61 (mpmath, dps 40) — the returned value is wrong by ~46 orders of magnitude, and it fails the round-trip identity CDF(InvCDF(p)) == p by ~4 orders, which is enough to call it a bug without any external reference. It degrades with shape: a = 0.5 is fine, a = 0.3 is off ~36x, a = 0.1 is off ~34000x at p = 1e-6. ChiSquared inherits it for fractional degrees of freedom (dof < ~1.2).

Cause

Both entries delegate to SpecialFunctions.GammaLowerRegularizedInv, a Cephes igami port that resolves the root by linear-space bisection over [0, big]. For small shape the true quantile is far below the resolver's ~1e-15 floor, and the Wilson–Hilferty initial guess underflows negative, so the routine floors. There are two coupled defects sharing one root cause — a small value rounded to zero:

1. GammaLowerRegularized rounded any x <= ~2.22e-15 up to zero (x.AlmostEqual(0.0)) and returned P = 0, even though P(0.1, 6e-61) = 1e-6 is a perfectly representable double. Only exact zero should short-circuit; genuine underflow is already caught by the existing ax guard.
2. GammaLowerRegularizedInv now solves P(a,x) = p directly for small quantiles, seeding from the leading series inversion x ~ (p·Γ(a+1))^(1/a) and refining with Newton. GammaLowerRegularized keeps full relative precision as x -> 0, so the residual avoids the 1 - P cancellation the old path relied on. Moderate/large shape still take the original path unchanged.

Verification

• mpmath oracle (dps 30 and 35 in agreement): small-shape quantiles now match to ~1e-14 across a grid of shape {0.1 … 20} x probability {1e-9 … 0.999999}.
• Oracle-free self-consistency: CDF(InvCDF(p)) == p holds to ~1e-15, and the inverse is monotone in p.
• No regression: moderate and large shape are byte-for-byte the legacy path; full SpecialFunctions + Distributions suites pass (4719 tests).

New tests cover the forward function (GammaLowerRegularizedSmallArgument), the inverse (GammaLowerRegularizedInvSmallShapeLowerTail, including trigger-boundary cases where the Newton refinement is load-bearing), self-consistency/monotonicity, and both distributions.