feat: add const-generic Gram matrix construction

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

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<const D: usize>` for fixed-length `f64` vectors backed by `[f64; D]`
+- `gram_matrix(&[Vector<N>; M])` for allocation-free `Matrix<M>` 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<const D: usize>` for fixed-size square `f64` matrices backed by `[[f64; D]; D]`
 - `Interval` and `IntervalMatrix<const D: usize>` for outward-rounded,
   proof-bearing determinant filters through D=7

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)

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<const M: usize, const N: usize>(
+    vectors: &[Vector<N>; M],
+) -> Result<Matrix<M>, 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<const M: usize, const N: usize>(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();
+}

src/gram.rs

@@ -0,0 +1,73 @@
+#![forbid(unsafe_code)]
+
+//! Fixed-size Gram construction.
+
+use crate::{LaError, Matrix, Vector};
+
+/// Construct a stack-backed [`Matrix<M>`] 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<N>`] 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<const M: usize, const N: usize>(
+    vectors: &[Vector<N>; M],
+) -> Result<Matrix<M>, 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)
+}

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")]

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<Item = &'a BigRational>) -> 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
 }