From 7622c3b3307c708764b5b91c8bc0390e1181ca20 Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Mon, 21 Apr 2025 16:28:09 +0900 Subject: [PATCH] libm --- Cargo.toml | 6 +- proptest-regressions/gamma.txt | 3 + proptest-regressions/lib.txt | 7 ++ src/err.rs | 16 ++++- src/gamma.rs | 88 +++--------------------- src/lib.rs | 119 ++++++++++++++++++++++++++++++++- src/m.rs | 65 ++++++++++++++++++ src/test.rs | 28 ++++++++ 8 files changed, 250 insertions(+), 82 deletions(-) create mode 100644 proptest-regressions/lib.txt create mode 100644 src/m.rs create mode 100644 src/test.rs diff --git a/Cargo.toml b/Cargo.toml index e95f84d..5d3216f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "pymath" -version = "0.1.0" +version = "0.0.1" edition = "2024" [features] @@ -8,6 +8,10 @@ edition = "2024" # 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"] } diff --git a/proptest-regressions/gamma.txt b/proptest-regressions/gamma.txt index 150f6e6..77162bb 100644 --- a/proptest-regressions/gamma.txt +++ b/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 diff --git a/proptest-regressions/lib.txt b/proptest-regressions/lib.txt new file mode 100644 index 0000000..25c7cb0 --- /dev/null +++ b/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 diff --git a/src/err.rs b/src/err.rs index bbd3b58..2d8e0e9 100644 --- a/src/err.rs +++ b/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 = std::result::Result; + +impl TryFrom for Error { + type Error = libc::c_int; + + fn try_from(value: libc::c_int) -> std::result::Result { + match value { + 33 => Ok(Error::EDOM), + 34 => Ok(Error::ERANGE), + _ => Err(value), + } + } +} diff --git a/src/gamma.rs b/src/gamma.rs index db96333..534f7ff 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -117,7 +117,8 @@ const GAMMA_INTEGRAL: [f64; NGAMMA_INTEGRAL] = [ 1124000727777607680000.0, ]; -pub fn tgamma(x: f64) -> Result { +// tgamma +pub fn gamma(x: f64) -> crate::Result { // special cases if !x.is_finite() { if x.is_nan() || x > 0.0 { @@ -213,7 +214,7 @@ pub fn tgamma(x: f64) -> Result { // 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 { +pub fn lgamma(x: f64) -> crate::Result { // special cases if !x.is_finite() { if x.is_nan() { @@ -258,79 +259,10 @@ pub fn lgamma(x: f64) -> Result { 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>, - v: Result, - ) -> 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::(py) { - assert_eq!(v.err(), Some(Error::EDOM)); - } else if e.is_instance_of::(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 +); diff --git a/src/lib.rs b/src/lib.rs index 1a884f7..14363cb 100644 --- a/src/lib.rs +++ b/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 { +// 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; diff --git a/src/m.rs b/src/m.rs new file mode 100644 index 0000000..0dff12a --- /dev/null +++ b/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; + // }} +} diff --git a/src/test.rs b/src/test.rs new file mode 100644 index 0000000..7869f1e --- /dev/null +++ b/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>, + v: Result, +) -> 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::(py) { + assert_eq!(v.err(), Some(Error::EDOM)); + } else if e.is_instance_of::(py) { + assert_eq!(v.err(), Some(Error::ERANGE)); + } else { + panic!(); + } + None + } + } +}