Merge pull request #128 from acgetchell/fix/94-tolerance-invariants

fix(matrix): reject overflowed symmetry tolerance scaling

Author: acgetchell

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Feat!(matrix): enforce fallible matrix invariants [`e26c283`](https://github.com/acgetchell/la-stack/commit/e26c28358b2358100353b2895441b68892e92cd7)
+- Feat!(api): enforce fallible numeric invariants [`adfc33b`](https://github.com/acgetchell/la-stack/commit/adfc33b945b259721bd1067e797ed2e7d4ec0e6e)
 
 ### Changed
 

examples/exact_det_3x3.rs

@@ -9,7 +9,7 @@
 
 use la_stack::prelude::*;
 
-fn main() {
+fn main() -> Result<(), LaError> {
     // Base matrix: rows in arithmetic progression → exactly singular (det = 0).
     //   [[1, 2, 3],
     //    [4, 5, 6],
@@ -24,9 +24,11 @@ fn main() {
         [7.0, 8.0, 9.0],
     ]);
 
-    let det_f64_approx = m.det_direct().unwrap().unwrap();
-    let det_exact = m.det_exact().unwrap();
-    let det_exact_as_f64 = m.det_exact_f64().unwrap();
+    let Some(det_f64_approx) = m.det_direct()? else {
+        unreachable!("D=3 is supported by det_direct");
+    };
+    let det_exact = m.det_exact()?;
+    let det_exact_as_f64 = m.det_exact_f64()?;
 
     println!("Near-singular 3×3 matrix (perturbation = 2^-50 ≈ {perturbation:.2e}):");
     for r in 0..3 {
@@ -35,7 +37,7 @@ fn main() {
             if c > 0 {
                 print!(", ");
             }
-            print!("{:22.18}", m.get(r, c).unwrap());
+            print!("{:22.18}", m.get_checked(r, c)?);
         }
         println!("]");
     }
@@ -45,4 +47,5 @@ fn main() {
     println!("det_exact_f64()    = {det_exact_as_f64:+.6e}");
     println!();
     println!("The exact determinant is −3/2^50 ≈ −2.66e-15.");
+    Ok(())
 }

examples/exact_solve_3x3.rs

@@ -9,7 +9,7 @@
 
 use la_stack::prelude::*;
 
-fn main() {
+fn main() -> Result<(), LaError> {
     // Near-singular 3×3 system.
     //
     // The base matrix [[1,2,3],[4,5,6],[7,8,9]] is exactly singular (rows in
@@ -26,16 +26,11 @@ fn main() {
 
     // f64 LU solve (using zero pivot tolerance since the matrix is nearly singular
     // and would be rejected by DEFAULT_PIVOT_TOL).
-    let lu_x = a
-        .lu(Tolerance::new(0.0).unwrap())
-        .unwrap()
-        .solve_vec(b)
-        .unwrap()
-        .into_array();
+    let lu_x = a.lu(Tolerance::new(0.0)?)?.solve_vec(b)?.into_array();
 
     // Exact solve.
-    let exact_x = a.solve_exact(b).unwrap();
-    let exact_x_f64 = a.solve_exact_f64(b).unwrap().into_array();
+    let exact_x = a.solve_exact(b)?;
+    let exact_x_f64 = a.solve_exact_f64(b)?.into_array();
 
     println!("Near-singular 3×3 system (perturbation = 2^-50 ≈ {perturbation:.2e}):");
     for r in 0..3 {
@@ -44,7 +39,7 @@ fn main() {
             if c > 0 {
                 print!(", ");
             }
-            print!("{:22.18}", a.get(r, c).unwrap());
+            print!("{:22.18}", a.get_checked(r, c)?);
         }
         println!("]");
     }
@@ -67,4 +62,5 @@ fn main() {
         "solve_exact_f64():  x = [{:+.6e}, {:+.6e}, {:+.6e}]",
         exact_x_f64[0], exact_x_f64[1], exact_x_f64[2]
     );
+    Ok(())
 }

src/lib.rs

@@ -197,6 +197,9 @@ pub const MAX_STACK_MATRIX_DISPATCH_DIM: usize = 7;
 /// Construct with [`Tolerance::new`] when accepting raw caller input. Once
 /// constructed, the stored value is guaranteed to be finite and `>= 0`, so
 /// downstream algorithms do not need to revalidate the tolerance.
+///
+/// This is the crate-wide tolerance contract: raw negative, NaN, and infinite
+/// values are rejected with [`LaError::InvalidTolerance`] at construction time.
 #[must_use]
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct Tolerance {
@@ -488,6 +491,11 @@ impl LaError {
     /// use la_stack::prelude::*;
     ///
     /// assert_eq!(LaError::validate_tolerance(1e-12)?.get(), 1e-12);
+    ///
+    /// let raw = 0.0;
+    /// let tol = LaError::validate_tolerance(raw)?;
+    /// let _lu = Matrix::<2>::identity().lu(tol)?;
+    ///
     /// assert_eq!(
     ///     LaError::validate_tolerance(-1.0),
     ///     Err(LaError::InvalidTolerance { value: -1.0 })
@@ -667,6 +675,8 @@ pub mod prelude {
 
 #[cfg(test)]
 mod tests {
+    use core::assert_matches;
+
     use super::*;
 
     use approx::assert_abs_diff_eq;
@@ -681,6 +691,74 @@ mod tests {
         );
     }
 
+    #[test]
+    fn tolerance_new_accepts_finite_non_negative_values() {
+        assert_eq!(
+            Tolerance::new(0.0).unwrap().get().to_bits(),
+            0.0f64.to_bits()
+        );
+        assert_eq!(
+            Tolerance::new(1e-12).unwrap().get().to_bits(),
+            1e-12f64.to_bits()
+        );
+        assert_eq!(
+            Tolerance::new(f64::MAX).unwrap().get().to_bits(),
+            f64::MAX.to_bits()
+        );
+    }
+
+    #[test]
+    fn tolerance_new_rejects_negative_nan_and_infinity() {
+        assert_eq!(
+            Tolerance::new(-1.0),
+            Err(LaError::InvalidTolerance { value: -1.0 })
+        );
+        assert_matches!(
+            Tolerance::new(f64::NAN),
+            Err(LaError::InvalidTolerance { value }) if value.is_nan()
+        );
+        assert_eq!(
+            Tolerance::new(f64::INFINITY),
+            Err(LaError::InvalidTolerance {
+                value: f64::INFINITY,
+            })
+        );
+        assert_eq!(
+            Tolerance::new(f64::NEG_INFINITY),
+            Err(LaError::InvalidTolerance {
+                value: f64::NEG_INFINITY,
+            })
+        );
+    }
+
+    #[test]
+    fn validate_tolerance_matches_tolerance_new() {
+        for value in [0.0, 1e-12, f64::MAX] {
+            assert_eq!(LaError::validate_tolerance(value), Tolerance::new(value));
+        }
+
+        assert_eq!(
+            LaError::validate_tolerance(-1.0),
+            Err(LaError::InvalidTolerance { value: -1.0 })
+        );
+        assert_matches!(
+            LaError::validate_tolerance(f64::NAN),
+            Err(LaError::InvalidTolerance { value }) if value.is_nan()
+        );
+        assert_eq!(
+            LaError::validate_tolerance(f64::INFINITY),
+            Err(LaError::InvalidTolerance {
+                value: f64::INFINITY,
+            })
+        );
+        assert_eq!(
+            LaError::validate_tolerance(f64::NEG_INFINITY),
+            Err(LaError::InvalidTolerance {
+                value: f64::NEG_INFINITY,
+            })
+        );
+    }
+
     #[test]
     fn laerror_display_formats_singular() {
         let err = LaError::Singular { pivot_col: 3 };

src/matrix.rs

@@ -195,6 +195,10 @@ impl<const D: usize> Matrix<D> {
     /// Non-finite entries are rejected with source coordinates instead of
     /// silently propagating NaN or infinity through the norm.
     ///
+    /// Row sums are accumulated in `f64` with ordinary addition.  This method
+    /// checks for non-finite inputs and overflowed accumulators, but it does not
+    /// provide a certified absolute rounding bound for the returned norm.
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -263,9 +267,17 @@ impl<const D: usize> Matrix<D> {
     /// Use [`first_asymmetry`](Self::first_asymmetry) to locate the first
     /// offending pair when this returns `Ok(false)`.
     ///
+    /// The `rel_tol` argument is a [`Tolerance`], so raw caller input must be
+    /// finite and non-negative before it can reach this predicate. Use
+    /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
+    /// raw `f64`; negative, NaN, and infinite tolerances return
+    /// [`LaError::InvalidTolerance`].
+    ///
     /// # NaN / infinity handling
     /// Stored NaN or ±∞ entries return [`LaError::NonFinite`] with the
-    /// offending matrix coordinates.  If both stored entries are finite but
+    /// offending matrix coordinates.  A finite matrix can still return
+    /// [`LaError::NonFinite`] if computing the scaled symmetry tolerance
+    /// overflows to NaN or infinity.  If both stored entries are finite but
     /// their difference overflows to ±∞, the pair is reported as asymmetric.
     ///
     /// # Examples
@@ -284,7 +296,9 @@ impl<const D: usize> Matrix<D> {
     /// ```
     ///
     /// # Errors
-    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite.
+    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite,
+    /// or when computing the scaled symmetry tolerance overflows to NaN or
+    /// infinity.
     #[inline]
     pub fn is_symmetric(&self, rel_tol: Tolerance) -> Result<bool, LaError> {
         Ok(self.first_asymmetry(rel_tol)?.is_none())
@@ -299,6 +313,18 @@ impl<const D: usize> Matrix<D> {
     /// predicate is the same as [`is_symmetric`](Self::is_symmetric):
     /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`.
     ///
+    /// Stored NaN or ±∞ entries return [`LaError::NonFinite`] with the
+    /// offending matrix coordinates. A finite matrix can still return
+    /// [`LaError::NonFinite`] if computing the scaled symmetry tolerance
+    /// overflows to NaN or infinity. If both stored entries are finite but
+    /// their difference overflows to ±∞, the pair is reported as asymmetric.
+    ///
+    /// The `rel_tol` argument is a [`Tolerance`], so raw caller input must be
+    /// finite and non-negative before it can reach this predicate. Use
+    /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
+    /// raw `f64`; negative, NaN, and infinite tolerances return
+    /// [`LaError::InvalidTolerance`].
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -317,7 +343,9 @@ impl<const D: usize> Matrix<D> {
     /// ```
     ///
     /// # Errors
-    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite.
+    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite,
+    /// or when computing the scaled symmetry tolerance overflows to NaN or
+    /// infinity.
     #[inline]
     pub fn first_asymmetry(&self, rel_tol: Tolerance) -> Result<Option<(usize, usize)>, LaError> {
         let eps = self.symmetry_epsilon(rel_tol)?;
@@ -351,7 +379,9 @@ impl<const D: usize> Matrix<D> {
     /// off-diagonal mismatches.
     ///
     /// # Errors
-    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite.
+    /// Returns [`LaError::NonFinite`] when any matrix entry is NaN or infinite,
+    /// or when computing the scaled symmetry tolerance overflows to NaN or
+    /// infinity.
     fn symmetry_epsilon(&self, rel_tol: Tolerance) -> Result<f64, LaError> {
         let rel_tol = rel_tol.get();
         let mut eps = rel_tol;
@@ -365,6 +395,10 @@ impl<const D: usize> Matrix<D> {
                     return Err(LaError::non_finite_cell(r, c));
                 }
                 row_eps = rel_tol.mul_add(entry.abs(), row_eps);
+                if !row_eps.is_finite() {
+                    cold_path();
+                    return Err(LaError::non_finite_at(c));
+                }
             }
             if row_eps > eps {
                 eps = row_eps;
@@ -393,6 +427,12 @@ impl<const D: usize> Matrix<D> {
     /// # }
     /// ```
     ///
+    /// The `tol` argument is a [`Tolerance`], so raw caller input must be
+    /// finite and non-negative before it can reach factorization. Use
+    /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
+    /// raw `f64`; negative, NaN, and infinite tolerances return
+    /// [`LaError::InvalidTolerance`].
+    ///
     /// # Errors
     /// Returns [`LaError::Singular`] if, for some column `k`, the largest-magnitude candidate pivot
     /// in that column satisfies `|pivot| <= tol` (so no numerically usable pivot exists).
@@ -416,6 +456,12 @@ impl<const D: usize> Matrix<D> {
     /// general-purpose factorization that tolerates non-symmetric inputs, use
     /// [`lu`](Self::lu) instead.
     ///
+    /// The `tol` argument is a [`Tolerance`], so raw caller input must be
+    /// finite and non-negative before it can reach factorization. Use
+    /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
+    /// raw `f64`; negative, NaN, and infinite tolerances return
+    /// [`LaError::InvalidTolerance`].
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -569,6 +615,12 @@ impl<const D: usize> Matrix<D> {
     /// speedup (see [`det_direct`](Self::det_direct)). The `tol` parameter is only used
     /// by the LU fallback path for D ≥ 5.
     ///
+    /// The `tol` argument is a [`Tolerance`], so raw caller input must be
+    /// finite and non-negative before it can reach the determinant path. Use
+    /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
+    /// raw `f64`; negative, NaN, and infinite tolerances return
+    /// [`LaError::InvalidTolerance`].
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -1464,6 +1516,21 @@ mod tests {
         assert!(!a.is_symmetric(Tolerance::new(0.0).unwrap()).unwrap());
     }
 
+    #[test]
+    fn first_asymmetry_rejects_scaled_epsilon_overflow() {
+        let a = Matrix::<2>::from_rows([[2.0, 0.0], [0.0, 1.0]]);
+        let tol = Tolerance::new(f64::MAX).unwrap();
+
+        assert_eq!(
+            a.first_asymmetry(tol),
+            Err(LaError::NonFinite { row: None, col: 0 })
+        );
+        assert_eq!(
+            a.is_symmetric(tol),
+            Err(LaError::NonFinite { row: None, col: 0 })
+        );
+    }
+
     #[test]
     fn first_asymmetry_flags_overflowed_finite_difference() {
         let a = Matrix::<2>::from_rows([[1.0, f64::MAX], [-f64::MAX, 1.0]]);