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,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 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

proptest-regressions/lib.txt

@@ -20,3 +20,65 @@ 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;
+        }
+    }
+}
+
+/// 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()
+        }
+    }
+}
+
+/// 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
+            if x.abs() < 1.5 {
+                Ok(x)
+            } else {
+                Err(Error::ERANGE)
+            }
+        }
+        _ => Ok(x), // Unknown errno, just return the value
+    }
+}

proptest-regressions/math/exponential.txt

@@ -1,29 +1,25 @@
+// Public modules
+pub mod math;
+
+// Internal modules
 mod err;
-mod gamma;
 mod m;
 #[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 +30,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 +55,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 +95,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;

src/err.rs

@@ -4,35 +4,69 @@
 // or by `compiler-builtins` on unsupported platforms.
 #[allow(dead_code)]
 unsafe extern "C" {
+    // Trigonometric functions
     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 cos(n: f64) -> f64;
+    pub fn sin(n: f64) -> f64;
+    pub fn tan(n: f64) -> f64;
+
+    // Hyperbolic functions
+    pub fn acosh(n: f64) -> f64;
+    pub fn asinh(n: f64) -> f64;
+    pub fn atanh(n: f64) -> f64;
     pub fn cosh(n: f64) -> f64;
+    pub fn sinh(n: f64) -> f64;
+    pub fn tanh(n: f64) -> f64;
+
+    // Exponential and logarithmic functions
+    pub fn exp(n: f64) -> f64;
+    pub fn exp2(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;
+    pub fn log(n: f64) -> f64;
+    pub fn log10(n: f64) -> f64;
+    pub fn log1p(n: f64) -> f64;
+    pub fn log1pf(n: f32) -> f32;
+    pub fn log2(n: f64) -> f64;
+
+    // Power functions
+    pub fn cbrt(n: f64) -> f64;
+    pub fn cbrtf(n: 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 pow(x: f64, y: f64) -> f64;
+    pub fn sqrt(n: f64) -> f64;
+
+    // Floating-point manipulation functions
+    pub fn ceil(n: f64) -> f64;
+    pub fn copysign(x: f64, y: f64) -> f64;
+    pub fn fabs(n: f64) -> f64;
+    pub fn fdim(a: f64, b: f64) -> f64;
+    pub fn fdimf(a: f32, b: f32) -> f32;
+    pub fn floor(n: f64) -> f64;
+    pub fn fmod(x: f64, y: f64) -> f64;
+    pub fn frexp(n: f64, exp: *mut i32) -> f64;
+    pub fn ldexp(x: f64, n: i32) -> f64;
+    pub fn modf(n: f64, iptr: *mut f64) -> f64;
+    pub fn nextafter(x: f64, y: f64) -> f64;
+    pub fn remainder(x: f64, y: f64) -> f64;
+    pub fn trunc(n: f64) -> f64;
+
+    // Special functions
     pub fn erf(n: f64) -> f64;
-    pub fn erff(n: f32) -> f32;
     pub fn erfc(n: f64) -> f64;
+    pub fn erff(n: f32) -> f32;
     pub fn erfcf(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 tgamma(n: f64) -> f64;
+    pub fn tgammaf(n: f32) -> f32;
 
     // pub fn acosf128(n: f128) -> f128;
     // pub fn asinf128(n: f128) -> f128;