From bd26b8e3bf6c34b03b9f4f38c63644c1fc2759fb Mon Sep 17 00:00:00 2001 From: youdie006 Date: Mon, 28 Sep 2026 22:31:20 +0900 Subject: [PATCH 1/3] fix: return 0 or 1 from Chi and Gamma cdf/sf when the argument under/overflows Chi passes x * x / 2 and Gamma passes x * rate to gamma_lr/gamma_ur, which panic when their argument is 0 or infinite. The existing guards only checked x itself, so a finite x whose product overflowed or underflowed panicked. ChiSquared and Erlang go through Gamma. --- src/distribution/chi.rs | 28 ++++++++++++++++++++-------- src/distribution/gamma.rs | 16 ++++++++++++++-- 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index 208047db..ead08266 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -119,12 +119,13 @@ impl ContinuousCDF for Chi { /// where `k` is the degrees of freedom and `P` is /// the regularized lower incomplete Gamma function fn cdf(&self, x: f64) -> f64 { - if x == f64::INFINITY { - 1.0 - } else if x <= 0.0 { + let u = x * x / 2.0; + if x <= 0.0 || u == 0.0 { 0.0 + } else if u == f64::INFINITY { + 1.0 } else { - gamma::gamma_lr(self.freedom() as f64 / 2.0, x * x / 2.0) + gamma::gamma_lr(self.freedom() as f64 / 2.0, u) } } @@ -140,12 +141,13 @@ impl ContinuousCDF for Chi { /// where `k` is the degrees of freedom and `P` is /// the regularized upper incomplete Gamma function fn sf(&self, x: f64) -> f64 { - if x == f64::INFINITY { - 0.0 - } else if x <= 0.0 { + let u = x * x / 2.0; + if x <= 0.0 || u == 0.0 { 1.0 + } else if u == f64::INFINITY { + 0.0 } else { - gamma::gamma_ur(self.freedom() as f64 / 2.0, x * x / 2.0) + gamma::gamma_ur(self.freedom() as f64 / 2.0, u) } } @@ -503,6 +505,11 @@ mod tests { test_exact(2, 1.0, cdf(f64::INFINITY)); test_exact(2, 0.0, cdf(0.0)); test_exact(2, 1.0, cdf(f64::INFINITY)); + // x * x / 2 underflows to 0 or overflows to infinity + test_exact(1, 0.0, cdf(1e-200)); + test_exact(3, 0.0, cdf(1e-200)); + test_exact(1, 1.0, cdf(1e160)); + test_exact(3, 1.0, cdf(f64::MAX)); } #[test] @@ -520,18 +527,23 @@ mod tests { test_exact(2, 0.0, sf(f64::INFINITY)); test_exact(2, 1.0, sf(0.0)); test_exact(2, 0.0, sf(f64::INFINITY)); + test_exact(1, 1.0, sf(1e-200)); + test_exact(1, 0.0, sf(1e160)); + test_exact(3, 0.0, sf(f64::MAX)); } #[test] fn test_neg_cdf() { let cdf = |arg: f64| move |x: Chi| x.cdf(arg); test_exact(1, 0.0, cdf(-1.0)); + test_exact(1, 0.0, cdf(f64::NEG_INFINITY)); } #[test] fn test_neg_sf() { let sf = |arg: f64| move |x: Chi| x.sf(arg); test_exact(1, 1.0, sf(-1.0)); + test_exact(1, 1.0, sf(f64::NEG_INFINITY)); } #[test] diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs index 6df26e38..a4d2da93 100644 --- a/src/distribution/gamma.rs +++ b/src/distribution/gamma.rs @@ -152,8 +152,10 @@ impl ContinuousCDF for Gamma { 1.0 } else if self.rate.is_infinite() { 0.0 - } else if x.is_infinite() { + } else if x.is_infinite() || (x * self.rate).is_infinite() { 1.0 + } else if x * self.rate == 0.0 { + 0.0 } else { gamma::gamma_lr(self.shape, x * self.rate) } @@ -177,8 +179,10 @@ impl ContinuousCDF for Gamma { 0.0 } else if self.rate.is_infinite() { 1.0 - } else if x.is_infinite() { + } else if x.is_infinite() || (x * self.rate).is_infinite() { 0.0 + } else if x * self.rate == 0.0 { + 1.0 } else { gamma::gamma_ur(self.shape, x * self.rate) } @@ -745,6 +749,14 @@ mod tests { test_relative(1.0, 0.1, 0.0, |x| x.cdf(0.0)); } + #[test] + fn test_cdf_sf_when_x_times_rate_overflows_or_underflows() { + test_exact(2.5, 1.5, 1.0, |g| g.cdf(f64::MAX)); + test_exact(2.5, 1.5, 0.0, |g| g.sf(f64::MAX)); + test_exact(1.0, 0.5, 0.0, |g| g.cdf(5e-324)); + test_exact(1.0, 0.5, 1.0, |g| g.sf(5e-324)); + } + #[test] fn test_infinite_rate_special_point_is_exact() { let dist = create_ok(10.0, f64::INFINITY); From 5ff489dcd59ed2017996147f44435ecc15a19815 Mon Sep 17 00:00:00 2001 From: youdie006 Date: Wed, 28 Oct 2026 16:32:34 +0900 Subject: [PATCH 2/3] fix: accept x = 0 in gamma_lr/gamma_ur and keep Chi::cdf accurate for tiny x gamma_lr(a, 0) and gamma_ur(a, 0) now return their limits, 0 and 1, instead of XInvalid; an infinite x is still rejected. The Gamma distribution no longer needs its own guard for x * rate underflowing to 0. Chi::cdf passes x * x / 2, which loses precision below f64::MIN_POSITIVE (Chi::new(1).cdf(1e-200) gave 0, the true value is about 7.98e-201). There it now uses the leading series term (x^2 / 2)^(k/2) / Gamma(k/2 + 1) in log space. --- src/distribution/chi.rs | 17 ++++++++++++----- src/distribution/gamma.rs | 4 ---- src/function/gamma.rs | 32 ++++++++++++++++++++------------ 3 files changed, 32 insertions(+), 21 deletions(-) diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index ead08266..c087b42d 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -119,13 +119,17 @@ impl ContinuousCDF for Chi { /// where `k` is the degrees of freedom and `P` is /// the regularized lower incomplete Gamma function fn cdf(&self, x: f64) -> f64 { + let k = self.freedom() as f64; let u = x * x / 2.0; - if x <= 0.0 || u == 0.0 { + if x <= 0.0 { 0.0 } else if u == f64::INFINITY { 1.0 + } else if u < f64::MIN_POSITIVE { + // x * x / 2 has lost precision; use the leading term (x^2 / 2)^(k/2) / Gamma(k/2 + 1) + (k * x.ln() - k / 2.0 * f64_consts::LN_2 - gamma::ln_gamma(k / 2.0 + 1.0)).exp() } else { - gamma::gamma_lr(self.freedom() as f64 / 2.0, u) + gamma::gamma_lr(k / 2.0, u) } } @@ -142,7 +146,7 @@ impl ContinuousCDF for Chi { /// the regularized upper incomplete Gamma function fn sf(&self, x: f64) -> f64 { let u = x * x / 2.0; - if x <= 0.0 || u == 0.0 { + if x <= 0.0 { 1.0 } else if u == f64::INFINITY { 0.0 @@ -505,8 +509,11 @@ mod tests { test_exact(2, 1.0, cdf(f64::INFINITY)); test_exact(2, 0.0, cdf(0.0)); test_exact(2, 1.0, cdf(f64::INFINITY)); - // x * x / 2 underflows to 0 or overflows to infinity - test_exact(1, 0.0, cdf(1e-200)); + // x * x / 2 underflows or overflows; for k = 1, cdf(x) = erf(x / sqrt(2)) ~ x * sqrt(2 / pi) + for x in [1e-160, 1e-200] { + let expected = x * (2.0 / f64_consts::PI).sqrt(); + crate::prec::assert_relative_eq!(create_ok(1).cdf(x), expected, epsilon = 0.0, max_relative = 1e-12); + } test_exact(3, 0.0, cdf(1e-200)); test_exact(1, 1.0, cdf(1e160)); test_exact(3, 1.0, cdf(f64::MAX)); diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs index a4d2da93..6b92cc20 100644 --- a/src/distribution/gamma.rs +++ b/src/distribution/gamma.rs @@ -154,8 +154,6 @@ impl ContinuousCDF for Gamma { 0.0 } else if x.is_infinite() || (x * self.rate).is_infinite() { 1.0 - } else if x * self.rate == 0.0 { - 0.0 } else { gamma::gamma_lr(self.shape, x * self.rate) } @@ -181,8 +179,6 @@ impl ContinuousCDF for Gamma { 1.0 } else if x.is_infinite() || (x * self.rate).is_infinite() { 0.0 - } else if x * self.rate == 0.0 { - 1.0 } else { gamma::gamma_ur(self.shape, x * self.rate) } diff --git a/src/function/gamma.rs b/src/function/gamma.rs index 51d5fed1..ff314163 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -14,7 +14,7 @@ pub enum GammaFuncError { /// `a` is infinite, zero or less than zero. AInvalid, - /// `x` is infinite, zero or less than zero. + /// `x` is infinite or less than zero. XInvalid, } @@ -22,7 +22,7 @@ impl core::fmt::Display for GammaFuncError { fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result { match self { GammaFuncError::AInvalid => write!(f, "a is infinite, zero or less than zero"), - GammaFuncError::XInvalid => write!(f, "x is infinite, zero or less than zero"), + GammaFuncError::XInvalid => write!(f, "x is infinite or less than zero"), } } } @@ -113,7 +113,7 @@ pub fn gamma(x: f64) -> f64 { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_ui(a: f64, x: f64) -> f64 { checked_gamma_ui(a, x).unwrap() } @@ -125,7 +125,7 @@ pub fn gamma_ui(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_ui(a: f64, x: f64) -> Result { checked_gamma_ur(a, x).map(|x| x * gamma(a)) } @@ -138,7 +138,7 @@ pub fn checked_gamma_ui(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_li(a: f64, x: f64) -> f64 { checked_gamma_li(a, x).unwrap() } @@ -151,7 +151,7 @@ pub fn gamma_li(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_li(a: f64, x: f64) -> Result { checked_gamma_lr(a, x).map(|x| x * gamma(a)) } @@ -167,7 +167,7 @@ pub fn checked_gamma_li(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_ur(a: f64, x: f64) -> f64 { checked_gamma_ur(a, x).unwrap() } @@ -183,7 +183,7 @@ pub fn gamma_ur(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_ur(a: f64, x: f64) -> Result { if a.is_nan() || x.is_nan() { return Ok(f64::NAN); @@ -191,7 +191,7 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result { if a <= 0.0 || a == f64::INFINITY { return Err(GammaFuncError::AInvalid); } - if x <= 0.0 || x == f64::INFINITY { + if x < 0.0 || x == f64::INFINITY { return Err(GammaFuncError::XInvalid); } @@ -261,7 +261,7 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_lr(a: f64, x: f64) -> f64 { checked_gamma_lr(a, x).unwrap() } @@ -277,7 +277,7 @@ pub fn gamma_lr(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_lr(a: f64, x: f64) -> Result { if a.is_nan() || x.is_nan() { return Ok(f64::NAN); @@ -285,7 +285,7 @@ pub fn checked_gamma_lr(a: f64, x: f64) -> Result { if a <= 0.0 || a == f64::INFINITY { return Err(GammaFuncError::AInvalid); } - if x <= 0.0 || x == f64::INFINITY { + if x < 0.0 || x == f64::INFINITY { return Err(GammaFuncError::XInvalid); } @@ -915,6 +915,14 @@ mod tests { assert!(super::checked_gamma_lr(1.0, f64::INFINITY).is_err()); } + #[test] + fn test_gamma_lr_ur_at_x_zero() { + assert_eq!(gamma_lr(2.5, 0.0), 0.0); + assert_eq!(gamma_ur(2.5, 0.0), 1.0); + assert_eq!(super::checked_gamma_lr(0.5, 0.0), Ok(0.0)); + assert_eq!(super::checked_gamma_ur(0.5, 0.0), Ok(1.0)); + } + #[test] fn test_gamma_li() { assert!(gamma_li(f64::NAN, f64::NAN).is_nan()); From 73f14f41987ea026f1d01b1391ea4ebb615bbf34 Mon Sep 17 00:00:00 2001 From: youdie006 Date: Thu, 29 Oct 2026 13:44:12 +0900 Subject: [PATCH 3/3] fix: return 0 from gamma_li at x = 0 when gamma(a) overflows --- src/function/gamma.rs | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/function/gamma.rs b/src/function/gamma.rs index ff314163..eaf793ad 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -153,7 +153,8 @@ pub fn gamma_li(a: f64, x: f64) -> f64 { /// /// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_li(a: f64, x: f64) -> Result { - checked_gamma_lr(a, x).map(|x| x * gamma(a)) + // gamma(a) can overflow, and 0 * inf is NaN + checked_gamma_lr(a, x).map(|p| if x == 0.0 { 0.0 } else { p * gamma(a) }) } /// Computes the upper incomplete regularized gamma function @@ -916,11 +917,13 @@ mod tests { } #[test] - fn test_gamma_lr_ur_at_x_zero() { + fn test_incomplete_gamma_at_x_zero() { assert_eq!(gamma_lr(2.5, 0.0), 0.0); assert_eq!(gamma_ur(2.5, 0.0), 1.0); assert_eq!(super::checked_gamma_lr(0.5, 0.0), Ok(0.0)); assert_eq!(super::checked_gamma_ur(0.5, 0.0), Ok(1.0)); + assert_eq!(gamma_li(172.0, 0.0), 0.0); + assert_eq!(super::checked_gamma_li(172.0, 0.0), Ok(0.0)); } #[test]