From 8ff18847006efe1c66127808a890d8db46ee21dc Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Tue, 28 Jul 2026 18:17:23 +0200 Subject: [PATCH] Fix NaN in gammpapprox lower branch for large a For a >= ASWITCH (100) the Gauss-Legendre quadrature path took the log of the integration width (xu - x), which is negative in the x <= a-1 branch, so (xu - x).ln() evaluated to NaN. That poisoned gammp/gammq (chi-square / gamma / Poisson CDF returned NaN over the whole 0.2a <= x <= a-1 region) and made invgammp panic for lower-tail p < 0.5 when the NaN reached gammp's assert. The log-space rewrite that introduced this kept only the x > a1 case of the Numerical Recipes sign convention and dropped the lower-branch fold. Take the log of |xu - x| and restore the per-branch sign on return: in the lower branch the signed integral is negative, so P = ans and Q = 1 - ans; the upper branch is unchanged (P = 1 - ans, Q = ans). With gammpapprox correct, invgammp converges for large a on both tails without any extra guard. Values verified against mpmath; the small-a gser/gcf path is untouched. --- src/gamma.rs | 19 ++++-- tests/gammpapprox_large_a_test.rs | 110 ++++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+), 5 deletions(-) create mode 100644 tests/gammpapprox_large_a_test.rs diff --git a/src/gamma.rs b/src/gamma.rs index 62a4dc9..efb67a9 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -397,9 +397,12 @@ fn gammpapprox(a: f64, x: f64, psig: IncGamma) -> f64 { } } - // Compute final answer with overflow/underflow protection + // Compute final answer with overflow/underflow protection. + // In the x <= a1 branch xu < x, so the integration width (xu - x) is + // negative; take the log of its magnitude and fold the sign into the + // return below, matching the signed `ans` in Numerical Recipes. let log_scale = a1 * (lna1 - 1f64) - gln; - let log_ans = log_max + sum.ln() + (xu - x).ln() + log_scale; + let log_ans = log_max + sum.ln() + (xu - x).abs().ln() + log_scale; let ans = if log_ans > 700.0 { f64::INFINITY @@ -409,13 +412,19 @@ fn gammpapprox(a: f64, x: f64, psig: IncGamma) -> f64 { log_ans.exp() }; - // Clamp results to [0, 1] and handle the P vs Q case + // ans is the magnitude |integral|. The signed NR value is negative when + // xu < x (the lower branch), so P = ans and Q = 1 - ans there; for the + // upper branch (xu > x) the value is positive, giving P = 1 - ans, Q = ans. + let lower = x <= a1; match psig { IncGamma::P => { - let result = 1f64 - ans; + let result = if lower { ans } else { 1f64 - ans }; + result.clamp(0.0, 1.0) + } + IncGamma::Q => { + let result = if lower { 1f64 - ans } else { ans }; result.clamp(0.0, 1.0) } - IncGamma::Q => ans.clamp(0.0, 1.0), } } diff --git a/tests/gammpapprox_large_a_test.rs b/tests/gammpapprox_large_a_test.rs new file mode 100644 index 0000000..326dceb --- /dev/null +++ b/tests/gammpapprox_large_a_test.rs @@ -0,0 +1,110 @@ +// Regression tests for the large-a (a >= ASWITCH = 100) Gauss-Legendre +// quadrature branch of gammp/gammq/invgammp. +// +// Reference values are mpmath (dps 30, stable at dps 35): +// P = mp.gammainc(a, 0, x, regularized=True) +// Q = mp.gammainc(a, x, mp.inf, regularized=True) +// The lower branch x <= a-1 previously returned NaN, so these tables cover +// x spanning {x << a-1, x ~ a-1, x = a, x > a} for several a >= 100. + +use approx::assert_relative_eq; +use puruspe::{gammp, gammq, invgammp}; + +// (a, x, P, Q); lower branch (x <= a-1) is the region that used to be NaN. +const LARGE_A_TABLE: [(f64, f64, f64, f64); 10] = [ + (100.0, 50.0, 3.20006532459e-10, 0.99999999968), // deep lower tail + (100.0, 90.0, 0.158220989186, 0.841779010814), // mid lower + (100.0, 99.0, 0.473304330399, 0.526695669601), // boundary x = a-1 + (100.0, 99.1, 0.477310410184, 0.522689589816), // just into upper branch + (100.0, 120.0, 0.972136260109, 0.0278637398905), // upper + (150.0, 120.0, 0.00456344130415, 0.995436558696), // lower + (150.0, 150.0, 0.510858229749, 0.489141770251), // x = a (upper) + (200.0, 100.0, 9.34315007299e-19, 1.0), // extreme lower tail + (500.0, 300.0, 3.80742341629e-26, 1.0), // extreme lower tail + (500.0, 600.0, 0.999987744058, 1.22559423306e-5), // upper +]; + +#[test] +fn test_gammp_gammq_large_a_quadrature() { + for &(a, x, p_ref, q_ref) in LARGE_A_TABLE.iter() { + let p = gammp(a, x); + let q = gammq(a, x); + eprintln!( + "a={a} x={x}: P computed={p:.12e} ref={p_ref:.12e} | Q computed={q:.12e} ref={q_ref:.12e}" + ); + assert!(p.is_finite() && q.is_finite(), "P/Q must be finite, got P={} Q={}", p, q); + assert!((0.0..=1.0).contains(&p), "P out of [0,1]: {}", p); + assert!((0.0..=1.0).contains(&q), "Q out of [0,1]: {}", q); + // Mid-range values (0.158, 0.00456, 0.472, 0.972) are asserted with + // full precision so a vacuous near-0/near-1 match cannot pass. + assert_relative_eq!(p, p_ref, epsilon = 1e-14, max_relative = 1e-9); + assert_relative_eq!(q, q_ref, epsilon = 1e-14, max_relative = 1e-9); + // Complement property across the whole surface. + assert_relative_eq!(p + q, 1.0, epsilon = 1e-12); + } +} + +#[test] +fn test_gammp_monotonic_increasing_large_a() { + // P(a, .) is a CDF: strictly increasing in x, spanning the former NaN region. + for &a in &[100.0_f64, 150.0, 500.0] { + let mut prev = -1.0; + let mut x = 0.2 * a; // start of the quadrature branch + while x <= a + 6.0 * a.sqrt() { + let p = gammp(a, x); + assert!(p.is_finite(), "gammp({},{}) not finite", a, x); + assert!(p > prev, "gammp not monotonic at a={} x={}: {} <= {}", a, x, p, prev); + prev = p; + x += 2.0; + } + } +} + +#[test] +fn test_small_a_path_unchanged() { + // a < ASWITCH uses gser/gcf, never the quadrature branch; these must be + // byte-identical to pre-fix behaviour (mpmath references). + let cases = [ + (5.0_f64, 3.0_f64, 0.18473675547942, 0.81526324452058), + (50.0, 40.0, 0.07033506665593, 0.92966493334407), + (5.0, 10.0, 0.97074731192252, 0.02925268807748), + ]; + for &(a, x, p_ref, q_ref) in cases.iter() { + let p = gammp(a, x); + let q = gammq(a, x); + eprintln!("small-a a={a} x={x}: P={p:.14} ref={p_ref:.14}"); + assert_relative_eq!(p, p_ref, epsilon = 1e-13, max_relative = 1e-11); + assert_relative_eq!(q, q_ref, epsilon = 1e-13, max_relative = 1e-11); + } +} + +#[test] +fn test_invgammp_round_trip_large_a() { + // invgammp previously panicked for a >= 100, lower-tail p < 0.5, because + // gammp fed a NaN into its assert. Round-trip both tails with no panic. + for &a in &[100.0_f64, 150.0, 500.0] { + for &p in &[1e-6, 1e-3, 0.1, 0.45, 0.5, 0.55, 0.9, 0.999] { + let x = invgammp(p, a); + assert!(x.is_finite() && x > 0.0, "invgammp({},{}) = {}", p, a, x); + let back = gammp(a, x); + eprintln!("a={a} p={p}: x={x:.10} gammp(a,x)={back:.14}"); + assert_relative_eq!(back, p, epsilon = 1e-10, max_relative = 1e-8); + } + } +} + +#[test] +fn test_invgammp_matches_reference_large_a() { + // Direct check against mpmath findroot values (not just self-consistent). + let cases = [ + (100.0_f64, 0.1_f64, 87.4176365), + (150.0, 0.001, 114.9817176), + (500.0, 0.45, 496.8631868), + (100.0, 0.9, 113.0105239), + ]; + for &(a, p, x_ref) in cases.iter() { + let x = invgammp(p, a); + eprintln!("invgammp a={a} p={p}: x={x:.10} ref={x_ref:.7}"); + assert_relative_eq!(x, x_ref, epsilon = 1e-6, max_relative = 1e-7); + } +}