Skip to content

fix: return 0 or 1 from Chi and Gamma cdf/sf when the argument under/overflows #475

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
youdie006 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

fix: return 0 or 1 from Chi and Gamma cdf/sf when the argument under/overflows #475

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bd26b8e
fix: return 0 or 1 from Chi and Gamma cdf/sf when the argument under/…
youdie006 Sep 28, 2026
5ff489d
fix: accept x = 0 in gamma_lr/gamma_ur and keep Chi::cdf accurate for…
youdie006 Oct 28, 2026
73f14f4
fix: return 0 from gamma_li at x = 0 when gamma(a) overflows
youdie006 Oct 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 27 additions & 8 deletions src/distribution/chi.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,12 +119,17 @@ impl ContinuousCDF<f64, f64> 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 k = self.freedom() as f64;
let u = x * x / 2.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, x * x / 2.0)
gamma::gamma_lr(k / 2.0, u)
}
}

Comment thread
coderabbitai[bot] marked this conversation as resolved.
Resolved
Expand All @@ -140,12 +145,13 @@ impl ContinuousCDF<f64, f64> 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 {
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)
}
}

Expand Down Expand Up @@ -503,6 +509,14 @@ 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 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));
}

#[test]
Expand All @@ -520,18 +534,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]
Expand Down
12 changes: 10 additions & 2 deletions src/distribution/gamma.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ impl ContinuousCDF<f64, f64> 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 {
gamma::gamma_lr(self.shape, x * self.rate)
Expand All @@ -177,7 +177,7 @@ impl ContinuousCDF<f64, f64> 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 {
gamma::gamma_ur(self.shape, x * self.rate)
Expand Down Expand Up @@ -745,6 +745,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);
Expand Down
37 changes: 24 additions & 13 deletions src/function/gamma.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,15 @@ 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,
}

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"),
}
}
}
Expand Down Expand Up @@ -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()
}
Expand All @@ -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<f64, GammaFuncError> {
checked_gamma_ur(a, x).map(|x| x * gamma(a))
}
Expand All @@ -138,7 +138,7 @@ pub fn checked_gamma_ui(a: f64, x: f64) -> Result<f64, GammaFuncError> {
///
/// # 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()
}
Expand All @@ -151,9 +151,10 @@ 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<f64, GammaFuncError> {
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
Expand All @@ -167,7 +168,7 @@ pub fn checked_gamma_li(a: f64, x: f64) -> Result<f64, GammaFuncError> {
///
/// # 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()
}
Expand All @@ -183,15 +184,15 @@ 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<f64, GammaFuncError> {
if a.is_nan() || x.is_nan() {
return Ok(f64::NAN);
}
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);
}

Expand Down Expand Up @@ -261,7 +262,7 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result<f64, GammaFuncError> {
///
/// # 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()
}
Expand All @@ -277,15 +278,15 @@ 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<f64, GammaFuncError> {
if a.is_nan() || x.is_nan() {
return Ok(f64::NAN);
}
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);
}

Expand Down Expand Up @@ -915,6 +916,16 @@ mod tests {
assert!(super::checked_gamma_lr(1.0, f64::INFINITY).is_err());
}

#[test]
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]
fn test_gamma_li() {
assert!(gamma_li(f64::NAN, f64::NAN).is_nan());
Expand Down
Loading