start a crate with tgamma

Author: youknowone

.gitignore

@@ -0,0 +1,2 @@
+/target
+Cargo.lock

Cargo.toml

@@ -0,0 +1,8 @@
+[package]
+name = "pymath"
+version = "0.1.0"
+edition = "2024"
+
+[dev-dependencies]
+proptest = "1.6.0"
+pyo3 = "0.23.4"

proptest-regressions/gamma.txt

@@ -0,0 +1,9 @@
+# Seeds for failure cases proptest has generated in the past. It is
+# automatically read and these particular cases re-run before any
+# novel cases are generated.
+#
+# It is recommended to check this file in to source control so that
+# everyone who runs the test benefits from these saved cases.
+cc e8ed768221998086795d95c68921437e80c4b7fe68fe9da15ca40faa216391b5 # shrinks to x = 0.0
+cc 23c7f86ab299daa966772921d8c615afda11e1b77944bed40e88264a68e62ac3 # shrinks to x = -19.80948467648103
+cc f57954d91904549b9431755f196b630435a43cbefd558b932efad487a403c6c8 # shrinks to x = 0.003585187864492183

src/err.rs

@@ -0,0 +1,6 @@
+// defined in libc
+#[derive(Debug, PartialEq, Eq)]
+pub enum Error {
+    EDOM = 33,
+    ERANGE = 34,
+}

src/gamma.rs

@@ -0,0 +1,263 @@
+use crate::Error;
+use std::f64::consts::PI;
+
+const LOG_PI: f64 = 1.144729885849400174143427351353058711647;
+
+const LANCZOS_N: usize = 13;
+const LANCZOS_G: f64 = 6.024680040776729583740234375;
+const LANCZOS_G_MINUS_HALF: f64 = 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,
+];
+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,
+];
+
+fn lanczos_sum(x: f64) -> f64 {
+    let mut num = 0.0;
+    let mut den = 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 in (0..LANCZOS_N).rev() {
+            num = num * x + LANCZOS_NUM_COEFFS[i];
+            den = den * x + LANCZOS_DEN_COEFFS[i];
+        }
+    } else {
+        for i in 0..LANCZOS_N {
+            num = num / x + LANCZOS_NUM_COEFFS[i];
+            den = den / x + LANCZOS_DEN_COEFFS[i];
+        }
+    }
+    num / den
+}
+
+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 r = match n {
+        0 => (PI * y).sin(),
+        1 => (PI * (y - 0.5)).cos(),
+        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()
+        }
+        3 => -(PI * (y - 1.5)).cos(),
+        4 => (PI * (y - 2.0)).sin(),
+        _ => 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() {
+        if x.is_nan() || x > 0.0 {
+            // tgamma(nan) = nan, tgamma(inf) = inf
+            return Ok(x);
+        } else {
+            // tgamma(-inf) = nan, invalid
+            return Err((f64::NAN, Error::EDOM).1);
+        }
+    }
+    if x == 0.0 {
+        // tgamma(+-0.0) = +-inf, divide-by-zero
+        let v = if x.is_sign_positive() {
+            f64::INFINITY
+        } else {
+            f64::NEG_INFINITY
+        };
+        return Err((v, Error::EDOM).1);
+    }
+    // integer arguments
+    if x == x.floor() {
+        if x < 0.0 {
+            // tgamma(n) = nan, invalid for
+            return Err((f64::NAN, Error::EDOM).1);
+        }
+        if x < NGAMMA_INTEGRAL as f64 {
+            return Ok(GAMMA_INTEGRAL[x as usize - 1]);
+        }
+    }
+    let absx = x.abs();
+    // tiny arguments:  tgamma(x) ~ 1/x for x near 0
+    if absx < 1e-20 {
+        let r = 1.0 / x;
+        if r.is_infinite() {
+            return Err((f64::INFINITY, Error::ERANGE).1);
+        } else {
+            return Ok(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 Ok(0.0 / m_sinpi(x));
+        } else {
+            return Err((f64::INFINITY, Error::ERANGE).1);
+        }
+    }
+
+    let y = absx + LANCZOS_G_MINUS_HALF;
+    let z = 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.
+        let q = y - absx;
+        q - LANCZOS_G_MINUS_HALF
+    } else {
+        let q = y - LANCZOS_G_MINUS_HALF;
+        q - absx
+    };
+    let z = z * LANCZOS_G / y;
+    let 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 {
+            r /= y.powf(absx - 0.5);
+        } else {
+            let sqrtpow = y.powf(absx / 2.0 - 0.25);
+            r /= sqrtpow;
+            r /= sqrtpow;
+        }
+        r
+    } else {
+        let mut r = lanczos_sum(absx) / y.exp();
+        r += z * r;
+        if absx < 140.0 {
+            r *= y.powf(absx - 0.5);
+        } else {
+            let sqrtpow = y.powf(absx / 2.0 - 0.25);
+            r *= sqrtpow;
+            r *= sqrtpow;
+        }
+        r
+    };
+    if r.is_infinite() {
+        return Err((f64::INFINITY, Error::ERANGE).1);
+    } else {
+        return Ok(r);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use pyo3::Python;
+    use pyo3::prelude::*;
+
+    use proptest::prelude::*;
+
+    fn unwrap<'a, T: 'a>(
+        py: Python,
+        py_v: PyResult<Bound<'a, PyAny>>,
+        v: Result<T, crate::Error>,
+    ) -> Option<(T, T)>
+    where
+        T: PartialEq + std::fmt::Debug + FromPyObject<'a>,
+    {
+        match py_v {
+            Ok(py_v) => {
+                let py_v: T = py_v.extract().unwrap();
+                Some((py_v, v.unwrap()))
+            }
+            Err(e) => {
+                if e.is_instance_of::<pyo3::exceptions::PyValueError>(py) {
+                    assert_eq!(v.err(), Some(Error::EDOM));
+                } else if e.is_instance_of::<pyo3::exceptions::PyOverflowError>(py) {
+                    assert_eq!(v.err(), Some(Error::ERANGE));
+                } else {
+                    panic!();
+                }
+                None
+            }
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn test_tgamma(x: f64) {
+            let rs_gamma = tgamma(x);
+
+            pyo3::prepare_freethreaded_python();
+            Python::with_gil(|py| {
+                let math = PyModule::import(py, "math").unwrap();
+                let py_gamma_func = math
+                    .getattr("gamma")
+                    .unwrap();
+                let r = py_gamma_func.call1((x,));
+                let Some((py_gamma, rs_gamma)) = unwrap(py, r, rs_gamma) else {
+                    return;
+                };
+                let py_gamma_repr = unsafe { std::mem::transmute::<f64, i64>(py_gamma) };
+                let rs_gamma_repr = unsafe { std::mem::transmute::<f64, i64>(rs_gamma) };
+                // assert_eq!(py_gamma_repr, rs_gamma_repr, "x = {x}, py_gamma = {py_gamma}, rs_gamma = {rs_gamma}");
+                // allow 1 bit error for now
+                assert!((py_gamma_repr - rs_gamma_repr).abs() <= 1, "x = {x} diff: {}, py_gamma = {py_gamma} ({py_gamma_repr:x}), rs_gamma = {rs_gamma} ({rs_gamma_repr:x})", py_gamma_repr ^ rs_gamma_repr);
+            });
+        }
+    }
+}

src/lib.rs

@@ -0,0 +1,5 @@
+mod err;
+mod gamma;
+
+pub use err::Error;
+pub use gamma::tgamma as gamma;