fix a lot

Author: youknowone

Cargo.toml

@@ -3,7 +3,7 @@ name = "pymath"
 author = "Jeong, YunWon <[email protected]>"
 repository = "https://github.com/RustPython/pymath"
 description = "A binary representation compatible Rust implementation of Python's math library."
-version = "0.1.1"
+version = "0.1.2"
 edition = "2024"
 license = "PSF-2.0"
 

README.md

@@ -57,13 +57,13 @@ Every function produces identical results to Python at the binary representation
 - [x] `nextafter`
 - [x] `pi`
 - [x] `pow`
-- [x] `prod`
+- [x] `prod` (`prod`, `prod_int`)
 - [x] `radians`
 - [x] `remainder`
 - [x] `sin`
 - [x] `sinh`
 - [x] `sqrt`
-- [x] `sumprod` (`sumprod`, `sumprod_int`, `sumprod_generic`)
+- [x] `sumprod` (`sumprod`, `sumprod_int`)
 - [x] `tan`
 - [x] `tanh`
 - [x] `tau`

proptest-regressions/cmath/trigonometric.txt

@@ -9,3 +9,4 @@ cc 409f359905a7e77eacdbd4643c3a1483cc128e46a1c06795c246f3d57a2e61b9 # shrinks to
 cc 5acbcf76c3bcbaf0b2b89a987e38c75265434f96fdc11598830e2221df72fc53 # shrinks to re = 0.00010260965539526095, im = 4.984721877290597e-19
 cc ac04191f916633de0408e071f5f6ada8f7fe7a1be02caca7a32f37ecb148a5ac # shrinks to re = 3.049837651806167e74, im = -2.222842222335753e-166
 cc 77ec3c08d2e057fb9467eb13c3b953e4e30f2800259d3653f7bcdfcbaf53614f # shrinks to re = 7.812038268590211e52, im = -2.623972069152808e-109
+cc c8fa2873ca664202a38d40ef3766b0f9f07dc82c1d5cd08dcd6d180564a013a1 # shrinks to re = -0.0016106222874309162, im = 0.10954881424244108

proptest-regressions/math/misc.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 46166a18f59569bb11107fd280e53c8f37902760e15afbe67d6951dcd3032c72 # shrinks to x = 0.0, y = 0.0, z = 0.0

src/cmath/exponential.rs

@@ -131,9 +131,10 @@ pub fn exp(z: Complex64) -> Result<Complex64> {
     Ok(Complex64::new(r_re, r_im))
 }
 
-/// Complex natural logarithm.
+/// Complex natural logarithm (ln, base e).
+/// TODO: consider to expose API
 #[inline]
-pub fn log(z: Complex64) -> Result<Complex64> {
+pub(crate) fn ln(z: Complex64) -> Result<Complex64> {
     special_value!(z, LOG_SPECIAL_VALUES);
 
     let ax = m::fabs(z.re);
@@ -167,10 +168,48 @@ pub fn log(z: Complex64) -> Result<Complex64> {
     Ok(Complex64::new(r_re, r_im))
 }
 
+/// Complex logarithm with optional base.
+///
+/// If base is None, returns the natural logarithm.
+/// If base is Some(b), returns log(z) / log(b).
+#[inline]
+pub fn log(z: Complex64, base: Option<Complex64>) -> Result<Complex64> {
+    let log_z = ln(z)?;
+    match base {
+        None => Ok(log_z),
+        Some(b) => {
+            let log_b = ln(b)?;
+            // Use _Py_c_quot-style division to preserve sign of zero
+            Ok(c_quot(log_z, log_b))
+        }
+    }
+}
+
+/// Complex division following _Py_c_quot algorithm.
+/// This preserves the sign of zero correctly.
+#[inline]
+fn c_quot(a: Complex64, b: Complex64) -> Complex64 {
+    let abs_breal = m::fabs(b.re);
+    let abs_bimag = m::fabs(b.im);
+
+    if abs_breal >= abs_bimag {
+        if abs_breal == 0.0 {
+            return Complex64::new(f64::NAN, f64::NAN);
+        }
+        let ratio = b.im / b.re;
+        let denom = b.re + b.im * ratio;
+        Complex64::new((a.re + a.im * ratio) / denom, (a.im - a.re * ratio) / denom)
+    } else {
+        let ratio = b.re / b.im;
+        let denom = b.re * ratio + b.im;
+        Complex64::new((a.re * ratio + a.im) / denom, (a.im * ratio - a.re) / denom)
+    }
+}
+
 /// Complex base-10 logarithm.
 #[inline]
 pub fn log10(z: Complex64) -> Result<Complex64> {
-    let r = log(z)?;
+    let r = ln(z)?;
     Ok(Complex64::new(r.re / M_LN10, r.im / M_LN10))
 }
 
@@ -191,13 +230,61 @@ mod tests {
     fn test_exp(re: f64, im: f64) {
         test_cmath_func("exp", exp, re, im);
     }
-    fn test_log(re: f64, im: f64) {
-        test_cmath_func("log", log, re, im);
+    fn test_log_without_base(re: f64, im: f64) {
+        test_cmath_func("log", |z| log(z, None), re, im);
     }
     fn test_log10(re: f64, im: f64) {
         test_cmath_func("log10", log10, re, im);
     }
 
+    /// Test log with base - compares with Python's cmath.log(z, base)
+    fn test_log(z_re: f64, z_im: f64, base_re: f64, base_im: f64) {
+        use pyo3::prelude::*;
+
+        let z = Complex64::new(z_re, z_im);
+        let base = Complex64::new(base_re, base_im);
+        let rs_result = log(z, Some(base));
+
+        pyo3::Python::attach(|py| {
+            let cmath = pyo3::types::PyModule::import(py, "cmath").unwrap();
+            let py_z = pyo3::types::PyComplex::from_doubles(py, z_re, z_im);
+            let py_base = pyo3::types::PyComplex::from_doubles(py, base_re, base_im);
+            let py_result = cmath.getattr("log").unwrap().call1((py_z, py_base));
+
+            match py_result {
+                Ok(result) => {
+                    use pyo3::types::PyComplexMethods;
+                    let c = result.cast::<pyo3::types::PyComplex>().unwrap();
+                    let py_re = c.real();
+                    let py_im = c.imag();
+                    match rs_result {
+                        Ok(rs) => {
+                            crate::cmath::tests::assert_complex_eq(
+                                py_re,
+                                py_im,
+                                rs,
+                                "log_with_base",
+                                z_re,
+                                z_im,
+                            );
+                        }
+                        Err(e) => {
+                            panic!(
+                                "log({z_re}+{z_im}j, {base_re}+{base_im}j): py=({py_re}, {py_im}) but rs returned error {e:?}"
+                            );
+                        }
+                    }
+                }
+                Err(_) => {
+                    assert!(
+                        rs_result.is_err(),
+                        "log({z_re}+{z_im}j, {base_re}+{base_im}j): py raised error but rs={rs_result:?}"
+                    );
+                }
+            }
+        });
+    }
+
     use crate::test::EDGE_VALUES;
 
     #[test]
@@ -222,7 +309,7 @@ mod tests {
     fn edgetest_log() {
         for &re in &EDGE_VALUES {
             for &im in &EDGE_VALUES {
-                test_log(re, im);
+                test_log_without_base(re, im);
             }
         }
     }
@@ -236,6 +323,25 @@ mod tests {
         }
     }
 
+    #[test]
+    fn edgetest_log_with_base() {
+        // Test log with various bases - sign preservation edge cases
+        let bases = [0.5, 2.0, 10.0];
+        let values = [0.01, 0.1, 0.5, 1.0, 2.0, 10.0, 100.0];
+        for &base in &bases {
+            for &v in &values {
+                test_log(v, 0.0, base, 0.0);
+            }
+        }
+        // Additional edge cases with imaginary parts
+        for &z_re in &EDGE_VALUES {
+            for &z_im in &EDGE_VALUES {
+                test_log(z_re, z_im, 2.0, 0.0);
+                test_log(z_re, z_im, 0.5, 0.0);
+            }
+        }
+    }
+
     proptest::proptest! {
         #[test]
         fn proptest_sqrt(re: f64, im: f64) {
@@ -249,7 +355,7 @@ mod tests {
 
         #[test]
         fn proptest_log(re: f64, im: f64) {
-            test_log(re, im);
+            test_log_without_base(re, im);
         }
 
         #[test]

src/cmath/trigonometric.rs

@@ -83,8 +83,8 @@ static SINH_SPECIAL_VALUES: [[Complex64; 7]; 7] = [
 static TANH_SPECIAL_VALUES: [[Complex64; 7]; 7] = [
     [c(-1.0, 0.0), c(U, U), c(-1.0, -0.0), c(-1.0, 0.0), c(U, U), c(-1.0, 0.0), c(-1.0, 0.0)],
     [c(N, N),      c(U, U), c(U, U),       c(U, U),      c(U, U), c(N, N),      c(N, N)],
-    [c(-0.0, N),   c(U, U), c(-0.0, -0.0), c(-0.0, 0.0), c(U, U), c(-0.0, N),   c(-0.0, N)],
-    [c(0.0, N),    c(U, U), c(0.0, -0.0),  c(0.0, 0.0),  c(U, U), c(0.0, N),    c(0.0, N)],
+    [c(N, N),      c(U, U), c(-0.0, -0.0), c(-0.0, 0.0), c(U, U), c(N, N),      c(N, N)],
+    [c(N, N),      c(U, U), c(0.0, -0.0),  c(0.0, 0.0),  c(U, U), c(N, N),      c(N, N)],
     [c(N, N),      c(U, U), c(U, U),       c(U, U),      c(U, U), c(N, N),      c(N, N)],
     [c(1.0, 0.0),  c(U, U), c(1.0, -0.0),  c(1.0, 0.0),  c(U, U), c(1.0, 0.0),  c(1.0, 0.0)],
     [c(N, N),      c(N, N), c(N, -0.0),    c(N, 0.0),    c(N, N), c(N, N),      c(N, N)],

src/math.rs

@@ -10,11 +10,11 @@ mod misc;
 mod trigonometric;
 
 // Re-export from submodules
-pub use aggregate::{dist, fsum, prod, sumprod};
+pub use aggregate::{dist, fsum, prod, prod_int, sumprod, sumprod_int};
 pub use exponential::{cbrt, exp, exp2, expm1, log, log1p, log2, log10, pow, sqrt};
 pub use gamma::{erf, erfc, gamma, lgamma};
 pub use misc::{
-    ceil, copysign, fabs, floor, fmod, frexp, isclose, isfinite, isinf, isnan, ldexp, modf,
+    ceil, copysign, fabs, floor, fma, fmod, frexp, isclose, isfinite, isinf, isnan, ldexp, modf,
     nextafter, remainder, trunc, ulp,
 };
 pub use trigonometric::{