From fa2470e82c4fe74dbc7df1a3fd9dfb9585bd8843 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 17:30:13 +0800 Subject: [PATCH 1/5] Improve numerical stability of `ln_gamma_approx` using Euler's reflection formula --- src/special/lanczos.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index 2e70b01..7142f59 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -23,6 +23,10 @@ const LG5N7: [f64; 7] = [ ]; pub fn ln_gamma_approx(z: f64) -> f64 { + if z < 0.5 { + return PI.ln() - (PI * z).sin().abs().ln() - ln_gamma_approx(1.0 - z); + } + let z = z - 1f64; let base = z + G + 0.5; let mut s = 0f64; From 931e421906e81542340556c8035bc5ed6be42fa7 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 17:46:19 +0800 Subject: [PATCH 2/5] Handle poles --- src/special/lanczos.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index 7142f59..92548ae 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -23,6 +23,10 @@ const LG5N7: [f64; 7] = [ ]; pub fn ln_gamma_approx(z: f64) -> f64 { + if z <= 0.0 && z.fract() == 0.0 { + return f64::INFINITY; + } + if z < 0.5 { return PI.ln() - (PI * z).sin().abs().ln() - ln_gamma_approx(1.0 - z); } From 59f985d8b3fa5f4af96311fb0b6c1288cb05c996 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 17:46:30 +0800 Subject: [PATCH 3/5] Handle poles --- src/special/lanczos.rs | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index 92548ae..e3ea028 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -26,7 +26,7 @@ pub fn ln_gamma_approx(z: f64) -> f64 { if z <= 0.0 && z.fract() == 0.0 { return f64::INFINITY; } - + if z < 0.5 { return PI.ln() - (PI * z).sin().abs().ln() - ln_gamma_approx(1.0 - z); } @@ -42,18 +42,19 @@ pub fn ln_gamma_approx(z: f64) -> f64 { } pub fn gamma_approx(z: f64) -> f64 { - if z > 1f64 { - let z_int = z as usize; - if z - (z_int as f64) == 0f64 { - return factorial(z_int - 1) as f64; + if z <= 0.0 && z.fract() == 0.0 { + if z == 0.0 { + return f64::INFINITY; + } else { + return f64::NAN; } } if z < 0.5 { - PI / ((PI * z).sin() * gamma_approx(1f64 - z)) - } else { - ln_gamma_approx(z).exp() + return PI / ((PI * z).sin() * gamma_approx(1f64 - z)); } + + ln_gamma_approx(z).exp() } /// Lanczos Approximation Coefficient From ef7068bba76f551461eb7461ffc9d52a17f49d02 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 18:03:11 +0800 Subject: [PATCH 4/5] Reimplement factorial fast-path --- src/special/lanczos.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index e3ea028..b4f859a 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -50,6 +50,16 @@ pub fn gamma_approx(z: f64) -> f64 { } } + if z > 1.0 && z.fract() == 0.0 { + let mut result = 1.0; + let n = (z - 1.0) as u64; + for i in 2..=n { + result *= i as f64; + } + + return result; + } + if z < 0.5 { return PI / ((PI * z).sin() * gamma_approx(1f64 - z)); } From 1bd8c449bc37a767e10b7808abef34b39f3a5f8d Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 18:18:11 +0800 Subject: [PATCH 5/5] Add tests --- src/special/lanczos.rs | 4 +- tests/special.rs | 101 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 2 deletions(-) diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index b4f859a..067380a 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -1,6 +1,6 @@ //! Lanczos approximation Coefficient generator -use crate::statistics::ops::{double_factorial, factorial, C}; +use crate::statistics::ops::{double_factorial, C}; use crate::structure::matrix::Matrix; use crate::traits::matrix::MatrixTrait; use crate::traits::pointer::{Oxide, RedoxCommon}; @@ -50,7 +50,7 @@ pub fn gamma_approx(z: f64) -> f64 { } } - if z > 1.0 && z.fract() == 0.0 { + if z > 0.0 && z.fract() == 0.0 { let mut result = 1.0; let n = (z - 1.0) as u64; for i in 2..=n { diff --git a/tests/special.rs b/tests/special.rs index acec5c6..610c9ba 100644 --- a/tests/special.rs +++ b/tests/special.rs @@ -1,7 +1,108 @@ use peroxide::fuga::{LambertWAccuracyMode::*, *}; +use std::f64::consts::PI; #[test] fn lambert_w_test() { assert_eq!(lambert_w0(1.0, Precise), 0.567143290409784); assert!(nearly_eq(lambert_w0(1.0, Simple), 0.567143290409784)); } + +#[test] +fn test_gamma_poles_and_undefined() { + // Gamma(0) approaches infinity + assert!(gamma(0.0).is_infinite()); + assert!(gamma(0.0).is_sign_positive()); + + // Gamma for negative integers is mathematically undefined (diverges) + assert!(gamma(-1.0).is_nan()); + assert!(gamma(-2.0).is_nan()); + assert!(gamma(-10.0).is_nan()); + + // Log-Gamma goes to positive infinity for all poles + assert!(ln_gamma(0.0).is_infinite()); + assert!(ln_gamma(-1.0).is_infinite()); + assert!(ln_gamma(-10.0).is_infinite()); +} + +#[test] +fn test_gamma_integer_fast_path() { + // Standard small factorials: Gamma(n) = (n-1)! + assert_eq!(gamma(1.0), 1.0); // 0! + assert_eq!(gamma(2.0), 1.0); // 1! + assert_eq!(gamma(4.0), 6.0); // 3! + assert_eq!(gamma(5.0), 24.0); // 4! + assert_eq!(gamma(10.0), 362_880.0); // 9! + + // Wolfram Alpha high-precision check (21!) + // f64 can exactly represent this without precision loss + assert_eq!(gamma(22.0), 51_090_942_171_709_440_000.0); + + // Maximum limit of f64 float representation (~171.6) + // Ensure it doesn't panic on overflow, but correctly yields Infinity + assert!(gamma(172.0).is_infinite()); +} + +#[test] +fn test_gamma_positive_floats() { + let sqrt_pi = PI.sqrt(); + + // Gamma(0.5) = sqrt(PI) + assert!(nearly_eq(gamma(0.5), sqrt_pi)); + + // Gamma(1.5) = 0.5 * sqrt(PI) + assert!(nearly_eq(gamma(1.5), 0.5 * sqrt_pi)); + + // Gamma(2.5) = 1.329340388179... + assert!(nearly_eq(gamma(2.5), 0.75 * sqrt_pi)); +} + +#[test] +fn test_gamma_negative_floats_reflection() { + let sqrt_pi = PI.sqrt(); + + // Gamma(-0.5) = -2 * sqrt(PI) + // This validates that .abs() is NOT used on the sine in gamma_approx + assert!(nearly_eq(gamma(-0.5), -2.0 * sqrt_pi)); + assert!(gamma(-0.5).is_sign_negative()); + + // Gamma(-1.5) = (4/3) * sqrt(PI) + assert!(nearly_eq(gamma(-1.5), (4.0 / 3.0) * sqrt_pi)); + assert!(gamma(-1.5).is_sign_positive()); + + // Gamma(-2.5) = -(8/15) * sqrt(PI) + assert!(nearly_eq(gamma(-2.5), -(8.0 / 15.0) * sqrt_pi)); + assert!(gamma(-2.5).is_sign_negative()); +} + +#[test] +fn test_ln_gamma_consistency() { + // ln_gamma(x) should equal ln(|Gamma(x)|) across the board + let test_values = vec![0.5, 1.5, 2.5, 10.5]; + + for &val in &test_values { + let expected = gamma(val).ln(); + let actual = ln_gamma(val); + assert!( + nearly_eq(expected, actual), + "Failed at positive float: val={}, expected={}, actual={}", + val, + expected, + actual + ); + } + + // Test Negative Floats to ensure `.abs()` prevents NaN + let negative_test_values = vec![-0.5, -1.5, -2.5, -10.5]; + + for &val in &negative_test_values { + let expected = gamma(val).abs().ln(); + let actual = ln_gamma(val); + assert!( + nearly_eq(expected, actual), + "Failed at negative float: val={}, expected={}, actual={}", + val, + expected, + actual + ); + } +}