From c34f7bd6eae7b207dca82fcb3155fa3fb22ecde8 Mon Sep 17 00:00:00 2001
From: Jeong YunWon <jeong@youknowone.org>
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<f64, Error> {
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<f64, Error> {
}
r
};
+
if r.is_infinite() {
return Err((f64::INFINITY, Error::ERANGE).1);
} else {
@@ -241,7 +242,14 @@ pub fn lgamma(x: f64) -> Result<f64, Error> {
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::<f64>()
+ .unwrap();
+
+ println!("x = {}", x);
+ println!("Python lgamma = {} ({:x})", py_lgamma, unsafe {
+ std::mem::transmute::<f64, u64>(py_lgamma)
+ });
+ println!("Rust lgamma = {} ({:x})", rs_lgamma, unsafe {
+ std::mem::transmute::<f64, u64>(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::<f64>()
+ .unwrap();
+ println!("Python sinpi = {}", py_sinpi);
+
+ let py_lgamma_repr = unsafe { std::mem::transmute::<f64, u64>(py_lgamma) };
+ let rs_lgamma_repr = unsafe { std::mem::transmute::<f64, u64>(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 <jeong@youknowone.org>
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 <stdio.h>
+#include <stdint.h>
+#include <string.h> // For memcpy
+#include <math.h> // 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<f64, Error> {
// special cases
if !x.is_finite() {
@@ -130,7 +145,7 @@ pub fn tgamma(x: f64) -> Result<f64, Error> {
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<f64, Error> {
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<f64, Error> {
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<f64, Error> {
}
// 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<f64, Error> {
}
}
- 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::<f64, u64>(py_lgamma) };
let rs_lgamma_repr = unsafe { std::mem::transmute::<f64, u64>(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 <jeong@youknowone.org>
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<f64, Error> {
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<f64, Error> {
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<f64, Error> {
// 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<f64, Error> {
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 <jeong@youknowone.org>
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 <math.h>
+#include <errno.h>
+#include <float.h>
+#include <assert.h>
+#include <stdbool.h>
+
+/* 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 <stdio.h>
+#include <stdint.h>
+
+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 <jeong@youknowone.org>
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: