feat(interval): add certified determinant sign enclosures

Cargo.toml

@@ -81,6 +81,11 @@ name = "exact"
 harness = false
 required-features = [ "bench", "exact" ]
 
+[[bench]]
+name = "interval"
+harness = false
+required-features = [ "bench" ]
+
 [profile.release]
 lto = "fat"
 codegen-units = 1

README.md

@@ -23,6 +23,8 @@ while keeping the API intentionally small and explicit.
 
 - `Vector<const D: usize>` for fixed-length `f64` vectors backed by `[f64; D]`
 - `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
 - `RationalVector<const D: usize>` and `RationalMatrix<const D: usize>` for
   exact rational inputs behind the optional `"exact"` feature
 - `Lu<const D: usize>` for LU factorization with partial pivoting (solve + det)
@@ -40,6 +42,13 @@ factorization tolerances are rejection thresholds, not accuracy guarantees. For
 D≤4, direct determinants can be paired with a conservative absolute roundoff
 bound when its range preconditions hold.
 
+Derived binary64 expressions can instead be assembled with `Interval`
+subtraction, addition, multiplication, negation, and square. The resulting
+`IntervalMatrix<D>` determinant sign is certified through D=7 when its enclosure
+separates zero; the singleton `[0, 0]` also certifies exact zero. Every other
+overlap with zero is explicitly inconclusive. This default-feature surface is
+distinct from arbitrary-precision exact arithmetic.
+
 With `features = ["exact"]`, callers can either lift stored binary64 inputs
 losslessly or supply already-exact rational inputs for exact determinant signs,
 determinant values, and solves. Exactness over binary64 input starts at the
@@ -61,6 +70,8 @@ for the algorithms, validity boundaries, and supporting references.
   semantics
 - ✅ Error-bounded f64 determinant filtering plus optional exact signs
   (`det_errbound`, `det_sign_exact`)
+- ✅ Outward-rounded interval expressions and division-free determinant signs
+  through D=7, with explicit inconclusive evidence
 - ✅ Exact determinant values and linear solves via optional arbitrary-precision
   arithmetic (`det_exact`, `solve_exact`, strict/rounded f64 conversions)
 - ✅ Explicit algorithms (LU, solve, determinant)
@@ -94,17 +105,24 @@ for current release planning.
 - You want explicit LU / LDLT / determinant APIs rather than a broad algebra toolkit
 - You need exact determinants, exact determinant signs, or exact linear solves
   for fixed-size systems
+- You need a cheap, sound interval filter for determinant expressions assembled
+  from rounded binary64 operations
 - Robust predicates matter for geometry-style workloads near degeneracy
 - You prefer a default build with no runtime dependencies
 
-## 🔢 Scalar types
+## 🔢 Scalar and bounded-value types
 
-The public scalar model deliberately has exactly two input domains:
+The public point-value scalar model deliberately has two input domains:
 
 - finite `f64` through `Matrix<D>` and `Vector<D>` for floating-point work;
 - arbitrary-precision `BigRational` through `RationalMatrix<D>` and
   `RationalVector<D>` behind the optional `"exact"` feature.
 
+`Interval` is a separate bounded-value layer over finite `f64` endpoints. It
+encloses exact-real values during a small set of outward-rounded operations and
+feeds `IntervalMatrix<D>` determinant proofs; it does not make `Matrix` generic
+over alternate scalars or provide a general interval package.
+
 This is not a generic scalar-parameterized API. Exact support intentionally
 covers the robustness-sensitive operations that require it: determinant sign,
 determinant value, and linear solve, followed by explicit strict or rounded
@@ -128,7 +146,8 @@ la-stack = "0.4.5"
 
 ### Feature flags
 
-- `default`: no runtime dependencies
+- `default`: no runtime dependencies; includes outward-rounded `Interval` and
+  `IntervalMatrix` APIs
 - `exact`: exact determinant signs, determinant values, and solves over stored
   `f64` values or caller-supplied `BigRational` inputs
 - `bench`: repository-development gate used only by benchmark targets and
@@ -269,6 +288,57 @@ non-zero pivot; it does not misreport that numerical failure as an exact zero.
 Use `lu()` directly when you need a different tolerance policy, and use the
 exact determinant APIs when exact singularity classification matters.
 
+## 📦 Outward-rounded interval determinants
+
+`Interval` encloses expression construction that has not yet been reduced to a
+single stored `f64`. Point intervals preserve finite binary64 values exactly;
+`try_from_subtraction`, `try_add`, `try_mul`, `negate`, and `try_square` enclose
+the corresponding exact-real operations. `IntervalMatrix<D>::det_sign()` then
+uses a division-free subset expansion through D=7, returning positive,
+negative, zero, or inconclusive evidence.
+
+```rust
+use la_stack::prelude::*;
+
+fn main() -> Result<(), LaError> {
+    // Relative coordinates and the lifted norm retain their construction error.
+    let x = Interval::try_from_subtraction(0.1, 0.0)?;
+    let y = Interval::try_from_subtraction(0.1, 0.0)?;
+    let z = Interval::try_from_subtraction(0.1, 0.0)?;
+    let lifted = x
+        .try_square()?
+        .try_add(&y.try_square()?)?
+        .try_add(&z.try_square()?)?;
+
+    let matrix = IntervalMatrix::<4>::from_rows([
+        [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE],
+        [x, y, z, lifted],
+    ]);
+    assert_eq!(
+        matrix.det_sign()?,
+        IntervalDeterminantSign::Negative,
+    );
+    Ok(())
+}
+```
+
+Every successful interval keeps finite ordered endpoints. Subnormal bounds are
+preserved, both signed zeros are treated as real zero and canonicalized to
+`+0.0`, and underflowed nonzero products widen toward the least subnormal value.
+If an exact result range cannot fit between finite binary64 endpoints, the
+operation returns `LaError::IntervalRangeExhausted` with its interval operation
+recorded in `ArithmeticOperation`.
+
+`Positive`, `Negative`, and `Zero` are proofs. `Inconclusive` only means that
+the determinant enclosure overlaps zero; it must not be converted to equality
+or singularity. A filtered-exact caller should rebuild the same derived
+expression with `RationalMatrix` and call `det_sign()` when the interval result
+is inconclusive or reports range failure. Lifting a finished `Matrix` with
+`IntervalMatrix::from_matrix` encloses its stored entries, but cannot recover
+rounding that occurred while those entries were assembled.
+
 ## 🔬 Exact arithmetic (`"exact"` feature)
 
 The default build has **zero runtime dependencies**. Enable the optional
@@ -567,6 +637,9 @@ out of the common prelude.
 |---|---|---|---|
 | `Vector<D>` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `norm2_sq` |
 | `Matrix<D>` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below |
+| `Interval` | Two finite ordered `f64` bounds | Outward-rounded exact-real enclosure | `try_new`, `point`, `try_from_subtraction`, `try_add`, `try_mul`, `negate`, `try_square` |
+| `IntervalMatrix<D>` | `[[Interval; D]; D]` | Division-free determinant enclosure and sign proof through D=7 | `from_rows`, `try_from_point_rows`, `from_matrix`, `det`, `det_sign` |
+| `IntervalDeterminantSign` | enum | Positive, negative, zero, or inconclusive determinant evidence | — |
 | `RationalVector<D>`¹ | `[BigRational; D]` | Exact rational right-hand side and solution | `try_new`, `try_from_fn`, `as_array`, `into_array`, `get` |
 | `RationalMatrix<D>`¹ | `[[BigRational; D]; D]` | Exact rational matrix for determinant and solve operations | `try_from_rows`, `try_from_fn`, `as_rows`, `det_sign`, `det`, `solve` |
 | `DeterminantWithErrorBound` | Opaque validated pair | Paired direct determinant and certified absolute bound | `determinant`, `absolute_error_bound` |
@@ -578,8 +651,15 @@ out of the common prelude.
 | `ExactF64Conversion`¹ | trait | Strict or explicitly rounded conversion of exact results to `f64` | `try_to_f64`, `to_rounded_f64` |
 
 `Matrix<D>` and `Vector<D>` use the intentional inline `f64` scalar model.
-The exact-feature rational types retain fixed-size outer arrays while their
-`BigRational` scalars use arbitrary-precision integer storage.
+`IntervalMatrix<D>` retains inline fixed-size storage and uses a fixed 128-entry
+stack workspace for its supported determinant dimensions. The exact-feature
+rational types retain fixed-size outer arrays while their `BigRational` scalars
+use arbitrary-precision integer storage.
+
+For a runtime-selected interval dimension from 0 through
+`MAX_INTERVAL_MATRIX_DIM` (7), `try_with_interval_matrix!` dispatches to a
+concrete `IntervalMatrix<N>`. This is useful when stable Rust cannot express a
+derived const dimension such as `D + 1`.
 
 For a runtime dimension from 0 through `MAX_STACK_MATRIX_DISPATCH_DIM` (7),
 `try_with_stack_matrix!` dispatches to a concrete `Matrix<N>` while preserving
@@ -620,6 +700,9 @@ observed pivot magnitude, and tolerance, while exact-arithmetic singularity is
 identified separately. `LaError::NonFinite` retains the crate-wide non-finite
 contract but uses `NonFiniteOrigin`, `NonFiniteLocation`, and
 `ArithmeticOperation` to distinguish invalid inputs from computed overflow.
+`LaError::InvertedInterval` preserves rejected finite bounds when the lower
+endpoint exceeds the upper endpoint, and `LaError::IntervalRangeExhausted`
+distinguishes finite-input interval range loss from a non-finite value.
 `InvalidToleranceReason` distinguishes negative from non-finite tolerances, and
 `PositiveSemidefiniteViolation` distinguishes negative LDLT pivots from a zero
 pivot with nonzero coupling. Match these public enums with a wildcard and use
@@ -665,6 +748,11 @@ references. Releases produced with the rational-input harness include Criterion
 point estimates and confidence intervals for these rows; comparisons against a
 pre-API baseline retain them as explicit current-only measurements.
 
+The focused `interval` Criterion suite covers conclusive and inconclusive
+relative-coordinate lifted determinant signs at D=4 and the maximum supported
+D=7 workload. Run it with `just bench-interval`; fixture validation stays
+outside the timed closures.
+
 <!-- BENCH_TABLE:lu_solve:median:new:BEGIN -->
 
 | D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | reduction vs nalgebra (point est.) | reduction vs faer (point est.) |

REFERENCES.md

@@ -23,6 +23,27 @@ No generated content was used without human oversight.
 
 ## Linear algebra algorithms
 
+### Outward-rounded interval determinant sign
+
+`Interval` uses IEEE-754 round-to-nearest binary64 operations plus adjacent
+representable values to enclose exact-real addition, subtraction,
+multiplication, and square results (references 9–11). Addition and subtraction
+use an error-free `TwoSum` residual [8]; multiplication independently compares
+the exact integer-significand product with the rounded binary64 result,
+including gradual underflow to zero. Results whose exact range cannot fit
+between finite binary64 endpoints return a typed range failure rather than
+storing infinity. For the broader standardized interval arithmetic model, see
+[14]; this crate does not claim IEEE 1788 conformance.
+
+`IntervalMatrix::det()` evaluates the Leibniz expansion with a division-free
+column-subset dynamic program. It uses `2^D` inline interval states and
+`D × 2^(D-1)` coefficient products through D=7. A determinant interval strictly
+separated from zero certifies its sign; `[0, 0]` certifies zero; every other
+overlap is explicitly inconclusive. The determinant identity is standard
+linear algebra (reference 12); the interval evaluation and subset-DP
+organization are implemented specifically for this crate's small
+fixed-dimension scope.
+
 ### Absolute error bound for closed-form determinants
 
 `Matrix::det_errbound()` returns a conservative Shewchuk-style absolute error bound [8]
@@ -150,3 +171,7 @@ algorithmic background.
 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)
+14. IEEE Computer Society. "IEEE Standard for Interval Arithmetic."
+    *IEEE Std 1788-2015*, 2015: 1–97.
+    [DOI](https://doi.org/10.1109/IEEESTD.2015.7140721) ·
+    [IEEE record](https://standards.ieee.org/ieee/1788/4431/)

benches/interval.rs

@@ -0,0 +1,127 @@
+#![forbid(unsafe_code)]
+
+//! Criterion coverage for conclusive and inconclusive interval determinants.
+
+use std::hint::black_box;
+
+use criterion::Criterion;
+
+use la_stack::{Interval, IntervalDeterminantSign, IntervalMatrix, LaError};
+
+#[path = "common/bench_utils.rs"]
+mod bench_utils;
+use bench_utils::OrAbort;
+
+/// Assemble the relative-coordinate lifted matrix for a tetrahedral in-sphere
+/// predicate whose interval determinant is conclusively negative.
+fn conclusive_lifted_4x4() -> Result<IntervalMatrix<4>, LaError> {
+    let x = Interval::try_from_subtraction(0.1, 0.0)?;
+    let y = Interval::try_from_subtraction(0.1, 0.0)?;
+    let z = Interval::try_from_subtraction(0.1, 0.0)?;
+    let lifted = x
+        .try_square()?
+        .try_add(&y.try_square()?)?
+        .try_add(&z.try_square()?)?;
+
+    Ok(IntervalMatrix::from_rows([
+        [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE],
+        [x, y, z, lifted],
+    ]))
+}
+
+/// Assemble a lifted boundary case whose final coefficient retains one ULP of
+/// expression uncertainty, forcing an inconclusive interval sign.
+fn inconclusive_lifted_4x4() -> Result<IntervalMatrix<4>, LaError> {
+    Ok(IntervalMatrix::from_rows([
+        [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE],
+        [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE],
+        [
+            Interval::ONE,
+            Interval::ONE,
+            Interval::ONE,
+            Interval::try_new(3.0_f64.next_down(), 3.0_f64.next_up())?,
+        ],
+    ]))
+}
+
+/// Assemble a six-coordinate lifted predicate matrix to exercise the maximum
+/// supported D=7 subset-DP workload.
+fn conclusive_lifted_7x7() -> Result<IntervalMatrix<7>, LaError> {
+    let mut matrix = IntervalMatrix::zero();
+    for index in 0..6 {
+        matrix.set(index, index, Interval::ONE)?;
+        matrix.set(index, 6, Interval::ONE)?;
+    }
+
+    let relative = Interval::try_from_subtraction(0.125, 0.0)?;
+    let square = relative.try_square()?;
+    let mut lifted = Interval::ZERO;
+    for column in 0..6 {
+        matrix.set(6, column, relative)?;
+        lifted = lifted.try_add(&square)?;
+    }
+    matrix.set(6, 6, lifted)?;
+    Ok(matrix)
+}
+
+fn main() {
+    let conclusive_4 =
+        conclusive_lifted_4x4().or_abort("conclusive D=4 interval fixture construction");
+    let inconclusive_4 =
+        inconclusive_lifted_4x4().or_abort("inconclusive D=4 interval fixture construction");
+    let conclusive_7 =
+        conclusive_lifted_7x7().or_abort("conclusive D=7 interval fixture construction");
+
+    assert_eq!(
+        conclusive_4
+            .det_sign()
+            .or_abort("conclusive D=4 interval fixture validation"),
+        IntervalDeterminantSign::Negative,
+    );
+    assert_eq!(
+        inconclusive_4
+            .det_sign()
+            .or_abort("inconclusive D=4 interval fixture validation"),
+        IntervalDeterminantSign::Inconclusive,
+    );
+    assert_eq!(
+        conclusive_7
+            .det_sign()
+            .or_abort("conclusive D=7 interval fixture validation"),
+        IntervalDeterminantSign::Negative,
+    );
+
+    let mut criterion = Criterion::default().configure_from_args();
+    {
+        let mut group = criterion.benchmark_group("interval_det_sign");
+        group.bench_function("d4_conclusive_lifted", |bencher| {
+            bencher.iter(|| {
+                let sign = black_box(&conclusive_4)
+                    .det_sign()
+                    .or_abort("D=4 conclusive interval determinant");
+                let _ = black_box(sign);
+            });
+        });
+        group.bench_function("d4_inconclusive_lifted", |bencher| {
+            bencher.iter(|| {
+                let sign = black_box(&inconclusive_4)
+                    .det_sign()
+                    .or_abort("D=4 inconclusive interval determinant");
+                let _ = black_box(sign);
+            });
+        });
+        group.bench_function("d7_conclusive_lifted", |bencher| {
+            bencher.iter(|| {
+                let sign = black_box(&conclusive_7)
+                    .det_sign()
+                    .or_abort("D=7 conclusive interval determinant");
+                let _ = black_box(sign);
+            });
+        });
+        group.finish();
+    }
+    criterion.final_summary();
+}

docs/BENCHMARKING.md

@@ -28,6 +28,7 @@ the commands measure and where their outputs go.
 | Non-exact release-signal check against tags | `just performance-local-non-exact v0.4.5 v0.4.4` |
 | Fast saved-baseline loop | `just bench-save-baseline <name> <suite>` then `just bench-compare <name> <suite> all-benches` |
 | Full crate comparison | `just bench-vs-linalg` |
+| Interval determinant filter | `just bench-interval` |
 | README table and plot | `just performance-release` then `just performance-readme` |
 | Release report | `just performance-release v0.4.5 v0.4.4` |
 | Build docs from retained release inputs | `just performance-doc` |
@@ -63,7 +64,7 @@ promotion in one command.
 
 ## Benchmark Suites
 
-`la-stack` has two Criterion benchmark suites.
+`la-stack` has three Criterion benchmark suites.
 
 Newly rendered reports use one table per selected suite. Dimension and
 adversarial-input group appear in a `Case` column instead of creating a separate
@@ -88,6 +89,14 @@ suite compares row-cleared Bareiss operations with direct `BigRational` Gaussian
 operations over already-exact rational inputs across D=2-8. Use it to understand
 exact-arithmetic cost and track optimization progress.
 
+**`interval`** (`benches/interval.rs`) measures the default-feature,
+division-free interval determinant sign filter. Its fixtures cover a conclusive
+4×4 relative-coordinate lifted predicate, the corresponding inconclusive
+boundary regime, and a conclusive 7×7 lifted workload at the supported dimension
+limit. Fixture construction and expected-sign validation occur outside the
+timed closures. This suite is a focused kernel signal; it is not part of the
+release-to-release `vs_linalg` or `exact` report schema.
+
 ## Common Workflows
 
 ### Compare Current Code With The Latest Release