From c34f7bd6eae7b207dca82fcb3155fa3fb22ecde8 Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Fri, 18 Apr 2025 20:11:13 +0900 Subject: [PATCH 1/5] fix --- src/gamma.rs | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 2 deletions(-) diff --git a/src/gamma.rs b/src/gamma.rs index 7dd6f98..12ef13c 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -173,7 +173,7 @@ pub fn tgamma(x: f64) -> Result { q - absx }; let z = z * LANCZOS_G / y; - let r = if x < 0.0 { + let mut r = if x < 0.0 { let mut r = -PI / m_sinpi(absx) / absx * y.exp() / lanczos_sum(absx); r -= z * r; if absx < 140.0 { @@ -196,6 +196,7 @@ pub fn tgamma(x: f64) -> Result { } r }; + if r.is_infinite() { return Err((f64::INFINITY, Error::ERANGE).1); } else { @@ -241,7 +242,14 @@ pub fn lgamma(x: f64) -> Result { if x < 0.0 { // Use reflection formula to get value for negative x - r = LOG_PI - m_sinpi(absx).abs().ln() - absx.ln() - r; + + // Calculate the sin(pi * x) value using m_sinpi + let sinpi_val = m_sinpi(absx); + + // In CPython, the expression is: + // r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r; + // We'll match this order exactly + r = LOG_PI - sinpi_val.abs().ln() - absx.ln() - r; } if r.is_infinite() { return Err(Error::ERANGE); @@ -283,6 +291,60 @@ mod tests { } } + #[test] + fn test_specific_lgamma_value() { + let x = -3.8510064710745118; + let rs_lgamma = lgamma(x).unwrap(); + + pyo3::prepare_freethreaded_python(); + Python::with_gil(|py| { + let math = PyModule::import(py, "math").unwrap(); + let py_lgamma = math + .getattr("lgamma") + .unwrap() + .call1((x,)) + .unwrap() + .extract::() + .unwrap(); + + println!("x = {}", x); + println!("Python lgamma = {} ({:x})", py_lgamma, unsafe { + std::mem::transmute::(py_lgamma) + }); + println!("Rust lgamma = {} ({:x})", rs_lgamma, unsafe { + std::mem::transmute::(rs_lgamma) + }); + + // Print intermediate values + let absx = x.abs(); + let sinpi_val = m_sinpi(absx); + + println!("absx = {}", absx); + println!("m_sinpi = {}", sinpi_val); + + // Compare with Python's sin(pi * x) + let py_code = PyModule::from_code( + py, + c"import math\ndef sinpi(x): return math.sin(math.pi * x)\n", + c"", + c"", + ) + .unwrap(); + let py_sinpi = py_code + .getattr("sinpi") + .unwrap() + .call1((absx,)) + .unwrap() + .extract::() + .unwrap(); + println!("Python sinpi = {}", py_sinpi); + + let py_lgamma_repr = unsafe { std::mem::transmute::(py_lgamma) }; + let rs_lgamma_repr = unsafe { std::mem::transmute::(rs_lgamma) }; + println!("Bit difference: {}", py_lgamma_repr ^ rs_lgamma_repr); + }); + } + proptest! { #[test] fn test_tgamma(x: f64) { From e8ebe6af0cb58c4bbeafc61ae16fee7ef8dc0b1a Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Sun, 20 Apr 2025 14:30:05 +0900 Subject: [PATCH 2/5] working --- check_constants.c | 91 +++++++++++++++++++ src/gamma.rs | 223 +++++++++++++++++++++++++++------------------- 2 files changed, 222 insertions(+), 92 deletions(-) create mode 100644 check_constants.c diff --git a/check_constants.c b/check_constants.c new file mode 100644 index 0000000..08c78bb --- /dev/null +++ b/check_constants.c @@ -0,0 +1,91 @@ +#include +#include +#include // For memcpy +#include // For M_PI if needed, though PI is defined below + +// Function to print the bit representation of a double +void print_double_bits(const char *name, double val) { + uint64_t bits; + // Use memcpy to safely copy the bits, avoiding potential strict-aliasing issues + memcpy(&bits, &val, sizeof(double)); + printf("%s = %.17g (0x%016lx)\n", name, val, bits); +} + +int main() { + // Constants from gamma.rs + const double PI = 3.141592653589793238462643383279502884197; + const double LOG_PI = 1.144729885849400174143427351353058711647; + const double LANCZOS_G = 6.024680040776729583740234375; + const double LANCZOS_G_MINUS_HALF = 5.524680040776729583740234375; + + const int LANCZOS_N = 13; + const double LANCZOS_NUM_COEFFS[LANCZOS_N] = { + 23531376880.410759688572007674451636754734846804940, + 42919803642.649098768957899047001988850926355848959, + 35711959237.355668049440185451547166705960488635843, + 17921034426.037209699919755754458931112671403265390, + 6039542586.3520280050642916443072979210699388420708, + 1439720407.3117216736632230727949123939715485786772, + 248874557.86205415651146038641322942321632125127801, + 31426415.585400194380614231628318205362874684987640, + 2876370.6289353724412254090516208496135991145378768, + 186056.26539522349504029498971604569928220784236328, + 8071.6720023658162106380029022722506138218516325024, + 210.82427775157934587250973392071336271166969580291, + 2.5066282746310002701649081771338373386264310793408, + }; + const double LANCZOS_DEN_COEFFS[LANCZOS_N] = { + 0.0, + 39916800.0, + 120543840.0, + 150917976.0, + 105258076.0, + 45995730.0, + 13339535.0, + 2637558.0, + 357423.0, + 32670.0, + 1925.0, + 66.0, + 1.0, + }; + + const int NGAMMA_INTEGRAL = 23; + const double GAMMA_INTEGRAL[NGAMMA_INTEGRAL] = { + 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, + 362880.0, 3628800.0, 39916800.0, 479001600.0, 6227020800.0, + 87178291200.0, 1307674368000.0, 20922789888000.0, + 355687428096000.0, 6402373705728000.0, 121645100408832000.0, + 2432902008176640000.0, 51090942171709440000.0, + 1124000727777607680000.0, + }; + + printf("--- Single Constants ---\n"); + print_double_bits("PI", PI); // Added PI for completeness + print_double_bits("LOG_PI", LOG_PI); + print_double_bits("LANCZOS_G", LANCZOS_G); + print_double_bits("LANCZOS_G_MINUS_HALF", LANCZOS_G_MINUS_HALF); + + printf("\n--- LANCZOS_NUM_COEFFS ---\n"); + for (int i = 0; i < LANCZOS_N; ++i) { + char name[32]; + snprintf(name, sizeof(name), "LANCZOS_NUM_COEFFS[%d]", i); + print_double_bits(name, LANCZOS_NUM_COEFFS[i]); + } + + printf("\n--- LANCZOS_DEN_COEFFS ---\n"); + for (int i = 0; i < LANCZOS_N; ++i) { + char name[32]; + snprintf(name, sizeof(name), "LANCZOS_DEN_COEFFS[%d]", i); + print_double_bits(name, LANCZOS_DEN_COEFFS[i]); + } + + printf("\n--- GAMMA_INTEGRAL ---\n"); + for (int i = 0; i < NGAMMA_INTEGRAL; ++i) { + char name[32]; + snprintf(name, sizeof(name), "GAMMA_INTEGRAL[%d]", i); + print_double_bits(name, GAMMA_INTEGRAL[i]); + } + + return 0; +} diff --git a/src/gamma.rs b/src/gamma.rs index 12ef13c..97bace6 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -1,40 +1,82 @@ use crate::Error; -use std::f64::consts::PI; +// use std::f64::consts::PI; -const LOG_PI: f64 = 1.144729885849400174143427351353058711647; + +// Import C library functions properly +unsafe extern "C" { + fn exp(x: f64) -> f64; + fn log(x: f64) -> f64; + fn sin(x: f64) -> f64; + fn cos(x: f64) -> f64; + fn pow(x: f64, y: f64) -> f64; + fn floor(x: f64) -> f64; + fn fabs(x: f64) -> f64; + fn round(x: f64) -> f64; + fn fmod(x: f64, y: f64) -> f64; +} + +const PI: f64 = f64::from_bits(0x400921fb54442d18); // = 3.141592653589793238462643383279502884197; +const LOG_PI: f64 = f64::from_bits(0x3ff250d048e7a1bd); // = 1.144729885849400174143427351353058711647; const LANCZOS_N: usize = 13; -const LANCZOS_G: f64 = 6.024680040776729583740234375; -const LANCZOS_G_MINUS_HALF: f64 = 5.524680040776729583740234375; +const LANCZOS_G: f64 = f64::from_bits(0x40181945b9800000); // = 6.024680040776729583740234375; +const LANCZOS_G_MINUS_HALF: f64 = f64::from_bits(0x40161945b9800000); // = 5.524680040776729583740234375; const LANCZOS_NUM_COEFFS: [f64; LANCZOS_N] = [ - 23531376880.410759688572007674451636754734846804940, - 42919803642.649098768957899047001988850926355848959, - 35711959237.355668049440185451547166705960488635843, - 17921034426.037209699919755754458931112671403265390, - 6039542586.3520280050642916443072979210699388420708, - 1439720407.3117216736632230727949123939715485786772, - 248874557.86205415651146038641322942321632125127801, - 31426415.585400194380614231628318205362874684987640, - 2876370.6289353724412254090516208496135991145378768, - 186056.26539522349504029498971604569928220784236328, - 8071.6720023658162106380029022722506138218516325024, - 210.82427775157934587250973392071336271166969580291, - 2.5066282746310002701649081771338373386264310793408, + f64::from_bits(0x4215ea5143c1a49e), // 23531376880.410759 + f64::from_bits(0x4223fc7075f54c57), // 42919803642.649101 + f64::from_bits(0x4220a132818ab61a), // 35711959237.355667 + f64::from_bits(0x4210b0b522e8261a), // 17921034426.037209 + f64::from_bits(0x41f67fc1b3a5a1e8), // 6039542586.3520279 + f64::from_bits(0x41d57418f5d3f33f), // 1439720407.3117216 + f64::from_bits(0x41adab0c7bb95f2a), // 248874557.86205417 + f64::from_bits(0x417df876f95dcc98), // 31426415.585400194 + f64::from_bits(0x4145f1e95080f44c), // 2876370.6289353725 + f64::from_bits(0x4106b6421f8787eb), // 186056.26539522348 + f64::from_bits(0x40bf87ac0858d804), // 8071.6720023658163 + f64::from_bits(0x406a5a607bbc3b52), // 210.82427775157936 + f64::from_bits(0x40040d931ff62705), // 2.5066282746310002 ]; const LANCZOS_DEN_COEFFS: [f64; LANCZOS_N] = [ - 0.0, - 39916800.0, - 120543840.0, - 150917976.0, - 105258076.0, - 45995730.0, - 13339535.0, - 2637558.0, - 357423.0, - 32670.0, - 1925.0, - 66.0, - 1.0, + f64::from_bits(0x0000000000000000), // 0.0 + f64::from_bits(0x418308a800000000), // 39916800.0 + f64::from_bits(0x419cbd6980000000), // 120543840.0 + f64::from_bits(0x41a1fda6b0000000), // 150917976.0 + f64::from_bits(0x4199187170000000), // 105258076.0 + f64::from_bits(0x4185eeb690000000), // 45995730.0 + f64::from_bits(0x41697171e0000000), // 13339535.0 + f64::from_bits(0x41441f7b00000000), // 2637558.0 + f64::from_bits(0x4115d0bc00000000), // 357423.0 + f64::from_bits(0x40dfe78000000000), // 32670.0 + f64::from_bits(0x409e140000000000), // 1925.0 + f64::from_bits(0x4050800000000000), // 66.0 + f64::from_bits(0x3ff0000000000000), // 1.0 +]; + +const NGAMMA_INTEGRAL: usize = 23; +const GAMMA_INTEGRAL: [f64; NGAMMA_INTEGRAL] = [ + f64::from_bits(0x3ff0000000000000), // 1.0 + f64::from_bits(0x3ff0000000000000), // 1.0 + f64::from_bits(0x4000000000000000), // 2.0 + f64::from_bits(0x4018000000000000), // 6.0 + f64::from_bits(0x4038000000000000), // 24.0 + f64::from_bits(0x405e000000000000), // 120.0 + f64::from_bits(0x4086800000000000), // 720.0 + f64::from_bits(0x40b3b00000000000), // 5040.0 + f64::from_bits(0x40e3b00000000000), // 40320.0 + f64::from_bits(0x4116260000000000), // 362880.0 + f64::from_bits(0x414baf8000000000), // 3628800.0 + f64::from_bits(0x418308a800000000), // 39916800.0 + f64::from_bits(0x41bc8cfc00000000), // 479001600.0 + f64::from_bits(0x41f7328cc0000000), // 6227020800.0 + f64::from_bits(0x42344c3b28000000), // 87178291200.0 + f64::from_bits(0x4273077775800000), // 1307674368000.0 + f64::from_bits(0x42b3077775800000), // 20922789888000.0 + f64::from_bits(0x42f437eeecd80000), // 355687428096000.0 + f64::from_bits(0x4336beecca730000), // 6402373705728000.0 + f64::from_bits(0x437b02b930689000), // 1.21645100408832e+17 + f64::from_bits(0x43c0e1b3be415a00), // 2.43290200817664e+18 + f64::from_bits(0x4406283be9b5c620), // 5.109094217170944e+19 + f64::from_bits(0x444e77526159f06c), // 1.1240007277776077e+21 ]; fn lanczos_sum(x: f64) -> f64 { @@ -65,50 +107,23 @@ fn lanczos_sum(x: f64) -> f64 { fn m_sinpi(x: f64) -> f64 { // this function should only ever be called for finite arguments debug_assert!(x.is_finite()); - let y = x.abs() % 2.0; - let n = (2.0 * y).round() as i32; + let y = unsafe { fmod(fabs(x), 2.0) }; + let n = unsafe { round(2.0 * y) } as i32; let r = match n { - 0 => (PI * y).sin(), - 1 => (PI * (y - 0.5)).cos(), + 0 => unsafe { sin(PI * y) }, + 1 => unsafe { cos(PI * (y - 0.5)) }, 2 => { // N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give // -0.0 instead of 0.0 when y == 1.0. - (PI * (1.0 - y)).sin() + unsafe { sin(PI * (1.0 - y)) } } - 3 => -(PI * (y - 1.5)).cos(), - 4 => (PI * (y - 2.0)).sin(), + 3 => unsafe { -cos(PI * (y - 1.5)) }, + 4 => unsafe { sin(PI * (y - 2.0)) }, _ => unreachable!(), }; (1.0f64).copysign(x) * r } -const NGAMMA_INTEGRAL: usize = 23; -const GAMMA_INTEGRAL: [f64; NGAMMA_INTEGRAL] = [ - 1.0, - 1.0, - 2.0, - 6.0, - 24.0, - 120.0, - 720.0, - 5040.0, - 40320.0, - 362880.0, - 3628800.0, - 39916800.0, - 479001600.0, - 6227020800.0, - 87178291200.0, - 1307674368000.0, - 20922789888000.0, - 355687428096000.0, - 6402373705728000.0, - 121645100408832000.0, - 2432902008176640000.0, - 51090942171709440000.0, - 1124000727777607680000.0, -]; - pub fn tgamma(x: f64) -> Result { // special cases if !x.is_finite() { @@ -130,7 +145,7 @@ pub fn tgamma(x: f64) -> Result { return Err((v, Error::EDOM).1); } // integer arguments - if x == x.floor() { + if x == unsafe { floor(x) } { if x < 0.0 { // tgamma(n) = nan, invalid for return Err((f64::NAN, Error::EDOM).1); @@ -139,7 +154,7 @@ pub fn tgamma(x: f64) -> Result { return Ok(GAMMA_INTEGRAL[x as usize - 1]); } } - let absx = x.abs(); + let absx = unsafe { fabs(x) }; // tiny arguments: tgamma(x) ~ 1/x for x near 0 if absx < 1e-20 { let r = 1.0 / x; @@ -173,28 +188,38 @@ pub fn tgamma(x: f64) -> Result { q - absx }; let z = z * LANCZOS_G / y; - let mut r = if x < 0.0 { - let mut r = -PI / m_sinpi(absx) / absx * y.exp() / lanczos_sum(absx); + let r = if x < 0.0 { + // Using C's math functions through libc to match CPython + let term1 = -PI / m_sinpi(absx); + let term2 = term1 / absx; + let exp_y = unsafe { exp(y) }; + let term3 = term2 * exp_y; + let lanczos = lanczos_sum(absx); + let mut r = term3 / lanczos; r -= z * r; + if absx < 140.0 { - r /= y.powf(absx - 0.5); + unsafe { r / pow(y, absx - 0.5) } } else { - let sqrtpow = y.powf(absx / 2.0 - 0.25); + let sqrtpow = unsafe { pow(y, absx / 2.0 - 0.25) }; r /= sqrtpow; r /= sqrtpow; + r } - r } else { - let mut r = lanczos_sum(absx) / y.exp(); + let lanczos = lanczos_sum(absx); + let exp_y = unsafe { exp(y) }; + let mut r = lanczos / exp_y; r += z * r; + if absx < 140.0 { - r *= y.powf(absx - 0.5); + unsafe { r * pow(y, absx - 0.5) } } else { - let sqrtpow = y.powf(absx / 2.0 - 0.25); + let sqrtpow = unsafe { pow(y, absx / 2.0 - 0.25) }; r *= sqrtpow; r *= sqrtpow; + r } - r }; if r.is_infinite() { @@ -217,7 +242,7 @@ pub fn lgamma(x: f64) -> Result { } // integer arguments - if x == x.floor() && x <= 2.0 { + if x == unsafe { floor(x) } && x <= 2.0 { if x <= 0.0 { // lgamma(n) = inf, divide-by-zero for integers n <= 0 return Err(Error::EDOM); @@ -227,30 +252,42 @@ pub fn lgamma(x: f64) -> Result { } } - let absx = x.abs(); + let absx = unsafe { fabs(x) }; // tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x if absx < 1e-20 { - return Ok(-absx.ln()); + return Ok(-unsafe { log(absx) }); } - // Lanczos' formula. We could save a fraction of a ulp in accuracy by - // having a second set of numerator coefficients for lanczos_sum that - // absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g - // subtraction below; it's probably not worth it. - let mut r = lanczos_sum(absx).ln() - LANCZOS_G; - r += (absx - 0.5) * ((absx + LANCZOS_G - 0.5).ln() - 1.0); + // Using C's math functions through libc to match CPython + let lanczos_sum_val = lanczos_sum(absx); + let log_lanczos = unsafe { log(lanczos_sum_val) }; + + // Subtract lanczos_g as a separate step + let mut r = log_lanczos - LANCZOS_G; + + // Calculate (absx - 0.5) term + let factor = absx - 0.5; + + // Calculate log term + let log_term = unsafe { log(absx + LANCZOS_G - 0.5) }; + + // Calculate the multiplication and subtraction + let step2 = factor * (log_term - 1.0); + + // Combine the results + r += step2; if x < 0.0 { - // Use reflection formula to get value for negative x - - // Calculate the sin(pi * x) value using m_sinpi + // Calculate each component separately as in CPython let sinpi_val = m_sinpi(absx); - - // In CPython, the expression is: - // r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r; - // We'll match this order exactly - r = LOG_PI - sinpi_val.abs().ln() - absx.ln() - r; + let abs_sinpi = unsafe { fabs(sinpi_val) }; + let log_abs_sinpi = unsafe { log(abs_sinpi) }; + let log_absx = unsafe { log(absx) }; + + // Combine in exactly the same order as CPython + r = LOG_PI - log_abs_sinpi - log_absx - r; } + if r.is_infinite() { return Err(Error::ERANGE); } @@ -342,6 +379,8 @@ mod tests { let py_lgamma_repr = unsafe { std::mem::transmute::(py_lgamma) }; let rs_lgamma_repr = unsafe { std::mem::transmute::(rs_lgamma) }; println!("Bit difference: {}", py_lgamma_repr ^ rs_lgamma_repr); + + assert_eq!(py_lgamma_repr, rs_lgamma_repr); }); } From 7fb326e441cc2cf9dd31f826eb0e5436368880af Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Mon, 21 Apr 2025 00:26:04 +0900 Subject: [PATCH 3/5] test_lietral --- src/gamma.rs | 85 +++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 74 insertions(+), 11 deletions(-) diff --git a/src/gamma.rs b/src/gamma.rs index 97bace6..a8ca12d 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -1,7 +1,6 @@ use crate::Error; // use std::f64::consts::PI; - // Import C library functions properly unsafe extern "C" { fn exp(x: f64) -> f64; @@ -197,7 +196,7 @@ pub fn tgamma(x: f64) -> Result { let lanczos = lanczos_sum(absx); let mut r = term3 / lanczos; r -= z * r; - + if absx < 140.0 { unsafe { r / pow(y, absx - 0.5) } } else { @@ -211,7 +210,7 @@ pub fn tgamma(x: f64) -> Result { let exp_y = unsafe { exp(y) }; let mut r = lanczos / exp_y; r += z * r; - + if absx < 140.0 { unsafe { r * pow(y, absx - 0.5) } } else { @@ -261,19 +260,19 @@ pub fn lgamma(x: f64) -> Result { // Using C's math functions through libc to match CPython let lanczos_sum_val = lanczos_sum(absx); let log_lanczos = unsafe { log(lanczos_sum_val) }; - + // Subtract lanczos_g as a separate step let mut r = log_lanczos - LANCZOS_G; - + // Calculate (absx - 0.5) term let factor = absx - 0.5; - + // Calculate log term let log_term = unsafe { log(absx + LANCZOS_G - 0.5) }; - + // Calculate the multiplication and subtraction let step2 = factor * (log_term - 1.0); - + // Combine the results r += step2; @@ -283,11 +282,11 @@ pub fn lgamma(x: f64) -> Result { let abs_sinpi = unsafe { fabs(sinpi_val) }; let log_abs_sinpi = unsafe { log(abs_sinpi) }; let log_absx = unsafe { log(absx) }; - + // Combine in exactly the same order as CPython r = LOG_PI - log_abs_sinpi - log_absx - r; } - + if r.is_infinite() { return Err(Error::ERANGE); } @@ -328,9 +327,73 @@ mod tests { } } + #[test] + fn test_literal() { + // Verify single constants + assert_eq!(PI, 3.141592653589793238462643383279502884197); + assert_eq!(LOG_PI, 1.144729885849400174143427351353058711647); + assert_eq!(LANCZOS_G, 6.024680040776729583740234375); + assert_eq!(LANCZOS_G_MINUS_HALF, 5.524680040776729583740234375); + + // Verify LANCZOS_NUM_COEFFS + assert_eq!(LANCZOS_NUM_COEFFS[0], 23531376880.410759); + assert_eq!(LANCZOS_NUM_COEFFS[1], 42919803642.649101); + assert_eq!(LANCZOS_NUM_COEFFS[2], 35711959237.355667); + assert_eq!(LANCZOS_NUM_COEFFS[3], 17921034426.037209); + assert_eq!(LANCZOS_NUM_COEFFS[4], 6039542586.3520279); + assert_eq!(LANCZOS_NUM_COEFFS[5], 1439720407.3117216); + assert_eq!(LANCZOS_NUM_COEFFS[6], 248874557.86205417); + assert_eq!(LANCZOS_NUM_COEFFS[7], 31426415.585400194); + assert_eq!(LANCZOS_NUM_COEFFS[8], 2876370.6289353725); + assert_eq!(LANCZOS_NUM_COEFFS[9], 186056.26539522348); + assert_eq!(LANCZOS_NUM_COEFFS[10], 8071.6720023658163); + assert_eq!(LANCZOS_NUM_COEFFS[11], 210.82427775157936); + assert_eq!(LANCZOS_NUM_COEFFS[12], 2.5066282746310002); + + // Verify LANCZOS_DEN_COEFFS + assert_eq!(LANCZOS_DEN_COEFFS[0], 0.0); + assert_eq!(LANCZOS_DEN_COEFFS[1], 39916800.0); + assert_eq!(LANCZOS_DEN_COEFFS[2], 120543840.0); + assert_eq!(LANCZOS_DEN_COEFFS[3], 150917976.0); + assert_eq!(LANCZOS_DEN_COEFFS[4], 105258076.0); + assert_eq!(LANCZOS_DEN_COEFFS[5], 45995730.0); + assert_eq!(LANCZOS_DEN_COEFFS[6], 13339535.0); + assert_eq!(LANCZOS_DEN_COEFFS[7], 2637558.0); + assert_eq!(LANCZOS_DEN_COEFFS[8], 357423.0); + assert_eq!(LANCZOS_DEN_COEFFS[9], 32670.0); + assert_eq!(LANCZOS_DEN_COEFFS[10], 1925.0); + assert_eq!(LANCZOS_DEN_COEFFS[11], 66.0); + assert_eq!(LANCZOS_DEN_COEFFS[12], 1.0); + + // Verify GAMMA_INTEGRAL + assert_eq!(GAMMA_INTEGRAL[0], 1.0); + assert_eq!(GAMMA_INTEGRAL[1], 1.0); + assert_eq!(GAMMA_INTEGRAL[2], 2.0); + assert_eq!(GAMMA_INTEGRAL[3], 6.0); + assert_eq!(GAMMA_INTEGRAL[4], 24.0); + assert_eq!(GAMMA_INTEGRAL[5], 120.0); + assert_eq!(GAMMA_INTEGRAL[6], 720.0); + assert_eq!(GAMMA_INTEGRAL[7], 5040.0); + assert_eq!(GAMMA_INTEGRAL[8], 40320.0); + assert_eq!(GAMMA_INTEGRAL[9], 362880.0); + assert_eq!(GAMMA_INTEGRAL[10], 3628800.0); + assert_eq!(GAMMA_INTEGRAL[11], 39916800.0); + assert_eq!(GAMMA_INTEGRAL[12], 479001600.0); + assert_eq!(GAMMA_INTEGRAL[13], 6227020800.0); + assert_eq!(GAMMA_INTEGRAL[14], 87178291200.0); + assert_eq!(GAMMA_INTEGRAL[15], 1307674368000.0); + assert_eq!(GAMMA_INTEGRAL[16], 20922789888000.0); + assert_eq!(GAMMA_INTEGRAL[17], 355687428096000.0); + assert_eq!(GAMMA_INTEGRAL[18], 6402373705728000.0); + assert_eq!(GAMMA_INTEGRAL[19], 1.21645100408832e+17); + assert_eq!(GAMMA_INTEGRAL[20], 2.43290200817664e+18); + assert_eq!(GAMMA_INTEGRAL[21], 5.109094217170944e+19); + assert_eq!(GAMMA_INTEGRAL[22], 1.1240007277776077e+21); + } + #[test] fn test_specific_lgamma_value() { - let x = -3.8510064710745118; + let x = 0.003585187864492183; let rs_lgamma = lgamma(x).unwrap(); pyo3::prepare_freethreaded_python(); From ed79cc2a8e965ca611bad3f66495c7c307a51825 Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Mon, 21 Apr 2025 14:12:34 +0900 Subject: [PATCH 4/5] tgamma.c --- tgamma.c | 316 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 316 insertions(+) create mode 100644 tgamma.c diff --git a/tgamma.c b/tgamma.c new file mode 100644 index 0000000..691970b --- /dev/null +++ b/tgamma.c @@ -0,0 +1,316 @@ +/* + * Implementation of the gamma function based on CPython's implementation. + * This is a standalone version without dependencies on CPython. + */ + +#include +#include +#include +#include +#include + +/* Constants and definitions needed for gamma calculation */ +static const double pi = 3.141592653589793238462643383279502884197; +static const double logpi = 1.144729885849400174143427351353058711647; + +/* Check if a value is finite */ +#ifndef Py_IS_FINITE +#define Py_IS_FINITE(x) isfinite(x) +#endif + +/* Check if a value is infinity */ +#ifndef Py_IS_INFINITY +#define Py_IS_INFINITY(x) isinf(x) +#endif + +/* Check if a value is NaN */ +#ifndef Py_IS_NAN +#define Py_IS_NAN(x) isnan(x) +#endif + +/* Define HUGE_VAL if not available */ +#ifndef Py_HUGE_VAL +#define Py_HUGE_VAL HUGE_VAL +#endif + +/* Define mathematical constants required for Lanczos approximation */ +#define LANCZOS_N 13 +static const double lanczos_g = 6.024680040776729583740234375; +static const double lanczos_g_minus_half = 5.524680040776729583740234375; +static const double lanczos_num_coeffs[LANCZOS_N] = { + 23531376880.410759688572007674451636754734846804940, + 42919803642.649098768957899047001988850926355848959, + 35711959237.355668049440185451547166705960488635843, + 17921034426.037209699919755754458931112671403265390, + 6039542586.3520280050642916443072979210699388420708, + 1439720407.3117216736632230727949123939715485786772, + 248874557.86205415651146038641322942321632125127801, + 31426415.585400194380614231628318205362874684987640, + 2876370.6289353724412254090516208496135991145378768, + 186056.26539522349504029498971604569928220784236328, + 8071.6720023658162106380029022722506138218516325024, + 210.82427775157934587250973392071336271166969580291, + 2.5066282746310002701649081771338373386264310793408 +}; + +/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ +static const double lanczos_den_coeffs[LANCZOS_N] = { + 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, + 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; + +/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ +#define NGAMMA_INTEGRAL 23 +static const double gamma_integral[NGAMMA_INTEGRAL] = { + 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, + 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, + 1307674368000.0, 20922789888000.0, 355687428096000.0, + 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, + 51090942171709440000.0, 1124000727777607680000.0, +}; + +/* Helper function for sin(pi*x) used in gamma calculation */ +static double m_sinpi(double x) +{ + double y, r; + int n; + /* this function should only ever be called for finite arguments */ + assert(Py_IS_FINITE(x)); + y = fmod(fabs(x), 2.0); + n = (int)round(2.0*y); + assert(0 <= n && n <= 4); + switch (n) { + case 0: + r = sin(pi*y); + break; + case 1: + r = cos(pi*(y-0.5)); + break; + case 2: + /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give + -0.0 instead of 0.0 when y == 1.0. */ + r = sin(pi*(1.0-y)); + break; + case 3: + r = -cos(pi*(y-1.5)); + break; + case 4: + r = sin(pi*(y-2.0)); + break; + default: + /* Should never happen due to the assert above */ + r = 0.0; + } + return copysign(1.0, x)*r; +} + +/* Implementation of Lanczos sum for gamma function */ +static double lanczos_sum(double x) +{ + double num = 0.0, den = 0.0; + int i; + assert(x > 0.0); + /* evaluate the rational function lanczos_sum(x). For large + x, the obvious algorithm risks overflow, so we instead + rescale the denominator and numerator of the rational + function by x**(1-LANCZOS_N) and treat this as a + rational function in 1/x. This also reduces the error for + larger x values. The choice of cutoff point (5.0 below) is + somewhat arbitrary; in tests, smaller cutoff values than + this resulted in lower accuracy. */ + if (x < 5.0) { + for (i = LANCZOS_N; --i >= 0; ) { + num = num * x + lanczos_num_coeffs[i]; + den = den * x + lanczos_den_coeffs[i]; + } + } + else { + for (i = 0; i < LANCZOS_N; i++) { + num = num / x + lanczos_num_coeffs[i]; + den = den / x + lanczos_den_coeffs[i]; + } + } + return num/den; +} + +/** + * Computes the gamma function value for x. + * + * The implementation is based on the Lanczos approximation with parameters + * N=13 and g=6.024680040776729583740234375, which provides excellent accuracy + * across the domain of the function. + * + * @param x The input value + * @return The gamma function value + * + * Special cases: + * - If x is NaN, returns NaN + * - If x is +Inf, returns +Inf + * - If x is -Inf, sets errno to EDOM and returns NaN + * - If x is 0, sets errno to EDOM and returns +/-Inf (with the sign of x) + * - If x is a negative integer, sets errno to EDOM and returns NaN + */ +double tgamma(double x) { + double absx, r, y, z, sqrtpow; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x) || x > 0.0) + return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ + else { + errno = EDOM; + return NAN; /* tgamma(-inf) = nan, invalid */ + } + } + if (x == 0.0) { + errno = EDOM; + /* tgamma(+-0.0) = +-inf, divide-by-zero */ + return copysign(INFINITY, x); + } + + /* integer arguments */ + if (x == floor(x)) { + if (x < 0.0) { + errno = EDOM; /* tgamma(n) = nan, invalid for */ + return NAN; /* negative integers n */ + } + if (x <= NGAMMA_INTEGRAL) + return gamma_integral[(int)x - 1]; + } + absx = fabs(x); + + /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ + if (absx < 1e-20) { + r = 1.0/x; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; + } + + /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for + x > 200, and underflows to +-0.0 for x < -200, not a negative + integer. */ + if (absx > 200.0) { + if (x < 0.0) { + return 0.0/m_sinpi(x); + } + else { + errno = ERANGE; + return Py_HUGE_VAL; + } + } + + y = absx + lanczos_g_minus_half; + /* compute error in sum */ + if (absx > lanczos_g_minus_half) { + /* note: the correction can be foiled by an optimizing + compiler that (incorrectly) thinks that an expression like + a + b - a - b can be optimized to 0.0. This shouldn't + happen in a standards-conforming compiler. */ + double q = y - absx; + z = q - lanczos_g_minus_half; + } + else { + double q = y - lanczos_g_minus_half; + z = q - absx; + } + z = z * lanczos_g / y; + if (x < 0.0) { + r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx); + r -= z * r; + if (absx < 140.0) { + r /= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r /= sqrtpow; + r /= sqrtpow; + } + } + else { + r = lanczos_sum(absx) / exp(y); + r += z * r; + if (absx < 140.0) { + r *= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r *= sqrtpow; + r *= sqrtpow; + } + } + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; +} + +/** + * Computes the natural logarithm of the absolute value of the gamma function. + * + * @param x The input value + * @return The log of the absolute gamma function value + * + * Special cases: + * - If x is NaN, returns NaN + * - If x is +/-Inf, returns +Inf + * - If x is a non-positive integer, sets errno to EDOM and returns +Inf + * - If x is 1 or 2, returns 0.0 + */ +double lgamma(double x) { + double r; + double absx; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x)) + return x; /* lgamma(nan) = nan */ + else + return HUGE_VAL; /* lgamma(+-inf) = +inf */ + } + + /* integer arguments */ + if (x == floor(x) && x <= 2.0) { + if (x <= 0.0) { + errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ + return HUGE_VAL; /* integers n <= 0 */ + } + else { + return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ + } + } + + absx = fabs(x); + /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ + if (absx < 1e-20) + return -log(absx); + + /* Lanczos' formula. We could save a fraction of a ulp in accuracy by + having a second set of numerator coefficients for lanczos_sum that + absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g + subtraction below; it's probably not worth it. */ + r = log(lanczos_sum(absx)) - lanczos_g; + r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1); + if (x < 0.0) + /* Use reflection formula to get value for negative x. */ + r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; +} + +#include +#include + +union Result { + double d; + uint64_t u; +}; + +int main() { + union Result result; + result.d = tgamma(-3.8510064710745118); + printf("The result of tgamma(-3.8510064710745118) is: %f\n", result.d); + + // Print the binary representation of a double + printf("Bit representation of result: %llu\n", result.u); + return 0; +} \ No newline at end of file From 5cff60fd06d16d9f11edae06340ec555c6bddbff Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Mon, 21 Apr 2025 14:35:43 +0900 Subject: [PATCH 5/5] Add macos again --- .github/workflows/rust.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 98ec583..654952a 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -19,7 +19,7 @@ jobs: os: [ ubuntu-latest, windows-latest, - # macos-latest # disabled due to incompatibility. See issue #1 + macos-latest, ] rust: [stable] steps: