Skip to content

fix(matrix): reject overflowed symmetry tolerance scaling #128

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
acgetchell merged 1 commit into from
Jun 3, 2026
Merged

fix(matrix): reject overflowed symmetry tolerance scaling #128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot Jun 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major | ⚡ Quick win

Do not edit CHANGELOG.md manually.

This entry should be generated via the changelog tooling, not committed as a direct file edit. Please regenerate using just changelog (or just changelog-unreleased <version> for unreleased prepend) and commit the generated result instead.

Based on learnings: "Never edit CHANGELOG.md directly - it's auto-generated from git commits. Use just changelog to regenerate or just changelog-unreleased <version> to prepend unreleased changes."

🤖 Prompt for AI Agents 
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@CHANGELOG.md` at line 13, The CHANGELOG.md contains a manual entry
("Feat!(api): enforce fallible numeric invariants [`adfc33b`]") that must not be
edited by hand; revert the manual edit to CHANGELOG.md and regenerate the file
using the changelog tooling (run either `just changelog` or `just
changelog-unreleased <version>` to prepend unreleased changes), then commit the
regenerated CHANGELOG.md instead of the manual change so the entry for commit
adfc33b is produced by the tool.


### Changed

Expand Down
13 changes: 8 additions & 5 deletions examples/exact_det_3x3.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand All @@ -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 {
Expand All @@ -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!("]");
}
Expand All @@ -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(())
}
16 changes: 6 additions & 10 deletions examples/exact_solve_3x3.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 {
Expand All @@ -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!("]");
}
Expand All @@ -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(())
}
78 changes: 78 additions & 0 deletions src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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 })
Expand Down Expand Up @@ -667,6 +675,8 @@ pub mod prelude {

#[cfg(test)]
mod tests {
use core::assert_matches;

use super::*;

use approx::assert_abs_diff_eq;
Expand All @@ -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 };
Expand Down
75 changes: 71 additions & 4 deletions src/matrix.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::*;
Expand Down Expand Up @@ -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
Expand All @@ -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())
Expand All @@ -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::*;
Expand All @@ -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)?;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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).
Expand All @@ -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::*;
Expand Down Expand Up @@ -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::*;
Expand Down Expand Up @@ -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]]);
Expand Down
Loading
Loading