Align behaviors to cpython (#18)

* O(1) bit-level nextafter with steps parameter

Replace the O(n) loop in nextafter(x, y, steps) with O(1) IEEE 754
bit manipulation matching math_nextafter_impl in mathmodule.c.
Handles sign-crossing, saturation, and all edge cases.

Add pyo3 proptest and edge/extreme-step tests for nextafter.

* Fix cmath.log error propagation and bigint log precision

cmath: Rewrite log(z, base) to match cmath_log_impl errno
semantics. c_log always returns a value and separately reports
EDOM for zero. c_quot returns EDOM and (0,0) for zero
denominator instead of NaN.

bigint: Add sticky bit to frexp_bigint for correct IEEE
round-half-to-even, matching _PyLong_Frexp. Use mul_add for
log/log2/log10 bigint to match loghelper fma usage.

Add regression tests for remainder tie-to-even and signed zero.

* Document design decisions and improve sumprod

- math_1/math_2: document why errno handling differs from
  CPython (platform-specific unreliability, output checks
  sufficient, verified by proptest)
- math.log: document EDOM substitution for ZeroDivisionError
- math.remainder: document libm delegation rationale
- sumprod: return Result for length mismatch instead of panic,
  improve overflow fallback to continue from where the fast
  path stopped instead of restarting from scratch

Author: youknowone

Cargo.toml

@@ -3,7 +3,7 @@ name = "pymath"
 authors = ["Jeong, YunWon <jeong@youknowone.org>"]
 repository = "https://github.com/RustPython/pymath"
 description = "A binary representation compatible Rust implementation of Python's math library."
-version = "0.1.5"
+version = "0.2.0"
 edition = "2024"
 license = "PSF-2.0"
 

src/cmath/exponential.rs

@@ -177,40 +177,47 @@ pub(crate) fn ln(z: Complex64) -> Result<Complex64> {
 /// If base is Some(b), returns log(z) / log(b).
 #[inline]
 pub fn log(z: Complex64, base: Option<Complex64>) -> Result<Complex64> {
-    // c_log always returns a value, but sets errno for special cases.
-    // The error check happens at the end of cmath_log_impl.
-    // For log(z) without base: z=0 raises EDOM
-    // For log(z, base): z=0 doesn't raise because c_log(base) clears errno
-    let z_is_zero = z.re == 0.0 && z.im == 0.0;
+    let (log_z, mut err) = c_log(z);
     match base {
-        None => {
-            // No base: raise error if z=0
-            if z_is_zero {
-                return Err(Error::EDOM);
-            }
-            ln(z)
-        }
+        None => err.map_or(Ok(log_z), Err),
         Some(b) => {
-            // With base: z=0 is allowed (second ln clears the "errno")
-            let log_z = ln(z)?;
-            let log_b = ln(b)?;
-            // Use _Py_c_quot-style division to preserve sign of zero
-            Ok(c_quot(log_z, log_b))
+            // Like cmath_log_impl, the second c_log call overwrites
+            // any pending error from the first one.
+            let (log_b, base_err) = c_log(b);
+            err = base_err;
+            let (q, quot_err) = c_quot(log_z, log_b);
+            if let Some(e) = quot_err {
+                err = Some(e);
+            }
+            err.map_or(Ok(q), Err)
         }
     }
 }
 
+/// c_log behavior: always returns a value, but reports EDOM for zero.
+#[inline]
+fn c_log(z: Complex64) -> (Complex64, Option<Error>) {
+    let r = ln(z).expect("ln handles special values without failing");
+    if z.re == 0.0 && z.im == 0.0 {
+        (r, Some(Error::EDOM))
+    } else {
+        (r, None)
+    }
+}
+
 /// Complex division following _Py_c_quot algorithm.
 /// This preserves the sign of zero correctly and recovers infinities
 /// from NaN results per C11 Annex G.5.2.
 #[inline]
-fn c_quot(a: Complex64, b: Complex64) -> Complex64 {
+fn c_quot(a: Complex64, b: Complex64) -> (Complex64, Option<Error>) {
     let abs_breal = m::fabs(b.re);
     let abs_bimag = m::fabs(b.im);
+    let mut err = None;
 
     let mut r = if abs_breal >= abs_bimag {
         if abs_breal == 0.0 {
-            Complex64::new(f64::NAN, f64::NAN)
+            err = Some(Error::EDOM);
+            Complex64::new(0.0, 0.0)
         } else {
             let ratio = b.im / b.re;
             let denom = b.re + b.im * ratio;
@@ -244,7 +251,7 @@ fn c_quot(a: Complex64, b: Complex64) -> Complex64 {
         }
     }
 
-    r
+    (r, err)
 }
 
 /// Complex base-10 logarithm.
@@ -325,6 +332,53 @@ mod tests {
         });
     }
 
+    fn test_log_error(z: Complex64, base: Complex64) {
+        use pyo3::prelude::*;
+
+        let rs_result = log(z, Some(base));
+
+        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();
+                    panic!(
+                        "log({}+{}j, {}+{}j): expected ValueError, got ({}, {})",
+                        z.re,
+                        z.im,
+                        base.re,
+                        base.im,
+                        c.real(),
+                        c.imag()
+                    );
+                }
+                Err(err) => {
+                    assert!(
+                        err.is_instance_of::<pyo3::exceptions::PyValueError>(py),
+                        "log({}+{}j, {}+{}j): expected ValueError, got {err:?}",
+                        z.re,
+                        z.im,
+                        base.re,
+                        base.im,
+                    );
+                    assert!(
+                        matches!(rs_result, Err(crate::Error::EDOM)),
+                        "log({}+{}j, {}+{}j): expected Err(EDOM), got {rs_result:?}",
+                        z.re,
+                        z.im,
+                        base.re,
+                        base.im,
+                    );
+                }
+            }
+        });
+    }
+
     use crate::test::EDGE_VALUES;
 
     #[test]
@@ -382,6 +436,29 @@ mod tests {
         }
     }
 
+    #[test]
+    fn regression_c_quot_zero_denominator_sets_edom() {
+        let (q, err) = c_quot(Complex64::new(2.0, -3.0), Complex64::new(0.0, 0.0));
+        assert_eq!(err, Some(crate::Error::EDOM));
+        assert_eq!(q.re.to_bits(), 0.0f64.to_bits());
+        assert_eq!(q.im.to_bits(), 0.0f64.to_bits());
+    }
+
+    #[test]
+    fn regression_log_zero_quotient_denominator_raises_edom() {
+        let cases = [
+            (Complex64::new(2.0, 0.0), Complex64::new(1.0, 0.0)),
+            (Complex64::new(1.0, 0.0), Complex64::new(1.0, 0.0)),
+            (Complex64::new(2.0, 0.0), Complex64::new(0.0, 0.0)),
+            (Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)),
+            (Complex64::new(0.0, 0.0), Complex64::new(0.0, 0.0)),
+        ];
+
+        for (z, base) in cases {
+            test_log_error(z, base);
+        }
+    }
+
     proptest::proptest! {
         #[test]
         fn proptest_sqrt(re: f64, im: f64) {

src/math.rs

@@ -49,10 +49,26 @@ macro_rules! libm_simple {
 
 pub(crate) use libm_simple;
 
-/// math_1: wrapper for 1-arg functions
+/// Wrapper for 1-arg libm functions, corresponding to FUNC1/is_error in
+/// mathmodule.c.
+///
 /// - isnan(r) && !isnan(x) -> domain error
 /// - isinf(r) && isfinite(x) -> overflow (can_overflow=true) or domain error (can_overflow=false)
 /// - isfinite(r) && errno -> check errno (unnecessary on most platforms)
+///
+/// CPython's approach: clear errno, call libm, then inspect both the result
+/// and errno to classify errors. We rely primarily on output inspection
+/// (NaN/Inf checks) because:
+///
+/// - On macOS and Windows, libm functions do not reliably set errno for
+///   edge cases, so CPython's own is_error() skips the errno check there
+///   too (it only uses it as a fallback on other Unixes).
+/// - The NaN/Inf output checks are sufficient to detect all domain and
+///   range errors on every platform we test against (verified by proptest
+///   and edgetest against CPython via pyo3).
+/// - The errno-only branch (finite result with errno set) is kept for
+///   non-macOS/non-Windows Unixes where libm might signal an error
+///   without producing a NaN/Inf result.
 #[inline]
 pub(crate) fn math_1(x: f64, func: fn(f64) -> f64, can_overflow: bool) -> crate::Result<f64> {
     crate::err::set_errno(0);
@@ -75,9 +91,17 @@ pub(crate) fn math_1(x: f64, func: fn(f64) -> f64, can_overflow: bool) -> crate:
     Ok(r)
 }
 
-/// math_2: wrapper for 2-arg functions
+/// Wrapper for 2-arg libm functions, corresponding to FUNC2 in
+/// mathmodule.c.
+///
 /// - isnan(r) && !isnan(x) && !isnan(y) -> domain error
 /// - isinf(r) && isfinite(x) && isfinite(y) -> range error
+///
+/// Unlike math_1, this does not set/check errno at all. CPython's FUNC2
+/// does clear and check errno, but the NaN/Inf output checks already
+/// cover all error cases for the 2-arg functions we wrap (atan2, fmod,
+/// copysign, remainder, pow). This is verified by bit-exact proptest
+/// and edgetest against CPython.
 #[inline]
 pub(crate) fn math_2(x: f64, y: f64, func: fn(f64, f64) -> f64) -> crate::Result<f64> {
     let r = func(x, y);

src/math/aggregate.rs

@@ -231,6 +231,11 @@ pub fn vector_norm(vec: &[f64], max: f64, found_nan: bool) -> f64 {
 ///
 /// The points are given as sequences of coordinates.
 /// Uses high-precision vector_norm algorithm.
+///
+/// Panics if `p` and `q` have different lengths. CPython raises ValueError
+/// for mismatched dimensions, but in this Rust API the caller is expected
+/// to guarantee equal-length slices. A length mismatch is a programming
+/// error, not a runtime condition.
 pub fn dist(p: &[f64], q: &[f64]) -> f64 {
     assert_eq!(
         p.len(),
@@ -261,24 +266,52 @@ pub fn dist(p: &[f64], q: &[f64]) -> f64 {
 
 /// Return the sum of products of values from two sequences (float version).
 ///
-/// Uses TripleLength arithmetic for high precision.
-/// Equivalent to sum(p[i] * q[i] for i in range(len(p))).
-pub fn sumprod(p: &[f64], q: &[f64]) -> f64 {
-    assert_eq!(p.len(), q.len(), "Inputs are not the same length");
+/// Uses TripleLength arithmetic for the fast path, then falls back to
+/// ordinary floating-point multiply/add starting at the first unsupported
+/// pair, matching Python's staged `math.sumprod` behavior for float inputs.
+///
+/// CPython's math_sumprod_impl is a 3-stage state machine that handles
+/// int/float/generic Python objects. This function only covers the float
+/// path (`&[f64]`). The int accumulation and generic PyNumber fallback
+/// stages are Python type-system concerns and should be handled by the
+/// caller (e.g. RustPython) before delegating here.
+///
+/// Returns EDOM if the inputs are not the same length.
+pub fn sumprod(p: &[f64], q: &[f64]) -> crate::Result<f64> {
+    if p.len() != q.len() {
+        return Err(crate::Error::EDOM);
+    }
 
+    let mut total = 0.0;
     let mut flt_total = TL_ZERO;
+    let mut flt_path_enabled = true;
+    let mut i = 0;
 
-    for (&pi, &qi) in p.iter().zip(q.iter()) {
-        let new_flt_total = tl_fma(pi, qi, flt_total);
-        if new_flt_total.hi.is_finite() {
-            flt_total = new_flt_total;
-        } else {
-            // Overflow or special value, fall back to simple sum
-            return p.iter().zip(q.iter()).map(|(a, b)| a * b).sum();
+    while i < p.len() {
+        let pi = p[i];
+        let qi = q[i];
+
+        if flt_path_enabled {
+            let new_flt_total = tl_fma(pi, qi, flt_total);
+            if new_flt_total.hi.is_finite() {
+                flt_total = new_flt_total;
+                i += 1;
+                continue;
+            }
+
+            flt_path_enabled = false;
+            total += tl_to_d(flt_total);
         }
+
+        total += pi * qi;
+        i += 1;
     }
 
-    tl_to_d(flt_total)
+    Ok(if flt_path_enabled {
+        tl_to_d(flt_total)
+    } else {
+        total
+    })
 }
 
 /// Return the sum of products of values from two sequences (integer version).
@@ -427,14 +460,27 @@ mod tests {
         crate::test::with_py_math(|py, math| {
             let py_p = pyo3::types::PyList::new(py, p).unwrap();
             let py_q = pyo3::types::PyList::new(py, q).unwrap();
-            let py: f64 = math
-                .getattr("sumprod")
-                .unwrap()
-                .call1((py_p, py_q))
-                .unwrap()
-                .extract()
-                .unwrap();
-            crate::test::assert_f64_eq(py, rs, format_args!("sumprod({p:?}, {q:?})"));
+            let py_result = math.getattr("sumprod").unwrap().call1((py_p, py_q));
+            match py_result {
+                Ok(py_val) => {
+                    let py: f64 = py_val.extract().unwrap();
+                    let rs = rs.unwrap_or_else(|e| {
+                        panic!("sumprod({p:?}, {q:?}): py={py} but rs returned error {e:?}")
+                    });
+                    crate::test::assert_f64_eq(py, rs, format_args!("sumprod({p:?}, {q:?})"));
+                }
+                Err(e) => {
+                    if e.is_instance_of::<pyo3::exceptions::PyValueError>(py) {
+                        assert_eq!(
+                            rs.as_ref().err(),
+                            Some(&crate::Error::EDOM),
+                            "sumprod({p:?}, {q:?}): py raised ValueError but rs={rs:?}"
+                        );
+                    } else {
+                        panic!("sumprod({p:?}, {q:?}): py raised unexpected error {e}");
+                    }
+                }
+            }
         });
     }
 
@@ -444,6 +490,9 @@ mod tests {
         test_sumprod_impl(&[], &[]);
         test_sumprod_impl(&[1.0], &[2.0]);
         test_sumprod_impl(&[1e100, 1e100], &[1e100, -1e100]);
+        test_sumprod_impl(&[1.0, 1e308, -1e308], &[1.0, 2.0, 2.0]);
+        test_sumprod_impl(&[1e-16, 1e308, -1e308], &[1.0, 2.0, 2.0]);
+        test_sumprod_impl(&[1.0], &[]);
     }
 
     fn test_prod_impl(values: &[f64], start: Option<f64>) {