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.

Author: youknowone

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/bigint.rs

@@ -81,10 +81,17 @@ fn frexp_bigint(n: &BigInt) -> (f64, i64) {
         return (m, e as i64);
     }
 
-    // For large integers, extract top ~53 bits
-    // Shift right to keep DBL_MANT_DIG + 2 = 55 bits for rounding
+    // For large integers, extract top DBL_MANT_DIG + 2 = 55 bits for rounding
     let shift = bits - 55;
-    let mantissa_int = n >> shift as u64;
+    let mut mantissa_int = n >> shift as u64;
+
+    // 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 {
+        mantissa_int |= BigInt::from(1);
+    }
+
     let mut x = mantissa_int.to_f64().unwrap();
 
     // x is now approximately n / 2^shift, with ~55 bits of precision
@@ -119,7 +126,7 @@ pub fn log_bigint(n: &BigInt, base: Option<f64>) -> crate::Result<f64> {
     // Use frexp decomposition for large values
     // n ~= x * 2^e, so log(n) = log(x) + log(2) * e
     let (x, e) = frexp_bigint(n);
-    let log_n = crate::m::log(x) + std::f64::consts::LN_2 * (e as f64);
+    let log_n = crate::mul_add(crate::m::log(2.0), e as f64, crate::m::log(x));
 
     match base {
         None => Ok(log_n),
@@ -150,7 +157,11 @@ pub fn log2_bigint(n: &BigInt) -> crate::Result<f64> {
     // Use frexp decomposition for large values
     // n ~= x * 2^e, so log2(n) = log2(x) + e
     let (x, e) = frexp_bigint(n);
-    Ok(crate::m::log2(x) + (e as f64))
+    Ok(crate::mul_add(
+        crate::m::log2(2.0),
+        e as f64,
+        crate::m::log2(x),
+    ))
 }
 
 /// Return the base-10 logarithm of a BigInt.
@@ -171,7 +182,11 @@ pub fn log10_bigint(n: &BigInt) -> crate::Result<f64> {
     // Use frexp decomposition for large values
     // n ~= x * 2^e, so log10(n) = log10(x) + log10(2) * e
     let (x, e) = frexp_bigint(n);
-    Ok(crate::m::log10(x) + std::f64::consts::LOG10_2 * (e as f64))
+    Ok(crate::mul_add(
+        crate::m::log10(2.0),
+        e as f64,
+        crate::m::log10(x),
+    ))
 }
 
 /// Compute ldexp(x, exp) where exp is a BigInt.
@@ -333,6 +348,33 @@ mod tests {
         });
     }
 
+    fn assert_exact_log_bigint_bits(n: &BigInt, func_name: &str, rs: crate::Result<f64>) {
+        crate::test::with_py_math(|py, math| {
+            let n_str = n.to_string();
+            let builtins = pyo3::types::PyModule::import(py, "builtins").unwrap();
+            let py_n = builtins
+                .getattr("int")
+                .unwrap()
+                .call1((n_str.as_str(),))
+                .unwrap();
+            let py_f: f64 = math
+                .getattr(func_name)
+                .unwrap()
+                .call1((py_n,))
+                .unwrap()
+                .extract()
+                .unwrap();
+            let rs_f = rs.unwrap();
+            assert_eq!(
+                py_f.to_bits(),
+                rs_f.to_bits(),
+                "{func_name}({n}): py={py_f} ({:#x}) vs rs={rs_f} ({:#x})",
+                py_f.to_bits(),
+                rs_f.to_bits()
+            );
+        });
+    }
+
     #[test]
     fn edgetest_log_bigint() {
         // Small values
@@ -410,6 +452,45 @@ mod tests {
         }
     }
 
+    #[test]
+    fn regression_log_bigint_pow10_309() {
+        let n = BigInt::from(10).pow(309);
+        assert_exact_log_bigint_bits(&n, "log", log_bigint(&n, None));
+    }
+
+    #[test]
+    fn regression_log2_bigint_pow10_309() {
+        let n = BigInt::from(10).pow(309);
+        assert_exact_log_bigint_bits(&n, "log2", log2_bigint(&n));
+    }
+
+    #[test]
+    fn regression_log10_bigint_pow10_309() {
+        let n = BigInt::from(10).pow(309);
+        assert_exact_log_bigint_bits(&n, "log10", log10_bigint(&n));
+    }
+
+    #[test]
+    fn regression_log_bigint_sticky_bit_rounding() {
+        // If frexp_bigint drops the sticky bit, this rounds 1 ULP low.
+        let n = (BigInt::from(0x40000000000c02u64) << 970u32) + BigInt::from(1u8);
+        assert_exact_log_bigint_bits(&n, "log", log_bigint(&n, None));
+    }
+
+    #[test]
+    fn regression_log2_bigint_sticky_bit_rounding() {
+        // If frexp_bigint drops the sticky bit, this rounds 1 ULP low.
+        let n = (BigInt::from(0x400000000010a2u64) << 970u32) + BigInt::from(1u8);
+        assert_exact_log_bigint_bits(&n, "log2", log2_bigint(&n));
+    }
+
+    #[test]
+    fn regression_log10_bigint_sticky_bit_rounding() {
+        // If frexp_bigint drops the sticky bit, this rounds 1 ULP low.
+        let n = (BigInt::from(0x4000000000049au64) << 970u32) + BigInt::from(1u8);
+        assert_exact_log_bigint_bits(&n, "log10", log10_bigint(&n));
+    }
+
     // ldexp_bigint tests
 
     fn test_ldexp_bigint_impl(x: f64, exp: &BigInt) {

src/math/integer.rs

@@ -730,8 +730,8 @@ mod tests {
         128, // Table boundary + 1
         // Powers of 2 and boundaries
         64,
-        63,   // 2^6 - 1
-        65,   // 2^6 + 1
+        63, // 2^6 - 1
+        65, // 2^6 + 1
         1024,
         65535, // 2^16 - 1
         65536, // 2^16

src/math/misc.rs

@@ -390,6 +390,64 @@ mod tests {
         }
     }
 
+    #[test]
+    fn regression_remainder_halfway_even_cases() {
+        // These cases exercise the half-way branch in CPython's custom
+        // m_remainder implementation, where ties are resolved toward the
+        // even multiple of y.
+        let cases = [
+            ((6.0, 4.0), 0xc000_0000_0000_0000_u64),  // -2.0
+            ((3.0, 2.0), 0xbff0_0000_0000_0000_u64),  // -1.0
+            ((5.0, 2.0), 0x3ff0_0000_0000_0000_u64),  // 1.0
+            ((-6.0, 4.0), 0x4000_0000_0000_0000_u64), // 2.0
+            ((-5.0, 2.0), 0xbff0_0000_0000_0000_u64), // -1.0
+            ((1.5, 1.0), 0xbfe0_0000_0000_0000_u64),  // -0.5
+            ((2.5, 1.0), 0x3fe0_0000_0000_0000_u64),  // 0.5
+            ((3.5, 1.0), 0xbfe0_0000_0000_0000_u64),  // -0.5
+            ((4.5, 1.0), 0x3fe0_0000_0000_0000_u64),  // 0.5
+            ((5.5, 1.0), 0xbfe0_0000_0000_0000_u64),  // -0.5
+        ];
+
+        for &((x, y), expected_bits) in &cases {
+            let r = remainder(x, y).unwrap();
+            assert_eq!(
+                r.to_bits(),
+                expected_bits,
+                "remainder({x}, {y}) = {r} ({:#x}), expected {:#x}",
+                r.to_bits(),
+                expected_bits
+            );
+            test_remainder(x, y);
+        }
+    }
+
+    #[test]
+    fn regression_remainder_boundary_and_signed_zero_cases() {
+        // These cases pin the sign of zero and the behavior immediately
+        // around the half-way boundary.
+        let just_below_half = f64::from_bits(0x3fdf_ffff_ffff_ffff);
+        let just_above_half = f64::from_bits(0x3fe0_0000_0000_0001);
+        let cases = [
+            ((4.0, 2.0), 0x0000_0000_0000_0000_u64),  // +0.0
+            ((-4.0, 2.0), 0x8000_0000_0000_0000_u64), // -0.0
+            ((4.0, -2.0), 0x0000_0000_0000_0000_u64), // +0.0
+            ((just_below_half, 1.0), 0x3fdf_ffff_ffff_ffff_u64),
+            ((just_above_half, 1.0), 0xbfdf_ffff_ffff_fffe_u64),
+        ];
+
+        for &((x, y), expected_bits) in &cases {
+            let r = remainder(x, y).unwrap();
+            assert_eq!(
+                r.to_bits(),
+                expected_bits,
+                "remainder({x}, {y}) = {r} ({:#x}), expected {:#x}",
+                r.to_bits(),
+                expected_bits
+            );
+            test_remainder(x, y);
+        }
+    }
+
     #[test]
     fn edgetest_copysign() {
         for &x in crate::test::EDGE_VALUES {