erf, erfc with libm

Cargo.toml

@@ -1,13 +1,17 @@
 [package]
 name = "pymath"
-version = "0.1.0"
+version = "0.0.1"
 edition = "2024"
 
 [features]
 # Turning on this feature on aarch64-apple-darwin helps bit representation compatibility
 # See also: https://github.com/python/cpython/issues/132763
 mul_add = []
 
+[dependencies]
+errno = "0.3"
+libc = "0.2"
+
 [dev-dependencies]
 proptest = "1.6.0"
 pyo3 = { version = "0.24", features = ["abi3"] }

proptest-regressions/gamma.txt

@@ -8,3 +8,6 @@ cc e8ed768221998086795d95c68921437e80c4b7fe68fe9da15ca40faa216391b5 # shrinks to
 cc 23c7f86ab299daa966772921d8c615afda11e1b77944bed40e88264a68e62ac3 # shrinks to x = -19.80948467648103
 cc f57954d91904549b9431755f196b630435a43cbefd558b932efad487a403c6c8 # shrinks to x = 0.003585187864492183
 cc 7a9a04aed4ed7e3d23eb7b32b748542b1062e349ae83cc1fad39672a5b2156cd # shrinks to x = -3.8510064710745118
+cc d884d4ef56bcd40d025660e0dec152754fd4fd4e48bc0bdf41e73ea001798fd8 # shrinks to x = 0.9882904125102558
+cc 3f1d36f364ce29810d0c37003465b186c07460861c7a3bf4b8962401b376f2d9 # shrinks to x = 1.402608516799205
+cc 4439ce674d91257d104063e2d5ade7908c83462d195f98a0c304ea25b022d0f4 # shrinks to x = 3.6215752811868267

proptest-regressions/lib.txt

@@ -0,0 +1,7 @@
+# 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 531a136f9fcde9d1da1ba5d173e62eee8ec8f7c877eb34abbc6d47611a641bc7 # shrinks to x = 0.0

src/err.rs

@@ -1,6 +1,20 @@
-// defined in libc
+// The values are defined in libc
 #[derive(Debug, PartialEq, Eq)]
 pub enum Error {
     EDOM = 33,
     ERANGE = 34,
 }
+
+pub type Result<T> = std::result::Result<T, Error>;
+
+impl TryFrom<libc::c_int> for Error {
+    type Error = libc::c_int;
+
+    fn try_from(value: libc::c_int) -> std::result::Result<Self, Self::Error> {
+        match value {
+            33 => Ok(Error::EDOM),
+            34 => Ok(Error::ERANGE),
+            _ => Err(value),
+        }
+    }
+}

src/gamma.rs

@@ -117,7 +117,8 @@ const GAMMA_INTEGRAL: [f64; NGAMMA_INTEGRAL] = [
     1124000727777607680000.0,
 ];
 
-pub fn tgamma(x: f64) -> Result<f64, Error> {
+// tgamma
+pub fn gamma(x: f64) -> crate::Result<f64> {
     // special cases
     if !x.is_finite() {
         if x.is_nan() || x > 0.0 {
@@ -213,7 +214,7 @@ pub fn tgamma(x: f64) -> Result<f64, Error> {
 
 // natural log of the absolute value of the Gamma function.
 // For large arguments, Lanczos' formula works extremely well here.
-pub fn lgamma(x: f64) -> Result<f64, Error> {
+pub fn lgamma(x: f64) -> crate::Result<f64> {
     // special cases
     if !x.is_finite() {
         if x.is_nan() {
@@ -258,79 +259,10 @@ pub fn lgamma(x: f64) -> Result<f64, Error> {
     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 = py_gamma.to_bits();
-                let rs_gamma_repr = rs_gamma.to_bits();
-                assert_eq!(py_gamma_repr, rs_gamma_repr, "x = {x}, py_gamma = {py_gamma}, rs_gamma = {rs_gamma}");
-            });
-        }
-
-        #[test]
-        fn test_lgamma(x: f64) {
-            let rs_lgamma = lgamma(x);
-
-            pyo3::prepare_freethreaded_python();
-            Python::with_gil(|py| {
-                let math = PyModule::import(py, "math").unwrap();
-                let py_lgamma_func = math
-                    .getattr("lgamma")
-                    .unwrap();
-                let r = py_lgamma_func.call1((x,));
-                let Some((py_lgamma, rs_lgamma)) = unwrap(py, r, rs_lgamma) else {
-                    return;
-                };
-                let py_lgamma_repr = py_lgamma.to_bits();
-                let rs_lgamma_repr = rs_lgamma.to_bits();
-                assert_eq!(py_lgamma_repr, rs_lgamma_repr, "x = {x}, py_lgamma = {py_lgamma}, rs_gamma = {rs_lgamma}");
-            });
-        }
-    }
-}
+super::pyo3_proptest!(gamma(Result<_>), test_gamma, proptest_gamma, fulltest_gamma);
+super::pyo3_proptest!(
+    lgamma(Result<_>),
+    test_lgamma,
+    proptest_lgamma,
+    fulltest_lgamma
+);

src/lib.rs

@@ -1,5 +1,120 @@
 mod err;
 mod gamma;
+mod m;
+#[cfg(test)]
+mod test;
 
-pub use err::Error;
-pub use gamma::{lgamma, tgamma as gamma};
+pub use err::{Error, Result};
+pub use gamma::{gamma, lgamma};
+
+macro_rules! libm {
+    // Reset errno and handle errno when return type contains Result
+    (fn $name:ident($arg:ident: $ty:ty) -> Result<$ret:ty>) => {
+        #[inline(always)]
+        pub fn $name($arg: $ty) -> Result<$ret> {
+            errno::set_errno(errno::Errno(0));
+            let r = unsafe { m::$name($arg) };
+            crate::is_error(r)
+        }
+    };
+    // Skip errno checking when return type is not Result
+    (fn $name:ident($arg:ident: $ty:ty) -> $ret:ty) => {
+        #[inline(always)]
+        pub fn $name($arg: $ty) -> $ret {
+            unsafe { m::$name($arg) }
+        }
+    };
+}
+
+macro_rules! pyo3_proptest {
+    ($fn_name:ident(Result<_>), $test_name:ident, $proptest_name:ident, $edgetest_name:ident) => {
+        #[cfg(test)]
+        fn $test_name(x: f64) {
+            use pyo3::prelude::*;
+
+            let rs_result = $fn_name(x);
+
+            pyo3::prepare_freethreaded_python();
+            Python::with_gil(|py| {
+                let math = PyModule::import(py, "math").unwrap();
+                let py_func = math
+                    .getattr(stringify!($fn_name))
+                    .unwrap();
+                let r = py_func.call1((x,));
+                let Some((py_result, rs_result)) = crate::test::unwrap(py, r, rs_result) else {
+                    return;
+                };
+                let py_result_repr = py_result.to_bits();
+                let rs_result_repr = rs_result.to_bits();
+                assert_eq!(py_result_repr, rs_result_repr, "x = {x}, py_result = {py_result}, rs_result = {rs_result}");
+            });
+        }
+
+        crate::pyo3_proptest!(@proptest, $test_name, $proptest_name);
+        crate::pyo3_proptest!(@edgetest, $test_name, $edgetest_name);
+    };
+    ($fn_name:ident(_), $test_name:ident, $proptest_name:ident, $edgetest_name:ident) => {
+        #[cfg(test)]
+        fn $test_name(x: f64) {
+            use pyo3::prelude::*;
+
+            let rs_result = Ok($fn_name(x));
+
+            pyo3::prepare_freethreaded_python();
+            Python::with_gil(|py| {
+                let math = PyModule::import(py, "math").unwrap();
+                let py_func = math
+                    .getattr(stringify!($fn_name))
+                    .unwrap();
+                let r = py_func.call1((x,));
+                let Some((py_result, rs_result)) = crate::test::unwrap(py, r, rs_result) else {
+                    return;
+                };
+                let py_result_repr = py_result.to_bits();
+                let rs_result_repr = rs_result.to_bits();
+                assert_eq!(py_result_repr, rs_result_repr, "x = {x}, py_result = {py_result}, rs_result = {rs_result}");
+            });
+        }
+        crate::pyo3_proptest!(@proptest, $test_name, $proptest_name);
+    };
+    (@proptest, $test_name:ident, $proptest_name:ident) => {
+        #[cfg(test)]
+        proptest::proptest! {
+            #[test]
+            fn $proptest_name(x: f64) {
+                $test_name(x);
+            }
+        }
+    };
+    (@edgetest, $test_name:ident, $edgetest_name:ident) => {
+        #[test]
+        fn $edgetest_name() {
+            $test_name(f64::MIN);
+            $test_name(-f64::MIN);
+            $test_name(f64::NAN);
+            $test_name(-f64::NAN);
+            $test_name(f64::INFINITY);
+            $test_name(-f64::NEG_INFINITY);
+            $test_name(0.0);
+            $test_name(-0.0);
+        }
+    };
+}
+
+libm!(fn erf(n: f64) -> f64);
+pyo3_proptest!(erf(_), test_erf, proptest_erf, edgetest_erf);
+
+libm!(fn erfc(n: f64) -> f64);
+pyo3_proptest!(erfc(_), test_erfc, proptest_erfc, edgetest_erfc);
+
+/// Call is_error when errno != 0, and where x is the result libm
+/// returned.  is_error will usually set up an exception and return
+/// true (1), but may return false (0) without setting up an exception.
+// fn is_error(x: f64) -> crate::Result<f64> {
+//     match errno::errno() {
+//         errno::Errno(0) => Ok(x),
+//         errno::Errno(libc::ERANGE) if x.abs() < 1.5 => Ok(0f64),
+//         errno::Errno(errno) => Err(errno.try_into().unwrap()),
+//     }
+// }
+use pyo3_proptest;

src/m.rs

@@ -0,0 +1,65 @@
+//! Partial copy of std::sys::_cmath
+
+// These symbols are all defined by `libm`,
+// or by `compiler-builtins` on unsupported platforms.
+#[allow(dead_code)]
+unsafe extern "C" {
+    pub fn acos(n: f64) -> f64;
+    pub fn asin(n: f64) -> f64;
+    pub fn atan(n: f64) -> f64;
+    pub fn atan2(a: f64, b: f64) -> f64;
+    pub fn cbrt(n: f64) -> f64;
+    pub fn cbrtf(n: f32) -> f32;
+    pub fn cosh(n: f64) -> f64;
+    pub fn expm1(n: f64) -> f64;
+    pub fn expm1f(n: f32) -> f32;
+    pub fn fdim(a: f64, b: f64) -> f64;
+    pub fn fdimf(a: f32, b: f32) -> f32;
+    #[cfg_attr(target_env = "msvc", link_name = "_hypot")]
+    pub fn hypot(x: f64, y: f64) -> f64;
+    #[cfg_attr(target_env = "msvc", link_name = "_hypotf")]
+    pub fn hypotf(x: f32, y: f32) -> f32;
+    pub fn log1p(n: f64) -> f64;
+    pub fn log1pf(n: f32) -> f32;
+    pub fn sinh(n: f64) -> f64;
+    pub fn tan(n: f64) -> f64;
+    pub fn tanh(n: f64) -> f64;
+    pub fn tgamma(n: f64) -> f64;
+    pub fn tgammaf(n: f32) -> f32;
+    pub fn lgamma_r(n: f64, s: &mut i32) -> f64;
+    #[cfg(not(target_os = "aix"))]
+    pub fn lgammaf_r(n: f32, s: &mut i32) -> f32;
+    pub fn erf(n: f64) -> f64;
+    pub fn erff(n: f32) -> f32;
+    pub fn erfc(n: f64) -> f64;
+    pub fn erfcf(n: f32) -> f32;
+
+    // pub fn acosf128(n: f128) -> f128;
+    // pub fn asinf128(n: f128) -> f128;
+    // pub fn atanf128(n: f128) -> f128;
+    // pub fn atan2f128(a: f128, b: f128) -> f128;
+    // pub fn cbrtf128(n: f128) -> f128;
+    // pub fn coshf128(n: f128) -> f128;
+    // pub fn expm1f128(n: f128) -> f128;
+    // pub fn hypotf128(x: f128, y: f128) -> f128;
+    // pub fn log1pf128(n: f128) -> f128;
+    // pub fn sinhf128(n: f128) -> f128;
+    // pub fn tanf128(n: f128) -> f128;
+    // pub fn tanhf128(n: f128) -> f128;
+    // pub fn tgammaf128(n: f128) -> f128;
+    // pub fn lgammaf128_r(n: f128, s: &mut i32) -> f128;
+    // pub fn erff128(n: f128) -> f128;
+    // pub fn erfcf128(n: f128) -> f128;
+
+    // cfg_if::cfg_if! {
+    // if #[cfg(not(all(target_os = "windows", target_env = "msvc", target_arch = "x86")))] {
+    //     pub fn acosf(n: f32) -> f32;
+    //     pub fn asinf(n: f32) -> f32;
+    //     pub fn atan2f(a: f32, b: f32) -> f32;
+    //     pub fn atanf(n: f32) -> f32;
+    //     pub fn coshf(n: f32) -> f32;
+    //     pub fn sinhf(n: f32) -> f32;
+    //     pub fn tanf(n: f32) -> f32;
+    //     pub fn tanhf(n: f32) -> f32;
+    // }}
+}

src/test.rs

@@ -0,0 +1,28 @@
+use crate::Error;
+use pyo3::{Python, prelude::*};
+
+pub(crate) 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
+        }
+    }
+}