From a62e72a0468ea13c7c5e851dbf195123d5bb0d38 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Fri, 4 Sep 2026 17:45:39 -0700 Subject: [PATCH] feat: add const-generic Gram matrix construction - Add allocation-free, const-evaluable `gram_matrix` with independent vector count and dimension - Preserve bitwise symmetry and typed dot-product overflow diagnostics - Document geometric uses, conditioning, and floating-point limitations - Add benchmarks for square and rectangular vector collections - Simplify exact rational scaling using canonical positive denominators - Deduplicate factorization property-test fixtures and assertions --- Cargo.toml | 5 + README.md | 9 ++ REFERENCES.md | 20 +++ benches/gram.rs | 113 ++++++++++++++++ src/gram.rs | 73 ++++++++++ src/lib.rs | 8 +- src/rational.rs | 28 +--- tests/proptest_factorizations.rs | 224 +++++++++---------------------- tests/proptest_gram.rs | 213 +++++++++++++++++++++++++++++ 9 files changed, 513 insertions(+), 180 deletions(-) create mode 100644 benches/gram.rs create mode 100644 src/gram.rs create mode 100644 tests/proptest_gram.rs diff --git a/Cargo.toml b/Cargo.toml index a07212f..25a6a0d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -91,6 +91,11 @@ name = "linear_form" harness = false required-features = [ "bench" ] +[[bench]] +name = "gram" +harness = false +required-features = [ "bench" ] + [profile.release] lto = "fat" codegen-units = 1 diff --git a/README.md b/README.md index 0d65a94..0d1338d 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,15 @@ while keeping the API intentionally small and explicit. `la-stack` provides a handful of const-generic, stack-backed building blocks: - `Vector` for fixed-length `f64` vectors backed by `[f64; D]` +- `gram_matrix(&[Vector; M])` for allocation-free `Matrix` construction + from pairwise vector inner products, with bit-for-bit symmetry. Gram matrices + encode lengths and angles and support simplex/facet volume calculations; see + [Gram matrices and geometric measures](REFERENCES.md#gram-matrices-and-geometric-measures). + Each independent dot product is checked once; + rounding has no certified error bound, and positive definiteness or affine + independence must still be established by factorization or the caller. + Benchmark square simplex and rectangular facet inputs through dimension 8 + with `cargo bench --locked --features bench --bench gram`. - `Matrix` for fixed-size square `f64` matrices backed by `[[f64; D]; D]` - `Interval` and `IntervalMatrix` for outward-rounded, proof-bearing determinant filters through D=7 diff --git a/REFERENCES.md b/REFERENCES.md index 49b0237..2368ad5 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -23,6 +23,23 @@ No generated content was used without human oversight. ## Linear algebra algorithms +### Gram matrices and geometric measures + +A Gram matrix collects pairwise inner products: `G[i,j] = v_i · v_j`. +Writing the vectors as rows of `V` gives `G = V Vᵀ`. Its diagonal contains +squared lengths; off-diagonal entries describe angles through +`v_i · v_j = ||v_i|| ||v_j|| cos(θ)` for nonzero vectors. + +In exact real arithmetic, `G` is positive semidefinite and is positive definite +exactly when the vectors are linearly independent. For `M ≤ N`, `det(G)` is +the squared M-dimensional volume spanned by the vectors. For simplex edge +vectors from a common vertex, the simplex volume is `sqrt(det(G)) / M!` [16]. +This applies to triangles embedded in 3D and to higher-dimensional facets. + +`gram_matrix` computes rounded binary64 entries with exact mirrored symmetry; +it does not certify rank, positive definiteness, or volume accuracy. See [9-12] +for floating-point and conditioning background. + ### Certified fixed-vector reductions `Vector::dot_with_errbound()` and `Vector::dot_difference_with_errbound()` use @@ -202,3 +219,6 @@ finite results from overflow. 15. Blue, James L. "A Portable Fortran Program to Find the Euclidean Norm of a Vector." *ACM Transactions on Mathematical Software* 4.1 (1978): 15–23. [DOI](https://doi.org/10.1145/355769.355771) +16. Kock, Anders. "Square-densities, and volume forms." Notes, December 10, 2020. + Introduction and §1.2 (Gram's formula). + [Author's PDF](https://math.au.dk/~kock/heron4.pdf) diff --git a/benches/gram.rs b/benches/gram.rs new file mode 100644 index 0000000..b9e27c9 --- /dev/null +++ b/benches/gram.rs @@ -0,0 +1,113 @@ +#![forbid(unsafe_code)] + +//! Gram construction versus checked hand-written assembly, excluding setup. + +use core::array::from_fn; +use std::hint::black_box; + +use criterion::Criterion; + +use la_stack::{LaError, Matrix, Vector, gram_matrix}; + +#[path = "common/bench_utils.rs"] +mod bench_utils; +use bench_utils::OrAbort; + +/// Fixture families whose labels and construction must agree. +enum Scenario { + Orthogonal, + Dependent, + NearDependent, + MixedScale, +} + +impl Scenario { + /// Stable label used in Criterion result paths. + const fn name(&self) -> &'static str { + match self { + Self::Orthogonal => "orthogonal", + Self::Dependent => "dependent", + Self::NearDependent => "near_dependent", + Self::MixedScale => "mixed_scale", + } + } + + /// Small integer entries keep the independent matrix-product oracle exact. + fn entry(&self, row: usize, coordinate: usize) -> i16 { + let diagonal = i16::from(row == coordinate); + match self { + Self::Orthogonal => diagonal, + Self::Dependent => 1, + Self::NearDependent => 256 + diagonal, + Self::MixedScale => 257 * diagonal - 1, + } + } +} + +fn hand_written( + vectors: &[Vector; M], +) -> Result, LaError> { + let mut matrix = Matrix::zero(); + for (i, left) in vectors.iter().enumerate() { + for (j, right) in vectors.iter().enumerate().skip(i) { + let value = left.dot(right)?; + matrix.set(i, j, value)?; + matrix.set(j, i, value)?; + } + } + Ok(matrix) +} + +fn register(c: &mut Criterion) { + for scenario in [ + Scenario::Orthogonal, + Scenario::Dependent, + Scenario::NearDependent, + Scenario::MixedScale, + ] { + let integers: [[i16; N]; M] = from_fn(|i| from_fn(|k| scenario.entry(i, k))); + let vectors = integers + .map(|row| Vector::try_new(row.map(f64::from)).or_abort("Gram benchmark input")); + let expected = Matrix::try_from_rows(from_fn(|i| { + from_fn(|j| { + let value: i32 = (0..N) + .map(|k| i32::from(integers[i][k]) * i32::from(integers[j][k])) + .sum(); + f64::from(value) + }) + })) + .or_abort("Gram oracle"); + assert_eq!(gram_matrix(&vectors).or_abort("Gram validation"), expected); + assert_eq!( + hand_written(&vectors).or_abort("hand-written validation"), + expected + ); + let mut group = c.benchmark_group(format!("gram/{M}x{N}/{}", scenario.name())); + group.bench_function("la_stack", |b| { + b.iter(|| black_box(gram_matrix(black_box(&vectors)).or_abort("Gram construction"))); + }); + group.bench_function("hand_written", |b| { + b.iter(|| black_box(hand_written(black_box(&vectors)).or_abort("Gram assembly"))); + }); + group.finish(); + } +} + +fn main() { + let mut c = Criterion::default().configure_from_args(); + register::<2, 2>(&mut c); + register::<1, 2>(&mut c); + register::<3, 3>(&mut c); + register::<2, 3>(&mut c); + register::<4, 4>(&mut c); + register::<3, 4>(&mut c); + register::<5, 5>(&mut c); + register::<4, 5>(&mut c); + register::<6, 6>(&mut c); + register::<5, 6>(&mut c); + register::<7, 7>(&mut c); + register::<6, 7>(&mut c); + register::<8, 8>(&mut c); + register::<7, 8>(&mut c); + c.final_summary(); +} diff --git a/src/gram.rs b/src/gram.rs new file mode 100644 index 0000000..f10bd4d --- /dev/null +++ b/src/gram.rs @@ -0,0 +1,73 @@ +#![forbid(unsafe_code)] + +//! Fixed-size Gram construction. + +use crate::{LaError, Matrix, Vector}; + +/// Construct a stack-backed [`Matrix`] of pairwise vector dot products. +/// +/// A Gram matrix records pairwise inner products: diagonal entries are squared +/// vector lengths, and off-diagonal entries encode their relative angles. +/// The input contains `M` finite-by-construction [`Vector`] values. +/// If `V` has these vectors as rows, the mathematical Gram matrix is `G = V Vᵀ`, +/// with `G[i,j] = vectors[i] · vectors[j]`. In exact real arithmetic and for +/// `M ≤ N`, its determinant is the squared volume of the spanned +/// parallelotope. For edges from one simplex vertex, the simplex volume is +/// `sqrt(det(G)) / M!`; this also handles facets embedded in higher dimensions. +/// See `REFERENCES.md` \[16\] for the Gram determinant and volume interpretation. +/// +/// Each upper-triangle dot product is computed once using [`Vector::dot`]'s +/// left-to-right fused multiply-add reduction and copied to the other triangle, +/// giving bit-for-bit symmetry. No absolute rounding-error bound is provided. +/// Rounding and underflow can destroy positive semidefiniteness or rank; this +/// operation proves neither positive definiteness nor affine independence. +/// [`Matrix::ldlt`] retains its symmetry and positive-definiteness preconditions. +/// Forming a Gram matrix squares the spectral condition number of an exact input +/// with full row rank. See the floating-point discussion in `REFERENCES.md` \[9-11\]. +/// +/// `M` and `N` are independent, with no dimension cap (including dimensions +/// through 8); storage is `O(M²)` and work is `O(M²(N + 1))`, including output +/// initialization when `N = 0`. `M = 0` returns an empty +/// matrix; `N = 0` returns an all-zero matrix. No optional feature is required. +/// +/// # Errors +/// Returns [`LaError::NonFinite`] if a dot-product accumulator overflows, even +/// when the exact result would be finite after cancellation. The error preserves +/// [`Vector::dot`]'s [`Computation`](crate::NonFiniteOrigin::Computation) origin, +/// [`VectorDotProduct`](crate::ArithmeticOperation::VectorDotProduct) operation, +/// and first failing reduction [`Step`](crate::NonFiniteLocation::Step). +/// The step indexes the vector coordinate, not the output matrix cell. +/// Pairs are visited in upper-triangle row order. +/// +/// # Examples +/// ``` +/// use la_stack::prelude::*; +/// # fn main() -> Result<(), LaError> { +/// let vectors = [Vector::try_new([1.0, 0.0, 0.0])?, +/// Vector::try_new([0.0, 2.0, 0.0])?]; +/// let gram = gram_matrix(&vectors)?; +/// assert_eq!(gram.as_rows(), &[[1.0, 0.0], [0.0, 4.0]]); +/// assert_eq!(gram.ldlt(Tolerance::try_new(0.0)?)?.det()?, 4.0); +/// # Ok(()) +/// # } +/// ``` +pub const fn gram_matrix( + vectors: &[Vector; M], +) -> Result, LaError> { + let mut rows = [[0.0; M]; M]; + let mut i = 0; + while i < M { + let mut j = i; + while j < M { + let value = match vectors[i].dot(&vectors[j]) { + Ok(value) => value, + Err(error) => return Err(error), + }; + rows[i][j] = value; + rows[j][i] = value; + j += 1; + } + i += 1; + } + Matrix::try_from_rows(rows) +} diff --git a/src/lib.rs b/src/lib.rs index 6744f9e..87b8cb9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -317,6 +317,7 @@ mod readme_doctests { mod error; #[cfg(feature = "exact")] mod exact; +mod gram; mod interval; mod ldlt; mod lu; @@ -507,6 +508,7 @@ pub use error::{ LaError, NonFiniteLocation, NonFiniteOrigin, PositiveSemidefiniteViolation, SingularityReason, UnrepresentableReason, }; +pub use gram::gram_matrix; pub use interval::{Interval, IntervalDeterminantSign, IntervalMatrix, MAX_INTERVAL_MATRIX_DIM}; pub use ldlt::Ldlt; pub use lu::Lu; @@ -770,7 +772,8 @@ macro_rules! try_with_rational_matrix { /// [`DeterminantWithErrorBound`], [`Interval`], [`IntervalMatrix`], /// [`IntervalDeterminantSign`], [`ScalarWithErrorBound`], [`Vector`], [`Lu`], /// [`Ldlt`], [`Tolerance`], -/// and [`LaError`]. Its typed +/// and [`LaError`]. It also includes [`gram_matrix`] for constructing a symmetric +/// matrix of pairwise vector inner products. Its typed /// error categories include [`ArithmeticOperation`], [`FactorizationKind`], /// [`IntervalBound`], [`IntervalOperand`], [`InvalidToleranceReason`], /// [`NonFiniteLocation`], [`NonFiniteOrigin`], [`PositiveSemidefiniteViolation`], @@ -834,7 +837,8 @@ pub mod prelude { InvalidToleranceReason, LaError, Ldlt, Lu, MAX_INTERVAL_MATRIX_DIM, MAX_STACK_MATRIX_DISPATCH_DIM, Matrix, NonFiniteLocation, NonFiniteOrigin, PositiveSemidefiniteViolation, ScalarWithErrorBound, SingularityReason, Tolerance, - UnrepresentableReason, Vector, try_with_interval_matrix, try_with_stack_matrix, + UnrepresentableReason, Vector, gram_matrix, try_with_interval_matrix, + try_with_stack_matrix, }; #[cfg(feature = "exact")] diff --git a/src/rational.rs b/src/rational.rs index 51f2c8b..02f6679 100644 --- a/src/rational.rs +++ b/src/rational.rs @@ -342,36 +342,22 @@ fn canonicalize_rational(value: BigRational) -> BigRational { BigRational::new(numerator, denominator) } -/// Return a positive least common multiple of all raw denominator magnitudes. +/// Return the least common multiple of canonical positive denominators. fn common_denominator<'a>(values: impl Iterator) -> BigInt { values.fold(BigInt::from(1), |scale, value| { - least_common_multiple(scale, denominator_magnitude(value)) + least_common_multiple(scale, value.denom()) }) } -/// Return the positive magnitude of a validated non-zero denominator. -fn denominator_magnitude(value: &BigRational) -> BigInt { - match value.denom().sign() { - Sign::Minus => -value.denom(), - Sign::Plus => value.denom().clone(), - Sign::NoSign => unreachable!("RationalMatrix and RationalVector validate denominators"), - } -} - -/// Convert a rational to an integer using a positive divisible scale. +/// Convert a canonical rational to an integer using a positive divisible scale. +/// +/// Matrix and vector construction already prove that the denominator is positive. fn integer_at_scale(value: &BigRational, scale: &BigInt) -> BigInt { - let denominator = denominator_magnitude(value); - let multiplier = scale / denominator; - let numerator = match value.denom().sign() { - Sign::Minus => -value.numer(), - Sign::Plus => value.numer().clone(), - Sign::NoSign => unreachable!("RationalMatrix and RationalVector validate denominators"), - }; - numerator * multiplier + value.numer() * (scale / value.denom()) } /// Return the positive least common multiple of two positive integers. -fn least_common_multiple(lhs: BigInt, rhs: BigInt) -> BigInt { +fn least_common_multiple(lhs: BigInt, rhs: &BigInt) -> BigInt { let gcd = greatest_common_divisor(lhs.clone(), rhs.clone()); (lhs / gcd) * rhs } diff --git a/tests/proptest_factorizations.rs b/tests/proptest_factorizations.rs index aac6c5f..2d13cae 100644 --- a/tests/proptest_factorizations.rs +++ b/tests/proptest_factorizations.rs @@ -5,6 +5,8 @@ //! These tests construct matrices from known factors so we have a reliable oracle for //! determinant and solve behavior. +use core::{array::from_fn, cmp::Ordering}; + use approx::assert_abs_diff_eq; use pastey::paste; use proptest::{array, prelude::*}; @@ -34,6 +36,62 @@ fn nonzero_diag_entry() -> impl Strategy { prop_oneof![(-20i16..=-1i16), (1i16..=20i16)].prop_map(|x| f64::from(x) / 10.0) } +fn unit_lower(raw: &[[f64; D]; D]) -> [[f64; D]; D] { + from_fn(|i| { + from_fn(|j| match i.cmp(&j) { + Ordering::Equal => 1.0, + Ordering::Greater => raw[i][j], + Ordering::Less => 0.0, + }) + }) +} + +/// Construct A = L * U with unit-lower L and the supplied diagonal of U. +fn lu_product( + l_raw: &[[f64; D]; D], + u_raw: &[[f64; D]; D], + u_diag: &[f64; D], +) -> [[f64; D]; D] { + let l = unit_lower(l_raw); + let u: [[f64; D]; D] = from_fn(|i| { + from_fn(|j| match i.cmp(&j) { + Ordering::Equal => u_diag[i], + Ordering::Less => u_raw[i][j], + Ordering::Greater => 0.0, + }) + }); + from_fn(|i| { + from_fn(|j| { + // L[i][k] is zero for k > i; U[k][j] is zero for k > j. + (0..=i.min(j)).fold(0.0, |sum, k| l[i][k].mul_add(u[k][j], sum)) + }) + }) +} + +fn matvec(rows: &[[f64; D]; D], vector: &[f64; D]) -> [f64; D] { + from_fn(|i| { + rows[i] + .iter() + .zip(vector) + .fold(0.0, |sum, (&coefficient, &value)| { + coefficient.mul_add(value, sum) + }) + }) +} + +fn check_lu(rows: &[[f64; D]; D], expected_det: f64, x_true: &[f64; D]) { + let b_arr = matvec(rows, x_true); + let a = Matrix::try_from_rows(*rows).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); + assert_abs_diff_eq!(lu.det().unwrap(), expected_det, epsilon = 1e-8); + + let b = Vector::try_new(b_arr).unwrap(); + let x = lu.solve(b).unwrap().into_array(); + for (actual, expected) in x.iter().zip(x_true) { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-8); + } +} + macro_rules! gen_factorization_proptests { ($d:literal) => { paste! { @@ -49,18 +107,7 @@ macro_rules! gen_factorization_proptests { x_true in array::[](small_f64()), ) { // Construct A = L * diag(D) * L^T, where L is unit-lower-triangular. - let mut l = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - l[i][j] = if i == j { - 1.0 - } else if i > j { - l_raw[i][j] - } else { - 0.0 - }; - } - } + let l = unit_lower(&l_raw); let mut a_rows = [[0.0f64; $d]; $d]; for i in 0..$d { @@ -83,14 +130,7 @@ macro_rules! gen_factorization_proptests { acc }; - let mut b_arr = [0.0f64; $d]; - for i in 0..$d { - let mut sum = 0.0; - for j in 0..$d { - sum = a_rows[i][j].mul_add(x_true[j], sum); - } - b_arr[i] = sum; - } + let b_arr = matvec(&a_rows, &x_true); let a = Matrix::<$d>::try_from_rows(a_rows).unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); @@ -118,73 +158,9 @@ macro_rules! gen_factorization_proptests { u_diag in array::[](nonzero_diag_entry()), x_true in array::[](small_f64()), ) { - // Construct A = L * U, where L is unit-lower-triangular and U is upper-triangular. - let mut l = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - l[i][j] = if i == j { - 1.0 - } else if i > j { - l_raw[i][j] - } else { - 0.0 - }; - } - } - - let mut u = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - u[i][j] = if i == j { - u_diag[i] - } else if i < j { - u_raw[i][j] - } else { - 0.0 - }; - } - } - - let mut a_rows = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - let mut sum = 0.0; - // L[i][k] is zero for k > i; U[k][j] is zero for k > j. - let k_max = if i < j { i } else { j }; - for k in 0..=k_max { - sum = l[i][k].mul_add(u[k][j], sum); - } - a_rows[i][j] = sum; - } - } - - let expected_det = { - let mut acc = 1.0; - for i in 0..$d { - acc *= u_diag[i]; - } - acc - }; - - let mut b_arr = [0.0f64; $d]; - for i in 0..$d { - let mut sum = 0.0; - for j in 0..$d { - sum = a_rows[i][j].mul_add(x_true[j], sum); - } - b_arr[i] = sum; - } - - let a = Matrix::<$d>::try_from_rows(a_rows).unwrap(); - let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); - - assert_abs_diff_eq!(lu.det().unwrap(), expected_det, epsilon = 1e-8); - - let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = lu.solve(b).unwrap().into_array(); - for i in 0..$d { - assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); - } + let a_rows = lu_product(&l_raw, &u_raw, &u_diag); + let expected_det = u_diag.iter().fold(1.0, |acc, &value| acc * value); + check_lu(&a_rows, expected_det, &x_true); } #[test] @@ -198,77 +174,11 @@ macro_rules! gen_factorization_proptests { u_diag in array::[](nonzero_diag_entry()), x_true in array::[](small_f64()), ) { - // Construct A = P^{-1} * L * U, where P swaps the first two rows. - // This ensures det(A) has an extra sign flip vs det(LU). - let mut l = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - l[i][j] = if i == j { - 1.0 - } else if i > j { - l_raw[i][j] - } else { - 0.0 - }; - } - } - - let mut u = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - u[i][j] = if i == j { - u_diag[i] - } else if i < j { - u_raw[i][j] - } else { - 0.0 - }; - } - } - - let mut lu_rows = [[0.0f64; $d]; $d]; - for i in 0..$d { - for j in 0..$d { - let mut sum = 0.0; - let k_max = if i < j { i } else { j }; - for k in 0..=k_max { - sum = l[i][k].mul_add(u[k][j], sum); - } - lu_rows[i][j] = sum; - } - } - - // Apply P^{-1}: swap rows 0 and 1. - let mut a_rows = lu_rows; + // A = P^{-1} * L * U: swapping rows 0 and 1 flips det(LU)'s sign. + let mut a_rows = lu_product(&l_raw, &u_raw, &u_diag); a_rows.swap(0, 1); - - let expected_det = { - let mut acc = 1.0; - for i in 0..$d { - acc *= u_diag[i]; - } - -acc - }; - - let mut b_arr = [0.0f64; $d]; - for i in 0..$d { - let mut sum = 0.0; - for j in 0..$d { - sum = a_rows[i][j].mul_add(x_true[j], sum); - } - b_arr[i] = sum; - } - - let a = Matrix::<$d>::try_from_rows(a_rows).unwrap(); - let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); - - assert_abs_diff_eq!(lu.det().unwrap(), expected_det, epsilon = 1e-8); - - let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = lu.solve(b).unwrap().into_array(); - for i in 0..$d { - assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); - } + let expected_det = -u_diag.iter().fold(1.0, |acc, &value| acc * value); + check_lu(&a_rows, expected_det, &x_true); } } } diff --git a/tests/proptest_gram.rs b/tests/proptest_gram.rs new file mode 100644 index 0000000..bac7ef3 --- /dev/null +++ b/tests/proptest_gram.rs @@ -0,0 +1,213 @@ +#![forbid(unsafe_code)] + +//! Independent integer matrix-product oracles for Gram construction. + +use core::array::from_fn; + +use approx::assert_abs_diff_eq; +use pastey::paste; +use proptest::prelude::*; + +use la_stack::prelude::*; + +#[path = "common/proptest_config.rs"] +mod proptest_config; +use proptest_config::with_default_cases; + +fn check_product(entries: &[i16]) { + let integers: [[i16; N]; M] = from_fn(|i| from_fn(|k| entries[i * N + k])); + let vectors = integers.map(|row| Vector::try_new(row.map(f64::from)).unwrap()); + let result = gram_matrix(&vectors).unwrap(); + for i in 0..M { + for j in 0..M { + // Separate integer multiplication and summation, not the FMA kernel. + let expected: i32 = (0..N) + .map(|k| i32::from(integers[i][k]) * i32::from(integers[j][k])) + .sum(); + assert_abs_diff_eq!(result.as_rows()[i][j], f64::from(expected), epsilon = 0.0); + assert_eq!( + result.as_rows()[i][j].to_bits(), + result.as_rows()[j][i].to_bits() + ); + } + } +} + +fn check_orthogonal_and_dependent() { + let vectors: [Vector; D] = + from_fn(|i| Vector::try_new(from_fn(|j| if i == j { 2.0 } else { 0.0 })).unwrap()); + let gram = gram_matrix(&vectors).unwrap(); + for i in 0..D { + for j in 0..D { + assert_abs_diff_eq!( + gram.as_rows()[i][j], + if i == j { 4.0 } else { 0.0 }, + epsilon = 0.0 + ); + } + } + let expected = (0..D).fold(1.0, |product, _| product * 4.0); + assert_abs_diff_eq!( + gram.ldlt(Tolerance::try_new(0.0).unwrap()) + .unwrap() + .det() + .unwrap(), + expected, + epsilon = 0.0 + ); + let dependent = [vectors[0]; D]; + let dependent_gram = gram_matrix(&dependent).unwrap(); + assert_eq!(dependent_gram.as_rows(), &[[4.0; D]; D]); + assert_eq!( + dependent_gram + .ldlt(Tolerance::try_new(0.0).unwrap()) + .unwrap_err(), + LaError::singular_numerical(1, FactorizationKind::Ldlt, 0.0, 0.0) + ); +} + +// A separately rounded product or reversed reduction loses the exact residual. +fn check_fused_reduction() { + let delta = 2.0_f64.powi(-27); + let left: Vector = Vector::try_new(from_fn(|k| match k { + 0 => 1.0, + 1 => 1.0 + delta, + _ => 0.0, + })) + .unwrap(); + let right = Vector::try_new(from_fn(|k| match k { + 0 => -1.0, + 1 => 1.0 - delta, + _ => 0.0, + })) + .unwrap(); + let gram = gram_matrix(&[left, right]).unwrap(); + // -1 + (1 + 2^-27)(1 - 2^-27) = -2^-54, exactly representable. + let expected = (-2.0_f64.powi(-54)).to_bits(); + assert_eq!(gram.as_rows()[0][1].to_bits(), expected); + assert_eq!(gram.as_rows()[1][0].to_bits(), expected); +} + +macro_rules! gram_cases { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(with_default_cases(64))] + #[test] + fn [](entries in prop::collection::vec(-100i16..=100, 64)) { + check_product::<$d, $d>(&entries); + if $d != 8 { + check_product::<$d, 8>(&entries); + check_product::<8, $d>(&entries); + } + } + } + + #[test] + fn []() { + check_orthogonal_and_dependent::<$d>(); + } + + #[test] + fn []() { + check_fused_reduction::<$d>(); + } + } + }; +} + +gram_cases!(2); +gram_cases!(3); +gram_cases!(4); +gram_cases!(5); +gram_cases!(8); + +#[test] +fn empty_dimensions_and_const_evaluation() { + const EMPTY: Result, LaError> = gram_matrix::<0, 8>(&[]); + const ZERO: Result, LaError> = gram_matrix(&[Vector::<0>::zero(); 3]); + assert_eq!(EMPTY.unwrap().as_rows(), &[] as &[[f64; 0]; 0]); + assert_eq!(ZERO.unwrap().as_rows(), &[[0.0; 3]; 3]); + assert_eq!( + gram_matrix::<0, 0>(&[]).unwrap().as_rows(), + &[] as &[[f64; 0]; 0] + ); +} + +#[test] +fn nonempty_const_evaluation() { + const VECTOR: Vector<2> = match Vector::try_new([3.0, 4.0]) { + Ok(vector) => vector, + Err(_) => panic!("finite literal input"), + }; + const GRAM: Result, LaError> = gram_matrix(&[VECTOR]); + assert_eq!(GRAM.unwrap().as_rows(), &[[25.0]]); +} + +#[test] +fn underflow_preserves_mirrored_negative_zero() { + let vectors = [ + Vector::try_new([2.0_f64.powi(-537)]).unwrap(), + Vector::try_new([-2.0_f64.powi(-538)]).unwrap(), + ]; + let gram = gram_matrix(&vectors).unwrap(); + assert_eq!(gram.as_rows()[0][0].to_bits(), 1); // Smallest positive subnormal. + assert_eq!(gram.as_rows()[0][1].to_bits(), (-0.0_f64).to_bits()); + assert_eq!(gram.as_rows()[1][0].to_bits(), (-0.0_f64).to_bits()); + assert_eq!(gram.as_rows()[1][1].to_bits(), 0); +} + +#[test] +fn overflow_checks_off_diagonal_before_later_diagonal() { + let vectors = [ + Vector::try_new([0.0, 2.0]).unwrap(), + Vector::try_new([f64::MAX, f64::MAX]).unwrap(), + ]; + // (0,0) succeeds; (0,1) fails at coordinate 1. Evaluating (1,1) first + // would instead report coordinate 0, violating upper-triangle row order. + assert_eq!( + gram_matrix(&vectors).unwrap_err(), + LaError::non_finite_computation_step(ArithmeticOperation::VectorDotProduct, 1) + ); + + let later_diagonal = [ + Vector::try_new([1.0, 0.0]).unwrap(), + Vector::try_new([0.0, f64::MAX]).unwrap(), + ]; + assert_eq!( + gram_matrix(&later_diagonal).unwrap_err(), + LaError::non_finite_computation_step(ArithmeticOperation::VectorDotProduct, 1) + ); +} + +#[test] +fn mixed_scale_and_signed_zero() { + let vectors = [ + Vector::try_new([2.0_f64.powi(500), 2.0_f64.powi(-500), -0.0]).unwrap(), + Vector::try_new([2.0_f64.powi(-500), -2.0_f64.powi(500), 0.0]).unwrap(), + ]; + let gram = gram_matrix(&vectors).unwrap(); + assert_eq!( + gram.as_rows(), + &[[2.0_f64.powi(1000), 0.0], [0.0, 2.0_f64.powi(1000)]] + ); + assert_eq!( + gram.as_rows()[0][1].to_bits(), + gram.as_rows()[1][0].to_bits() + ); +} + +#[test] +fn overflow_preserves_dot_diagnostics() { + for (data, index) in [([f64::MAX, 0.0], 0), ([1.0e154, 1.0e154], 1)] { + let vectors = [Vector::try_new(data).unwrap()]; + assert!(matches!( + gram_matrix(&vectors), + Err(LaError::NonFinite { + origin: NonFiniteOrigin::Computation { operation: ArithmeticOperation::VectorDotProduct, .. }, + location: NonFiniteLocation::Step { index: actual, .. }, + .. + }) if actual == index + )); + } +}