Fix NaN in gammpapprox lower branch for large a

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),
     }
 }
 

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);
+    }
+}