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/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>) {

src/math/bigint.rs

@@ -63,7 +63,14 @@ pub fn comb_bigint(n: &BigInt, k: u64) -> BigUint {
 /// - mantissa is in [0.5, 1.0) for positive n
 /// - n ~= mantissa * 2^exponent
 ///
-/// See: _PyLong_Frexp in CPython longobject.c
+/// `_PyLong_Frexp` extracts digits one-by-one into a fixed-size
+/// accumulator and applies a `half_even_correction` lookup table for
+/// rounding.  We instead extract the top 55 bits via a single right
+/// shift and use a sticky-bit to mark whether any discarded bits were
+/// non-zero, then delegate to `BigInt::to_f64()` which performs
+/// IEEE 754 round-half-to-even.  The two approaches are equivalent
+/// because the sticky bit preserves the same rounding information
+/// that the digit-by-digit extraction would.
 fn frexp_bigint(n: &BigInt) -> (f64, i64) {
     let bits = n.bits();
     if bits == 0 {
@@ -87,8 +94,15 @@ fn frexp_bigint(n: &BigInt) -> (f64, i64) {
 
     // Sticky bit: if any shifted-out bits were non-zero, set the LSB.
     // This ensures correct IEEE round-half-to-even when converting to f64.
-    // See _PyLong_Frexp in longobject.c.
-    if (&mantissa_int << shift as u64) != *n {
+    //
+    // `_PyLong_Frexp` checks the remainder from `v_rshift` first, then
+    // iterates shifted-out digits top-down.  We use `trailing_zeros()`
+    // which scans digits bottom-up instead.  The worst-case traversal
+    // order differs (e.g. exact powers of two), but for typical inputs
+    // both terminate in O(1).  If you observe a performance regression
+    // from this, please file a bug report.
+    let tz = n.magnitude().trailing_zeros().unwrap(); // n != 0 here
+    if tz < shift as u64 {
         mantissa_int |= BigInt::from(1);
     }
 

src/math/exponential.rs

@@ -113,7 +113,10 @@ pub fn log(x: f64, base: Option<f64>) -> Result<f64> {
             if den.is_infinite() && b.is_finite() {
                 return Err(crate::Error::EDOM);
             }
-            // log(x, 1) -> division by zero
+            // log(x, 1) -> division by zero.
+            // CPython raises ZeroDivisionError here (via PyNumber_TrueDivide),
+            // but we return EDOM since our error type has no ZeroDivisionError
+            // variant. The caller (e.g. RustPython) may remap this if needed.
             if den == 0.0 {
                 return Err(crate::Error::EDOM);
             }

src/math/misc.rs

@@ -10,6 +10,12 @@ super::libm_simple!(@1 ceil, floor, trunc);
 /// manipulation on the IEEE 754 representation. Steps that overshoot y
 /// are clamped so the result never passes y.
 ///
+/// CPython's math_nextafter_impl accepts a Python integer for steps,
+/// rejects negative values, and saturates overflows to UINT64_MAX. This
+/// Rust API takes `Option<u64>`, so negative rejection and big-int
+/// saturation are structurally unnecessary. The caller (e.g. RustPython)
+/// should handle Python int conversion and negative checks before calling.
+///
 /// See math_nextafter_impl in mathmodule.c.
 #[inline]
 pub fn nextafter(x: f64, y: f64, steps: Option<u64>) -> f64 {
@@ -219,6 +225,12 @@ pub fn fmod(x: f64, y: f64) -> Result<f64> {
 }
 
 /// Return the IEEE 754-style remainder of x with respect to y.
+///
+/// CPython implements this from scratch using fmod (m_remainder in
+/// mathmodule.c) rather than calling the C library's remainder().
+/// We delegate to libm's remainder() which is correct on all platforms
+/// where it conforms to IEEE 754. If you find a platform where the
+/// results differ from CPython, please file a bug.
 #[inline]
 pub fn remainder(x: f64, y: f64) -> Result<f64> {
     super::math_2(x, y, crate::m::remainder)