Fix hypot,isclose,integer

Author: youknowone

proptest-regressions/math/integer.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 79031914da5204bfc75b0d7cf66e7f76d2e455d6c2837b66cbd11ebf7225be4a # shrinks to n = 32, k = 30

src/cmath/misc.rs

@@ -98,22 +98,38 @@ pub fn abs(z: Complex64) -> f64 {
 }
 
 /// Determine whether two complex numbers are close in value.
+///
+/// Default tolerances: rel_tol = 1e-09, abs_tol = 0.0
+/// Returns Err(EDOM) if rel_tol or abs_tol is negative.
 #[inline]
-pub fn isclose(a: Complex64, b: Complex64, rel_tol: f64, abs_tol: f64) -> bool {
+pub fn isclose(
+    a: Complex64,
+    b: Complex64,
+    rel_tol: Option<f64>,
+    abs_tol: Option<f64>,
+) -> Result<bool> {
+    let rel_tol = rel_tol.unwrap_or(1e-09);
+    let abs_tol = abs_tol.unwrap_or(0.0);
+
+    // Tolerances must be non-negative
+    if rel_tol < 0.0 || abs_tol < 0.0 {
+        return Err(Error::EDOM);
+    }
+
     // short circuit exact equality
     if a.re == b.re && a.im == b.im {
-        return true;
+        return Ok(true);
     }
 
     // This catches the case of two infinities of opposite sign, or
     // one infinity and one finite number.
     if a.re.is_infinite() || a.im.is_infinite() || b.re.is_infinite() || b.im.is_infinite() {
-        return false;
+        return Ok(false);
     }
 
     // now do the regular computation
     let diff = abs(Complex64::new(a.re - b.re, a.im - b.im));
-    (diff <= rel_tol * abs(b)) || (diff <= rel_tol * abs(a)) || (diff <= abs_tol)
+    Ok((diff <= rel_tol * abs(b)) || (diff <= rel_tol * abs(a)) || (diff <= abs_tol))
 }
 
 #[cfg(test)]
@@ -315,39 +331,64 @@ mod tests {
     #[test]
     fn test_isclose_basic() {
         // Equal values
-        assert!(isclose(
-            Complex64::new(1.0, 2.0),
-            Complex64::new(1.0, 2.0),
-            1e-9,
-            0.0
-        ));
+        assert_eq!(
+            isclose(
+                Complex64::new(1.0, 2.0),
+                Complex64::new(1.0, 2.0),
+                Some(1e-9),
+                Some(0.0)
+            ),
+            Ok(true)
+        );
         // Close values
-        assert!(isclose(
-            Complex64::new(1.0, 2.0),
-            Complex64::new(1.0 + 1e-10, 2.0),
-            1e-9,
-            0.0
-        ));
+        assert_eq!(
+            isclose(
+                Complex64::new(1.0, 2.0),
+                Complex64::new(1.0 + 1e-10, 2.0),
+                Some(1e-9),
+                Some(0.0)
+            ),
+            Ok(true)
+        );
         // Not close
-        assert!(!isclose(
-            Complex64::new(1.0, 2.0),
-            Complex64::new(2.0, 2.0),
-            1e-9,
-            0.0
-        ));
+        assert_eq!(
+            isclose(
+                Complex64::new(1.0, 2.0),
+                Complex64::new(2.0, 2.0),
+                Some(1e-9),
+                Some(0.0)
+            ),
+            Ok(false)
+        );
         // Infinities
-        assert!(isclose(
-            Complex64::new(f64::INFINITY, 0.0),
-            Complex64::new(f64::INFINITY, 0.0),
-            1e-9,
-            0.0
-        ));
-        assert!(!isclose(
-            Complex64::new(f64::INFINITY, 0.0),
-            Complex64::new(f64::NEG_INFINITY, 0.0),
-            1e-9,
-            0.0
-        ));
+        assert_eq!(
+            isclose(
+                Complex64::new(f64::INFINITY, 0.0),
+                Complex64::new(f64::INFINITY, 0.0),
+                Some(1e-9),
+                Some(0.0)
+            ),
+            Ok(true)
+        );
+        assert_eq!(
+            isclose(
+                Complex64::new(f64::INFINITY, 0.0),
+                Complex64::new(f64::NEG_INFINITY, 0.0),
+                Some(1e-9),
+                Some(0.0)
+            ),
+            Ok(false)
+        );
+        // Negative tolerance
+        assert!(
+            isclose(
+                Complex64::new(1.0, 2.0),
+                Complex64::new(1.0, 2.0),
+                Some(-1.0),
+                Some(0.0)
+            )
+            .is_err()
+        );
     }
 
     proptest::proptest! {

src/math.rs

@@ -96,17 +96,33 @@ pub(crate) fn math_1a(x: f64, func: fn(f64) -> f64) -> crate::Result<f64> {
     crate::err::is_error(r)
 }
 
-/// Return the Euclidean distance, sqrt(x*x + y*y).
+/// Return the Euclidean norm of n-dimensional coordinates.
 ///
-/// Uses high-precision vector_norm algorithm instead of libm hypot()
-/// for consistent results across platforms and better handling of overflow/underflow.
+/// Equivalent to sqrt(sum(x**2 for x in coords)).
+/// Uses high-precision vector_norm algorithm for consistent results
+/// across platforms and better handling of overflow/underflow.
 #[inline]
-pub fn hypot(x: f64, y: f64) -> f64 {
-    let ax = x.abs();
-    let ay = y.abs();
-    let max = if ax > ay { ax } else { ay };
-    let found_nan = x.is_nan() || y.is_nan();
-    aggregate::vector_norm_2(ax, ay, max, found_nan)
+pub fn hypot(coords: &[f64]) -> f64 {
+    let n = coords.len();
+    if n == 0 {
+        return 0.0;
+    }
+
+    let mut max = 0.0_f64;
+    let mut found_nan = false;
+    let abs_coords: Vec<f64> = coords
+        .iter()
+        .map(|&x| {
+            let ax = x.abs();
+            found_nan |= ax.is_nan();
+            if ax > max {
+                max = ax;
+            }
+            ax
+        })
+        .collect();
+
+    aggregate::vector_norm(&abs_coords, max, found_nan)
 }
 
 // Mathematical constants
@@ -208,16 +224,43 @@ mod tests {
         assert!(NAN.is_nan());
     }
 
-    // hypot tests
-    fn test_hypot(x: f64, y: f64) {
-        crate::test::test_math_2(x, y, "hypot", |x, y| Ok(hypot(x, y)));
+    // hypot tests - n-dimensional
+    fn test_hypot(coords: &[f64]) {
+        let rs_result = hypot(coords);
+
+        pyo3::Python::attach(|py| {
+            let math = PyModule::import(py, "math").unwrap();
+            let py_func = math.getattr("hypot").unwrap();
+            let py_args = pyo3::types::PyTuple::new(py, coords).unwrap();
+            let py_result: f64 = py_func.call1(py_args).unwrap().extract().unwrap();
+
+            if py_result.is_nan() && rs_result.is_nan() {
+                return;
+            }
+            assert_eq!(
+                py_result.to_bits(),
+                rs_result.to_bits(),
+                "hypot({:?}): py={} vs rs={}",
+                coords,
+                py_result,
+                rs_result
+            );
+        });
+    }
+
+    #[test]
+    fn test_hypot_basic() {
+        test_hypot(&[3.0, 4.0]); // 5.0
+        test_hypot(&[1.0, 2.0, 2.0]); // 3.0
+        test_hypot(&[]); // 0.0
+        test_hypot(&[5.0]); // 5.0
     }
 
     #[test]
     fn edgetest_hypot() {
         for &x in &crate::test::EDGE_VALUES {
             for &y in &crate::test::EDGE_VALUES {
-                test_hypot(x, y);
+                test_hypot(&[x, y]);
             }
         }
     }
@@ -235,7 +278,7 @@ mod tests {
 
         #[test]
         fn proptest_hypot(x: f64, y: f64) {
-            test_hypot(x, y);
+            test_hypot(&[x, y]);
         }
     }
 }

src/math/aggregate.rs

@@ -170,80 +170,8 @@ pub fn fsum(iter: impl IntoIterator<Item = f64>) -> crate::Result<f64> {
 
 // VECTOR_NORM - for dist and hypot
 
-/// Compute the Euclidean norm of two values with high precision.
-/// Optimized version for hypot(x, y).
-pub(super) fn vector_norm_2(x: f64, y: f64, max: f64, found_nan: bool) -> f64 {
-    // Check for infinity first (inf wins over nan)
-    if x.is_infinite() || y.is_infinite() {
-        return f64::INFINITY;
-    }
-    if found_nan {
-        return f64::NAN;
-    }
-    if max == 0.0 {
-        return 0.0;
-    }
-    // n == 1 case: only one non-zero value
-    if x == 0.0 || y == 0.0 {
-        return max;
-    }
-
-    let mut max_e: i32 = 0;
-    crate::m::frexp(max, &mut max_e);
-
-    if max_e < -1023 {
-        // When max_e < -1023, ldexp(1.0, -max_e) would overflow
-        return f64::MIN_POSITIVE
-            * vector_norm_2(
-                x / f64::MIN_POSITIVE,
-                y / f64::MIN_POSITIVE,
-                max / f64::MIN_POSITIVE,
-                found_nan,
-            );
-    }
-
-    let scale = crate::m::ldexp(1.0, -max_e);
-    debug_assert!(max * scale >= 0.5);
-    debug_assert!(max * scale < 1.0);
-
-    let mut csum = 1.0;
-    let mut frac1 = 0.0;
-    let mut frac2 = 0.0;
-
-    // Process x
-    let xs = x * scale;
-    debug_assert!(xs.abs() < 1.0);
-    let pr = dl_mul(xs, xs);
-    debug_assert!(pr.hi <= 1.0);
-    let sm = dl_fast_sum(csum, pr.hi);
-    csum = sm.hi;
-    frac1 += pr.lo;
-    frac2 += sm.lo;
-
-    // Process y
-    let ys = y * scale;
-    debug_assert!(ys.abs() < 1.0);
-    let pr = dl_mul(ys, ys);
-    debug_assert!(pr.hi <= 1.0);
-    let sm = dl_fast_sum(csum, pr.hi);
-    csum = sm.hi;
-    frac1 += pr.lo;
-    frac2 += sm.lo;
-
-    let mut h = (csum - 1.0 + (frac1 + frac2)).sqrt();
-    let pr = dl_mul(-h, h);
-    let sm = dl_fast_sum(csum, pr.hi);
-    csum = sm.hi;
-    frac1 += pr.lo;
-    frac2 += sm.lo;
-    let x = csum - 1.0 + (frac1 + frac2);
-    h += x / (2.0 * h); // differential correction
-
-    h / scale
-}
-
 /// Compute the Euclidean norm of a vector with high precision.
-fn vector_norm(vec: &[f64], max: f64, found_nan: bool) -> f64 {
+pub fn vector_norm(vec: &[f64], max: f64, found_nan: bool) -> f64 {
     let n = vec.len();
 
     if max.is_infinite() {