diff --git a/CHANGELOG.md b/CHANGELOG.md index 31b845a..798be0d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/examples/exact_det_3x3.rs b/examples/exact_det_3x3.rs index 1bfc832..6ecadc3 100644 --- a/examples/exact_det_3x3.rs +++ b/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(()) } diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs index 88cc1ee..671226d 100644 --- a/examples/exact_solve_3x3.rs +++ b/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(()) } diff --git a/src/lib.rs b/src/lib.rs index cca1fbc..f9f3513 100644 --- a/src/lib.rs +++ b/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 }; diff --git a/src/matrix.rs b/src/matrix.rs index 2dae220..82f3597 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -195,6 +195,10 @@ impl Matrix { /// 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 Matrix { /// 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 Matrix { /// ``` /// /// # 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 { Ok(self.first_asymmetry(rel_tol)?.is_none()) @@ -299,6 +313,18 @@ impl Matrix { /// 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 Matrix { /// ``` /// /// # 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, LaError> { let eps = self.symmetry_epsilon(rel_tol)?; @@ -351,7 +379,9 @@ impl Matrix { /// 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 { let rel_tol = rel_tol.get(); let mut eps = rel_tol; @@ -365,6 +395,10 @@ impl Matrix { 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 Matrix { /// # } /// ``` /// + /// 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 Matrix { /// 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 Matrix { /// 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]]); diff --git a/src/vector.rs b/src/vector.rs index 3f14c34..1b7c2ad 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -71,6 +71,12 @@ impl Vector { /// Dot product. /// + /// Terms are accumulated in `f64` using [`f64::mul_add`] at each index. + /// Intermediate rounding occurs, and this method does not provide a + /// certified absolute rounding bound for the returned dot product. The + /// returned [`Result`] is still checked for non-finite inputs and for + /// non-finite accumulation. + /// /// # Examples /// ``` /// use la_stack::prelude::*; @@ -108,6 +114,13 @@ impl Vector { /// Squared Euclidean norm. /// + /// This is computed as `dot(self, self)`, so `norm2_sq` has the same + /// `f64` [`mul_add`](f64::mul_add) accumulation behavior as [`dot`](Self::dot). + /// Intermediate rounding occurs, and this method does not provide a + /// certified absolute rounding bound for the returned squared norm. The + /// returned [`Result`] is still checked for non-finite inputs and for + /// non-finite accumulation. + /// /// # Examples /// ``` /// use la_stack::prelude::*; @@ -271,6 +284,30 @@ mod tests { }) ); } + + #[test] + fn []() { + let a = Vector::<$d>::new([1.0; $d]); + let mut b_arr = [1.0f64; $d]; + b_arr[0] = f64::INFINITY; + let b = Vector::<$d>::new(b_arr); + + assert_eq!(a.dot(b), Err(LaError::NonFinite { row: None, col: 0 })); + } + + #[test] + fn []() { + let mut a_arr = [1.0f64; $d]; + a_arr[0] = f64::MAX; + let a = Vector::<$d>::new(a_arr); + + let mut b_arr = [1.0f64; $d]; + b_arr[0] = 2.0; + let b = Vector::<$d>::new(b_arr); + + assert_eq!(a.dot(b), Err(LaError::NonFinite { row: None, col: 0 })); + assert_eq!(a.norm2_sq(), Err(LaError::NonFinite { row: None, col: 0 })); + } } }; }