math module

Author: youknowone

.github/workflows/rust.yml

@@ -43,3 +43,7 @@ jobs:
     - name: Run tests with FMA (macOS)
       if: matrix.os == 'macos-latest'
       run: cargo test --verbose --features mul_add && cargo test --release --verbose --features mul_add
+    - name: Run tests with num-bigint
+      run: cargo test --verbose --features num-bigint
+    - name: Run tests with malachite-bigint
+      run: cargo test --verbose --features malachite-bigint

Cargo.toml

@@ -1,12 +1,18 @@
 [package]
 name = "pymath"
-version = "0.0.2"
+version = "0.0.3"
 edition = "2024"
 description = "A binary representation compatible Rust implementation of Python's math library."
 license = "PSF-2.0"
 
 
 [features]
+default = ["mul_add"]
+# complex = ["dep:num-complex"]
+num-bigint = ["_bigint", "dep:num-bigint"]
+malachite-bigint = ["_bigint", "dep:malachite-bigint"]
+_bigint = ["dep:num-traits", "dep:num-integer"]  # Internal feature. User must use num-bigint or malachite-bigint instead.
+
 # Do not enable this feature unless you really need it.
 # CPython didn't intend to use FMA for its math library.
 # This project uses this feature in CI to verify the code doesn't have additional bugs on aarch64-apple-darwin.
@@ -17,7 +23,12 @@ mul_add = []
 
 [dependencies]
 libc = "0.2"
+num-complex = { version = "0.4", optional = true }
+num-bigint = { version = "0.4", optional = true }
+num-traits = { version = "0.2", optional = true }
+num-integer = { version = "0.1", optional = true }
+malachite-bigint = { version = "0.2", optional = true }
 
 [dev-dependencies]
 proptest = "1.6.0"
-pyo3 = { version = "0.24", features = ["abi3"] }
+pyo3 = { version = "0.27", features = ["abi3", "auto-initialize"] }

README.md

@@ -11,10 +11,16 @@ A binary representation compatible Rust implementation of Python's math library.
 
 Each function has been carefully translated from CPython's C implementation to Rust, preserving the same algorithms, constants, and corner case handling. The code maintains the same numerical properties, but in Rust!
 
+## Module Structure
+
+- `pymath::math` - Real number math functions (Python's `math` module)
+- `pymath::cmath` - Complex number functions (Python's `cmath` module, requires `complex` feature)
+- `pymath::m` - Direct libm bindings
+
 ## Usage
 
 ```rust
-use pymath::{gamma, lgamma};
+use pymath::math::{gamma, lgamma};
 
 fn main() {
     // Get the same results as Python's math.gamma and math.lgamma

proptest-regressions/gamma.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 6f6becc96f663a83d20def559f516c7c5ce1a90b87c373d6c025dd3ab8f1fc39 # shrinks to x = 5.868849392888587e-309, y = 1.985586796867676e-308

proptest-regressions/math.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 f58d401e380fa9d3a8dd5b137a668be6bd437776f47186a9479aee07c0aa28b8 # shrinks to p1 = 0.0, p2 = 0.0, q1 = -1.156587418587806e301, q2 = -1.315804087909368e-150

proptest-regressions/math/aggregate.txt

@@ -0,0 +1,10 @@
+# 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 c756e015f8e09d21d37e7f805e54f9270cac137e4570d10ecd5302d752e2559f # shrinks to x = -2.9225692145627912e-248, y = 1.565985096226121e-308
+cc 2f4b02dd2c2019dee35dcf980b7b4f077899289a7005bba21f849776768eda7a # shrinks to x = 6.095938323843682e143
+cc b52b190a89a473a7e5269f9a3efa83735318ae3bd52f3340c32bc4a40f4cda88 # shrinks to x = -0.0, y = -9.787367203123051e54
+cc a34a3387ec2547faa88352529f9730d47c6202f9b418a0b4474c9ebb850081e7 # shrinks to x = -1.3916042894622981e-207

proptest-regressions/math/exponential.txt

@@ -4,4 +4,4 @@
 #
 # 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
+cc ac5c11a6ec8450aef2b03644e8b46a066f81c689b55f9657fcd0ae73d5829bc2 # shrinks to x = -1.8815128643365582

proptest-regressions/lib.txt proptest-regressions/math/gamma.txtproptest-regressions/lib.txt renamed to proptest-regressions/math/gamma.txt

@@ -20,3 +20,78 @@ impl TryFrom<c_int> for Error {
         }
     }
 }
+
+/// Set errno to the given value.
+#[inline]
+pub(crate) fn set_errno(value: i32) {
+    unsafe {
+        #[cfg(target_os = "linux")]
+        {
+            *libc::__errno_location() = value;
+        }
+        #[cfg(target_os = "macos")]
+        {
+            *libc::__error() = value;
+        }
+        #[cfg(target_os = "windows")]
+        {
+            unsafe extern "C" {
+                safe fn _errno() -> *mut i32;
+            }
+            *_errno() = value;
+        }
+        #[cfg(all(unix, not(any(target_os = "linux", target_os = "macos"))))]
+        {
+            // FreeBSD, NetBSD, OpenBSD, etc. use __error()
+            *libc::__error() = value;
+        }
+    }
+}
+
+/// Get the current errno value.
+#[inline]
+pub(crate) fn get_errno() -> i32 {
+    unsafe {
+        #[cfg(target_os = "linux")]
+        {
+            *libc::__errno_location()
+        }
+        #[cfg(target_os = "macos")]
+        {
+            *libc::__error()
+        }
+        #[cfg(target_os = "windows")]
+        {
+            unsafe extern "C" {
+                safe fn _errno() -> *mut i32;
+            }
+            *_errno()
+        }
+        #[cfg(all(unix, not(any(target_os = "linux", target_os = "macos"))))]
+        {
+            // FreeBSD, NetBSD, OpenBSD, etc. use __error()
+            *libc::__error()
+        }
+    }
+}
+
+/// Check errno after libm call and convert to Result.
+#[inline]
+pub(crate) fn is_error(x: f64) -> Result<f64> {
+    match get_errno() {
+        0 => Ok(x),
+        libc::EDOM => Err(Error::EDOM),
+        libc::ERANGE => {
+            // Underflow to zero is not an error.
+            // Use 1.5 threshold to handle subnormal results that don't underflow to zero
+            // (e.g., on Ubuntu/ia64) and to correctly detect underflows in expm1()
+            // which may underflow toward -1.0 rather than 0.0. (bpo-46018)
+            if x.abs() < 1.5 {
+                Ok(x)
+            } else {
+                Err(Error::ERANGE)
+            }
+        }
+        _ => Ok(x), // Unknown errno, just return the value
+    }
+}

src/err.rs

@@ -1,29 +1,26 @@
+// Public modules
+pub mod math;
+
+// Internal modules
 mod err;
-mod gamma;
-mod m;
+pub(crate) mod m;
+mod m_sys;
 #[cfg(test)]
 mod test;
 
+// Re-export error types at root level
 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) }
-        }
-    };
+/// Fused multiply-add operation.
+/// When `mul_add` feature is enabled, uses hardware FMA instruction.
+/// Otherwise, uses separate multiply and add operations.
+#[inline(always)]
+pub(crate) fn mul_add(a: f64, b: f64, c: f64) -> f64 {
+    if cfg!(feature = "mul_add") {
+        a.mul_add(b, c)
+    } else {
+        a * b + c
+    }
 }
 
 macro_rules! pyo3_proptest {
@@ -34,8 +31,7 @@ macro_rules! pyo3_proptest {
 
             let rs_result = $fn_name(x);
 
-            pyo3::prepare_freethreaded_python();
-            Python::with_gil(|py| {
+            pyo3::Python::attach(|py| {
                 let math = PyModule::import(py, "math").unwrap();
                 let py_func = math
                     .getattr(stringify!($fn_name))
@@ -60,8 +56,7 @@ macro_rules! pyo3_proptest {
 
             let rs_result = Ok($fn_name(x));
 
-            pyo3::prepare_freethreaded_python();
-            Python::with_gil(|py| {
+            pyo3::Python::attach(|py| {
                 let math = PyModule::import(py, "math").unwrap();
                 let py_func = math
                     .getattr(stringify!($fn_name))
@@ -101,20 +96,4 @@ macro_rules! pyo3_proptest {
     };
 }
 
-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;