docs(science): add mathematical basis and clarify numerical guarantees

AGENTS.md

@@ -430,7 +430,8 @@ When creating or updating issues:
   - `src/vector.rs`: `Vector<const D: usize>` (`[f64; D]`)
   - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `inf_norm`, `det`, `det_direct`)
   - `src/lu.rs`: `Lu<const D: usize>` factorization with partial pivoting (`solve`, `det`)
-  - `src/ldlt.rs`: `Ldlt<const D: usize>` factorization without pivoting for symmetric SPD/PSD matrices (`solve`, `det`)
+  - `src/ldlt.rs`: `Ldlt<const D: usize>` factorization without pivoting for exactly
+    symmetric positive-definite matrices (`solve`, `det`)
   - `src/exact.rs`: exact arithmetic behind `features = ["exact"]`:
     - Determinants: `det_exact()`, strict `det_exact_f64()`, rounded
       `det_exact_rounded_f64()`, and `det_sign_exact()` via a scaled `BigInt`

CITATION.cff

@@ -23,6 +23,7 @@ keywords:
   - "robust-predicates"
 abstract: >-
   la-stack is a Rust library providing fast, stack-allocated linear algebra primitives
-  and factorization routines (LU with partial pivoting and LDLT for symmetric SPD/PSD matrices)
+  and factorization routines (LU with partial pivoting and no-pivot LDLT for exactly
+  symmetric positive-definite matrices)
   for fixed small dimensions using const generics.
 license: BSD-3-Clause

CONTRIBUTING.md

@@ -69,6 +69,7 @@ Use the existing canonical documents instead of duplicating their guidance:
 |-------|---------------------|
 | Agent rules and repository invariants | [`AGENTS.md`](AGENTS.md) |
 | User-facing API, examples, and project scope | [`README.md`](README.md) |
+| Mathematical basis and numerical validity | [`docs/mathematical_basis.md`](docs/mathematical_basis.md) |
 | Package metadata, features, and dependencies | [`Cargo.toml`](Cargo.toml) |
 | Commands and validation workflow | [`justfile`](justfile), `just --list` |
 | Python support tooling | [`scripts/README.md`](scripts/README.md) |

REFERENCES.md

@@ -15,6 +15,7 @@ Tagged releases are archived on Zenodo under the all-versions concept DOI
 - CodeRabbit AI, Inc. "CodeRabbit." <https://coderabbit.ai/>.
 - KiloCode. "KiloCode AI Engineering Assistant." <https://kilocode.ai/>.
 - OpenAI. "ChatGPT." <https://openai.com/chatgpt>.
+- OpenAI. "Codex." <https://openai.com/codex/>.
 - Warp Dev, Inc. "WARP." <https://www.warp.dev/>.
 
 All AI-generated output was reviewed and/or edited by the maintainer.
@@ -62,50 +63,62 @@ matrix and RHS. After back-substitution, multiplying by the exact power-of-two s
 
 Both the determinant and solve paths convert their finite-by-construction entries via
 `decompose_proven_finite_f64`, which extracts the IEEE 754 binary64 sign, unbiased exponent,
-and significand [9] and strips trailing zeros from the significand so `|x| = m · 2^e` with
-`m` odd. The integer matrix is then assembled by shifting each mantissa left by
+and significand [9]. For nonzero `x`, it strips trailing zeros from the
+significand so `|x| = m · 2^e` with `m` odd; signed zeros use a separate zero
+component. The integer matrix is then assembled by shifting each mantissa left by
 `exp − e_min`, giving a GCD-free exact-integer starting point. Solves and D ≥ 5 determinants
 then apply Bareiss elimination; D ≤ 4 determinants use direct expansions. The test-only
 fallible wrapper `decompose_f64` verifies rejection of non-finite raw scalars, while the
 test-only `f64_to_big_rational` helper packages the same decomposition into a single
-`BigRational`. See Goldberg [10] for background on IEEE 754 representation and exact rational
-reconstruction.
+`BigRational`. See Goldberg [10] for background on floating-point representation and
+conversion.
 
-### LDL^T factorization (symmetric SPD/PSD)
+### LDLᵀ factorization (exactly symmetric positive-definite inputs)
 
-The LDL^T (often abbreviated "LDLT") implementation in `la-stack` is intended for symmetric positive
-definite (SPD) and positive semi-definite (PSD) matrices (e.g. Gram matrices), and does not perform
-pivoting.
+The no-pivot LDLT implementation targets `A = L D Lᵀ` for exactly symmetric
+positive-definite inputs [4-5, 11-12]. Successful construction requires every
+computed diagonal pivot to be positive and greater than the caller's tolerance.
+Computed zero and tolerance-small positive pivots are therefore part of the
+typed diagnostic domain, not returned in a usable factorization. Because the
+pivots are computed in binary64, a successful factorization is not an exact
+certificate that the represented matrix is positive definite.
 
-For background on the SPD/PSD setting, see [4-5]. For pivoted variants used for symmetric *indefinite*
-matrices, see [6].
+For pivoted variants used for symmetric *indefinite* matrices, see [6, 11-12].
 
 ### LU decomposition (Gaussian elimination with partial pivoting)
 
-The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization
-approach with partial pivoting for numerical stability.
+The LU implementation targets `P A = L U` and uses partial pivoting: each step
+selects the remaining entry of largest magnitude in the active column. Partial
+pivoting is a practical stability strategy, not an unconditional accuracy
+guarantee; worst-case growth and average-case behavior are distinct concerns.
 
-See references [1-3] below.
+See [1-3, 11-12] for stability analysis, finite-precision behavior, and standard
+algorithmic background.
 
 ## References
 
 1. Trefethen, Lloyd N., and Robert S. Schreiber. "Average-case stability of Gaussian elimination."
    *SIAM Journal on Matrix Analysis and Applications* 11.3 (1990): 335–360.
+   [DOI](https://doi.org/10.1137/0611023) ·
    [PDF](https://people.maths.ox.ac.uk/trefethen/publication/PDF/1990_44.pdf)
 2. Businger, P. A. "Monitoring the Numerical Stability of Gaussian Elimination."
-   *Numerische Mathematik* 16 (1970/71): 360–361.
-   [Full text](https://eudml.org/doc/132040)
+   *Numerische Mathematik* 16.4 (1971): 360–361.
+   [DOI](https://doi.org/10.1007/BF02165006) · [Full text](https://eudml.org/doc/132040)
 3. Huang, Han, and K. Tikhomirov. "Average-case analysis of the Gaussian elimination with partial pivoting."
    *Probability Theory and Related Fields* 189 (2024): 501–567.
-   [Open-access PDF](https://link.springer.com/article/10.1007/s00440-024-01276-2) (also: [arXiv:2206.01726](https://arxiv.org/abs/2206.01726))
-4. Cholesky, Andre-Louis. "On the numerical solution of systems of linear equations"
-   (manuscript dated 2 Dec 1910; published 2005).
-   Scan + English analysis: [BibNum](https://www.bibnum.education.fr/mathematiques/algebre/sur-la-resolution-numerique-des-systemes-d-equations-lineaires)
-5. Brezinski, Claude. "La methode de Cholesky." (2005).
-   [PDF](https://eudml.org/doc/252115)
-6. Bunch, J. R., L. Kaufman, and B. N. Parlett. "Decomposition of a Symmetric Matrix."
-   *Numerische Mathematik* 27 (1976/1977): 95–110.
-   [Full text](https://eudml.org/doc/132435)
+   [DOI](https://doi.org/10.1007/s00440-024-01276-2) ·
+   [Open-access article](https://link.springer.com/article/10.1007/s00440-024-01276-2) ·
+   [arXiv:2206.01726](https://arxiv.org/abs/2206.01726)
+4. Cholesky, André-Louis. "Sur la résolution numérique des systèmes d'équations linéaires."
+   *Bulletin de la Sabix* 39 (2005): 81–95. Manuscript dated 2 December 1910.
+   [DOI](https://doi.org/10.4000/sabix.529)
+5. Brezinski, Claude. "La méthode de Cholesky."
+   *Revue d'histoire des mathématiques* 11.2 (2005): 205–238.
+   [DOI](https://doi.org/10.24033/rhm.30) ·
+   [Full text](https://www.numdam.org/articles/10.24033/rhm.30/)
+6. Bunch, James R., Linda Kaufman, and Beresford N. Parlett. "Decomposition of a Symmetric Matrix."
+   *Numerische Mathematik* 27 (1976): 95–109.
+   [DOI](https://doi.org/10.1007/BF01399088) · [Full text](https://eudml.org/doc/132435)
 7. Bareiss, Erwin H. "Sylvester's Identity and Multistep Integer-Preserving Gaussian
    Elimination." *Mathematics of Computation* 22.103 (1968): 565–578.
    [DOI](https://doi.org/10.1090/S0025-5718-1968-0226829-0) ·
@@ -123,6 +136,15 @@ See references [1-3] below.
 10. Goldberg, David. "What Every Computer Scientist Should Know About Floating-Point
     Arithmetic." *ACM Computing Surveys* 23.1 (1991): 5–48.
     [DOI](https://doi.org/10.1145/103162.103163) ·
-    [PDF](https://www.validlab.com/goldberg/paper.pdf)
-    Comprehensive survey of IEEE 754 representation, rounding, and exact rational
-    reconstruction of floating-point values.
+    [Authorized HTML reprint](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)
+    Comprehensive survey of floating-point representation, rounding, and conversion.
+11. Higham, Nicholas J. *Accuracy and Stability of Numerical Algorithms*. 2nd ed.
+    Society for Industrial and Applied Mathematics, 2002.
+    [DOI](https://doi.org/10.1137/1.9780898718027)
+12. Golub, Gene H., and Charles F. Van Loan. *Matrix Computations*. 4th ed.
+    Johns Hopkins University Press, 2013.
+    [DOI](https://doi.org/10.56021/9781421407944) ·
+    [Publisher record](https://www.press.jhu.edu/books/title/10678/matrix-computations)
+13. Kalibera, Tomas, and Richard Jones. "Rigorous Benchmarking in Reasonable Time."
+    *Proceedings of the 2013 International Symposium on Memory Management* (ISMM '13),
+    2013: 63–74. [DOI](https://doi.org/10.1145/2464157.2464160)

benches/common/exact.rs

@@ -493,9 +493,9 @@ fn assert_approximate_determinant(actual: f64, exact: &BigRational, operation: &
 /// Validate the floating-point determinant operations used by Criterion.
 ///
 /// This runs during setup, outside timed closures. The deterministic `det` and
-/// `det_direct` results are compared with the independent Leibniz oracle; on
-/// current revisions the combined direct result is additionally checked
-/// against its certified absolute bound.
+/// `det_direct` results are compared with the independent Leibniz oracle, and
+/// `det_errbound` must certify the observed direct error. On current revisions,
+/// the paired result must also match the separate determinant and bound calls.
 ///
 /// # Panics
 ///
@@ -513,11 +513,24 @@ pub fn validate_f64_determinant_benchmarks<const D: usize>(input: &ValidatedExac
         .matrix()
         .det_direct()
         .or_abort("direct f64 determinant oracle check");
+    let standalone_bound = input
+        .matrix()
+        .det_errbound()
+        .or_abort("standalone determinant error-bound oracle check");
     if D <= 4 {
         let Some(direct) = direct else {
             panic!("det_direct must support benchmark dimension {D}");
         };
+        let Some(standalone_bound) = standalone_bound else {
+            panic!("det_errbound must support benchmark dimension {D}");
+        };
         assert_approximate_determinant(direct, &exact, "direct f64 determinant");
+        let observed_error = (rational_from_f64(direct) - &exact).abs();
+        let certified_bound = rational_from_f64(standalone_bound);
+        assert!(
+            observed_error <= certified_bound,
+            "direct determinant error {observed_error} exceeds standalone certified bound {certified_bound}",
+        );
 
         #[cfg(not(la_stack_v0_4_3_api))]
         {
@@ -529,19 +542,21 @@ pub fn validate_f64_determinant_benchmarks<const D: usize>(input: &ValidatedExac
                 panic!("the baseline fixture must have a certified D={D} determinant bound");
             };
             assert_eq!(estimate.determinant().to_bits(), direct.to_bits());
+            assert_eq!(
+                estimate.absolute_error_bound().to_bits(),
+                standalone_bound.to_bits(),
+            );
             assert!(
                 estimate.determinant().abs() > estimate.absolute_error_bound(),
                 "the headline D={D} det_sign_exact benchmark must exercise the fast filter",
             );
-            let observed_error = (rational_from_f64(direct) - &exact).abs();
-            let certified_bound = rational_from_f64(estimate.absolute_error_bound());
-            assert!(
-                observed_error <= certified_bound,
-                "direct determinant error {observed_error} exceeds certified bound {certified_bound}",
-            );
         }
     } else {
         assert!(direct.is_none(), "det_direct unexpectedly supports D={D}");
+        assert!(
+            standalone_bound.is_none(),
+            "det_errbound unexpectedly supports D={D}",
+        );
 
         #[cfg(not(la_stack_v0_4_3_api))]
         assert!(

benches/exact.rs

@@ -8,7 +8,8 @@
 //! 1. **General-case benches** (`exact_d{2..5}`) — a single
 //!    well-conditioned diagonally-dominant matrix per dimension.  These
 //!    measure typical-case performance and track regressions against a
-//!    reproducible input.
+//!    reproducible input. D=2..=4 also include the direct determinant and
+//!    certified error-bound f64 baselines.
 //! 2. **Adversarial / extreme-input benches** — matrices chosen to
 //!    stress specific corners of the exact-arithmetic pipeline:
 //!    near-singularity (forces the exact integer fallback), large f64 entries
@@ -200,9 +201,57 @@ fn bench_det_direct<const D: usize>(
     });
 }
 
-macro_rules! register_det_direct_benchmark {
+/// Add the standalone certified determinant-bound baseline.
+fn bench_det_errbound<const D: usize>(
+    group: &mut BenchmarkGroup<'_, WallTime>,
+    input: &ValidatedExactInput<D>,
+) {
+    let bound = input
+        .matrix()
+        .det_errbound()
+        .or_abort("determinant error-bound setup")
+        .or_abort("determinant error-bound setup");
+    black_box(bound);
+    group.bench_function("det_errbound", |bencher| {
+        bencher.iter(|| {
+            let bound = black_box(input.matrix())
+                .det_errbound()
+                .or_abort("f64 determinant error bound")
+                .or_abort("f64 determinant error bound");
+            black_box(bound);
+        });
+    });
+}
+
+/// Add the paired direct-determinant and certified-bound baseline.
+#[cfg(not(la_stack_v0_4_3_api))]
+fn bench_det_direct_with_errbound<const D: usize>(
+    group: &mut BenchmarkGroup<'_, WallTime>,
+    input: &ValidatedExactInput<D>,
+) {
+    let estimate = input
+        .matrix()
+        .det_direct_with_errbound()
+        .or_abort("paired determinant-bound setup")
+        .or_abort("paired determinant-bound setup");
+    black_box((estimate.determinant(), estimate.absolute_error_bound()));
+    group.bench_function("det_direct_with_errbound", |bencher| {
+        bencher.iter(|| {
+            let estimate = black_box(input.matrix())
+                .det_direct_with_errbound()
+                .or_abort("paired f64 determinant and error bound")
+                .or_abort("paired f64 determinant and error bound");
+            black_box((estimate.determinant(), estimate.absolute_error_bound()));
+        });
+    });
+}
+
+macro_rules! register_det_filter_benchmarks {
     ($group:expr, $input:expr, supported) => {{
         bench_det_direct(&mut $group, &$input);
+        #[cfg(not(la_stack_v0_4_3_api))]
+        bench_det_direct_with_errbound(&mut $group, &$input);
+        bench_det_errbound(&mut $group, &$input);
     }};
     ($group:expr, $matrix:expr, unsupported) => {};
 }
@@ -227,7 +276,7 @@ macro_rules! gen_exact_benches_for_dim {
             });
         });
 
-        register_det_direct_benchmark!(group, input, $direct);
+        register_det_filter_benchmarks!(group, input, $direct);
 
         for &operation in GENERAL_OPERATIONS {
             bench_exact_operation(&mut group, operation, &input);

benches/vs_linalg.rs

@@ -6,7 +6,8 @@
 //! fixed dimensions.
 //!
 //! Notes:
-//! - Determinant is benchmarked via LU on all sides (nalgebra uses closed-forms for 1×1/2×2/3×3).
+//! - Determinant groups distinguish factorization-inclusive LU, `Matrix::det`,
+//!   precomputed LU determinant queries, and precomputed LDLT/Cholesky queries.
 //! - Matrix infinity norm is the maximum absolute row sum on all sides.
 
 use std::hint::black_box;

docs/BENCHMARKING.md

@@ -56,8 +56,9 @@ factorization in the dependency version used here.
 **`exact`** (`benches/exact.rs`) measures exact-arithmetic methods
 (`det_exact`, `solve_exact`, `det_sign_exact`, strict `*_result` conversions,
 and lossy `*_rounded_f64` conversions) alongside the f64 `det` baseline across
-D=2-5 and `det_direct` across its supported D=2-4 range. Use this suite to
-understand exact-arithmetic cost and track optimization progress.
+D=2-5. Its supported D=2-4 range also includes `det_direct`, the paired
+`det_direct_with_errbound`, and the bound-only `det_errbound`. Use this suite
+to understand exact-arithmetic cost and track optimization progress.
 
 ## Common Workflows
 
@@ -91,6 +92,11 @@ at warning for both revisions because the current manifest's lint policy may
 reject historical source that predates a lint, even though that source remains
 valid benchmark input.
 
+The paired `det_direct_with_errbound` API postdates v0.4.3, so its D=2-4
+baseline rows remain explicitly unavailable in that comparison. The older
+bound-only `det_errbound` API is source-compatible and remains a required
+shared-harness baseline.
+
 The v0.4.3 LU/LDLT balanced-range determinant paths return an incorrect zero,
 so their two D=8 stress rows are deliberately not timed as baselines. Reports
 leave those baselines explicitly unavailable rather than presenting invalid
@@ -179,7 +185,7 @@ coverage or provenance aborts publication. Use `--allow-partial` only for
 exploratory CSV/SVG output; it cannot update README and its sidecar explicitly
 marks measurement provenance unavailable.
 
-See `scripts/criterion_dim_plot.py --help` for plotting options.
+See `uv run --locked criterion-dim-plot --help` for plotting options.
 
 ### Create The Release Performance Report
 
@@ -247,6 +253,9 @@ The README table uses `median.point_estimate` in nanoseconds. Lower is better,
 but point-estimate ratios alone are descriptive and do not establish a
 statistically supported performance difference. Preserve Criterion confidence
 intervals or repeat controlled runs when making a stronger claim.
+For experimental-design background on controlled repetitions and uncertainty,
+see [REFERENCES.md](../REFERENCES.md) \[13\]; these workflows do not claim to
+implement every recommendation in that study.
 
 All three crates receive equivalent deterministic inputs for a given dimension:
 
@@ -341,7 +350,8 @@ Before timing begins, every fixed, adversarial, and corpus input is consumed int
 a private-field `ValidatedExactInput` after checks by an independent exact
 oracle. Timed and registration helpers accept only that proof-bearing wrapper. A
 factorial-time Leibniz determinant over exact rational reconstructions verifies
-determinant values and signs; exact residuals verify `A x = b`; and
+determinant values and signs, including each direct determinant's certified
+absolute bound; exact residuals verify `A x = b`; and
 strict/rounded binary64 results are checked for their exact bits, typed reason,
 and first failing component. These checks run outside timed Criterion closures.
 Any disagreement or unexpected error fails setup instead of becoming an