From 9c0e3aede88d877e734cc6396c264695dc59e92a Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 04:46:48 +0900
Subject: [PATCH 01/19] feat(core): add covariance standardization contract

---
 .../src/covariance_standardization.rs         | 264 ++++++++++++++++++
 1 file changed, 264 insertions(+)
 create mode 100644 crates/mlsirm-core/src/covariance_standardization.rs

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
new file mode 100644
index 000000000..de4e458bd
--- /dev/null
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -0,0 +1,264 @@
+//! Domain-neutral covariance-to-correlation standardization.
+//!
+//! This module owns only reusable static numerical normalization. Product- or
+//! study-specific temporal admission rules (for example, requiring event time)
+//! belong to the consuming bounded context and must wrap this contract rather
+//! than being encoded here.
+//!
+//! For a covariance matrix `Σ` with strictly positive diagonal `D`, the
+//! standardized matrix is `R = D^{-1/2} Σ D^{-1/2}`. The scalar specialization
+//! is therefore `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`.
+//! The implementation evaluates divisions sequentially to avoid overflowing a
+//! representable correlation by first forming `sqrt(v_i) * sqrt(v_j)`.
+//!
+//! The TEPP migration that motivated this owner contract concerns ctsem's
+//! `TIPREDVARstd`, but ctsem names, clocks, state equations, and event semantics
+//! deliberately do not appear in this reusable kernel.
+//!
+//! # Research provenance
+//!
+//! Driver, C. C., Oud, J. H. L., & Voelkle, M. C. (2017). Continuous time
+//! structural equation modeling with R package ctsem. *Journal of Statistical
+//! Software, 77*(5), 1–35. https://doi.org/10.18637/jss.v077.i05
+//!
+//! The ctsem source and paper provide the motivating covariance-standardization
+//! use case; the matrix identity implemented here is the ordinary definition
+//! converting a covariance matrix to its correlation matrix.
+
+use std::error::Error;
+use std::fmt::{Display, Formatter};
+
+/// Versioned public contract for reusable covariance standardization.
+pub const COVARIANCE_STANDARDIZATION_CONTRACT_VERSION: &str =
+    "fast_mlsirm.covariance_standardization@1.0.0";
+
+/// Fail-closed input and arithmetic errors for covariance standardization.
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+pub enum CovarianceStandardizationError {
+    /// Matrix dimension is zero, overflows `usize`, or does not match the slice.
+    InvalidShape,
+    /// At least one matrix entry or scalar variance is NaN or infinite.
+    NonFiniteInput,
+    /// A variance on the diagonal is zero or negative and cannot be standardized.
+    NonPositiveVariance,
+    /// Mirrored covariance cells disagree beyond floating-point tolerance.
+    NonSymmetricCovariance,
+    /// A pairwise covariance implies an absolute correlation materially above one.
+    InvalidPairwiseCovariance,
+    /// A finite input produced a non-finite standardized result.
+    NonFiniteResult,
+}
+
+impl Display for CovarianceStandardizationError {
+    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
+        let message = match self {
+            Self::InvalidShape => "covariance matrix shape is invalid",
+            Self::NonFiniteInput => "covariance input must be finite",
+            Self::NonPositiveVariance => "covariance diagonal must be strictly positive",
+            Self::NonSymmetricCovariance => "covariance matrix must be symmetric",
+            Self::InvalidPairwiseCovariance => {
+                "covariance pair violates the correlation magnitude bound"
+            }
+            Self::NonFiniteResult => "covariance standardization produced a non-finite result",
+        };
+        formatter.write_str(message)
+    }
+}
+
+impl Error for CovarianceStandardizationError {}
+
+/// Standardize one finite, strictly positive variance against itself.
+///
+/// The arithmetic is evaluated rather than replaced with a hard-coded `1.0` so
+/// this scalar reference exercises the same normalization contract consumed by
+/// matrix standardization and downstream parity tests.
+///
+/// # Errors
+///
+/// Returns [`CovarianceStandardizationError::NonFiniteInput`] for NaN or
+/// infinity, [`CovarianceStandardizationError::NonPositiveVariance`] for zero
+/// or a negative value, and [`CovarianceStandardizationError::NonFiniteResult`]
+/// if the arithmetic cannot produce a finite value.
+pub fn standardize_variance(
+    variance: f64,
+) -> Result<f64, CovarianceStandardizationError> {
+    if !variance.is_finite() {
+        return Err(CovarianceStandardizationError::NonFiniteInput);
+    }
+    if variance <= 0.0 {
+        return Err(CovarianceStandardizationError::NonPositiveVariance);
+    }
+
+    let inverse_sd = 1.0 / variance.sqrt();
+    let standardized = (variance * inverse_sd) * inverse_sd;
+    if !standardized.is_finite() {
+        return Err(CovarianceStandardizationError::NonFiniteResult);
+    }
+    Ok(standardized)
+}
+
+/// Convert a finite symmetric covariance matrix to a correlation matrix.
+///
+/// `covariance` is row-major with shape `dimension × dimension`. Every
+/// diagonal variance must be strictly positive. Symmetry is checked with a
+/// scale-aware binary64 tolerance; off-diagonal pairs that only differ by
+/// rounding are averaged safely (`a/2 + b/2`) before standardization so the
+/// returned matrix is exactly symmetric. Pairwise correlations whose absolute
+/// value exceeds one by more than floating-point tolerance fail closed.
+///
+/// This routine validates the pairwise covariance bounds but does not claim a
+/// full positive-semidefinite proof. A caller that requires PSD admission must
+/// apply that model-specific invariant separately.
+///
+/// # Errors
+///
+/// Returns a typed error for invalid shape, non-finite input, non-positive
+/// diagonal variance, material asymmetry, an impossible pairwise covariance,
+/// or non-finite output arithmetic.
+pub fn standardize_covariance_matrix(
+    covariance: &[f64],
+    dimension: usize,
+) -> Result<Vec<f64>, CovarianceStandardizationError> {
+    let expected_len = dimension
+        .checked_mul(dimension)
+        .ok_or(CovarianceStandardizationError::InvalidShape)?;
+    if dimension == 0 || covariance.len() != expected_len {
+        return Err(CovarianceStandardizationError::InvalidShape);
+    }
+    if covariance.iter().any(|value| !value.is_finite()) {
+        return Err(CovarianceStandardizationError::NonFiniteInput);
+    }
+
+    let mut standard_deviations = Vec::with_capacity(dimension);
+    for index in 0..dimension {
+        let variance = covariance[index * dimension + index];
+        if variance <= 0.0 {
+            return Err(CovarianceStandardizationError::NonPositiveVariance);
+        }
+        standard_deviations.push(variance.sqrt());
+    }
+
+    let mut correlation = vec![0.0; expected_len];
+    for index in 0..dimension {
+        correlation[index * dimension + index] = standardize_variance(
+            covariance[index * dimension + index],
+        )?;
+    }
+
+    for row in 0..dimension {
+        for column in (row + 1)..dimension {
+            let upper = covariance[row * dimension + column];
+            let lower = covariance[column * dimension + row];
+            let symmetry_scale = upper.abs().max(lower.abs()).max(1.0);
+            let symmetry_tolerance = 64.0 * f64::EPSILON * symmetry_scale;
+            if (upper - lower).abs() > symmetry_tolerance {
+                return Err(CovarianceStandardizationError::NonSymmetricCovariance);
+            }
+
+            let symmetric_covariance = upper * 0.5 + lower * 0.5;
+            let standardized =
+                (symmetric_covariance / standard_deviations[row]) / standard_deviations[column];
+            if !standardized.is_finite() {
+                return Err(CovarianceStandardizationError::NonFiniteResult);
+            }
+            let bound_tolerance = 128.0 * f64::EPSILON;
+            if standardized.abs() > 1.0 + bound_tolerance {
+                return Err(CovarianceStandardizationError::InvalidPairwiseCovariance);
+            }
+            let bounded = standardized.clamp(-1.0, 1.0);
+            correlation[row * dimension + column] = bounded;
+            correlation[column * dimension + row] = bounded;
+        }
+    }
+
+    Ok(correlation)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::{
+        CovarianceStandardizationError, standardize_covariance_matrix, standardize_variance,
+    };
+
+    fn assert_close(actual: f64, expected: f64, tolerance: f64) {
+        assert!(
+            (actual - expected).abs() <= tolerance,
+            "actual={actual:?} expected={expected:?} tolerance={tolerance:?}"
+        );
+    }
+
+    #[test]
+    fn scalar_reference_recovers_one_across_positive_scales() {
+        for variance in [
+            f64::MIN_POSITIVE,
+            1.0e-200,
+            0.25,
+            1.0,
+            6.4,
+            1.0e200,
+            f64::MAX,
+        ] {
+            assert_close(standardize_variance(variance).expect("positive variance"), 1.0, 8.0e-15);
+        }
+    }
+
+    #[test]
+    fn scalar_reference_fails_closed_for_invalid_variance() {
+        for variance in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
+            assert!(standardize_variance(variance).is_err());
+        }
+        assert_eq!(
+            standardize_variance(0.0),
+            Err(CovarianceStandardizationError::NonPositiveVariance)
+        );
+    }
+
+    #[test]
+    fn matrix_standardization_recovers_expected_correlation() {
+        let covariance = [4.0, 2.0, 2.0, 9.0];
+        let correlation = standardize_covariance_matrix(&covariance, 2).expect("covariance");
+        assert_close(correlation[0], 1.0, 8.0e-15);
+        assert_close(correlation[1], 1.0 / 3.0, 8.0e-15);
+        assert_close(correlation[2], 1.0 / 3.0, 8.0e-15);
+        assert_close(correlation[3], 1.0, 8.0e-15);
+    }
+
+    #[test]
+    fn matrix_standardization_is_scale_invariant() {
+        let base = [4.0, -3.0, -3.0, 9.0];
+        let scaled = [4.0e100, -3.0e100, -3.0e100, 9.0e100];
+        let base_result = standardize_covariance_matrix(&base, 2).expect("base");
+        let scaled_result = standardize_covariance_matrix(&scaled, 2).expect("scaled");
+        for (left, right) in base_result.iter().zip(scaled_result.iter()) {
+            assert_close(*left, *right, 8.0e-15);
+        }
+    }
+
+    #[test]
+    fn matrix_standardization_fails_closed_for_shape_and_numeric_defects() {
+        assert_eq!(
+            standardize_covariance_matrix(&[], 0),
+            Err(CovarianceStandardizationError::InvalidShape)
+        );
+        assert_eq!(
+            standardize_covariance_matrix(&[1.0, 0.0, 0.0], 2),
+            Err(CovarianceStandardizationError::InvalidShape)
+        );
+        assert_eq!(
+            standardize_covariance_matrix(&[1.0, f64::NAN, f64::NAN, 1.0], 2),
+            Err(CovarianceStandardizationError::NonFiniteInput)
+        );
+        assert_eq!(
+            standardize_covariance_matrix(&[0.0], 1),
+            Err(CovarianceStandardizationError::NonPositiveVariance)
+        );
+        assert_eq!(
+            standardize_covariance_matrix(&[1.0, 0.2, 0.3, 1.0], 2),
+            Err(CovarianceStandardizationError::NonSymmetricCovariance)
+        );
+        assert_eq!(
+            standardize_covariance_matrix(&[1.0, 1.1, 1.1, 1.0], 2),
+            Err(CovarianceStandardizationError::InvalidPairwiseCovariance)
+        );
+    }
+}

From 300b19fff61b9dd8b6475f89f9970b2da09d8774 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 04:47:00 +0900
Subject: [PATCH 02/19] feat(core): expose covariance standardization contract

---
 crates/mlsirm-core/src/entrypoint.rs | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/crates/mlsirm-core/src/entrypoint.rs b/crates/mlsirm-core/src/entrypoint.rs
index e18996c53..2122eff39 100644
--- a/crates/mlsirm-core/src/entrypoint.rs
+++ b/crates/mlsirm-core/src/entrypoint.rs
@@ -1,11 +1,12 @@
 //! Crate entrypoint that preserves the established core while adding modular rotation.
 //!
-//! The historical core remains in `lib.rs`. Keeping the rotation implementation
-//! in its own module boundary avoids adding another several-thousand-line
-//! section to that file and gives future CPU/GPU backends a stable public home.
+//! The historical core remains in `lib.rs`. Keeping newer implementations in
+//! explicit module boundaries avoids adding another several-thousand-line
+//! section to that file and gives reusable numerical contracts stable public homes.
 
 include!("lib.rs");
 
+pub mod covariance_standardization;
 pub mod interaction_map;
 pub mod rotation;
 pub mod sampling_design;

From a323e4de5ac1b8131ab5f272e6c33930f612777b Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 04:47:25 +0900
Subject: [PATCH 03/19] docs: trace covariance standardization owner contract

---
 ...variance-standardization-owner-contract.md | 43 +++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 docs/papers/covariance-standardization-owner-contract.md

diff --git a/docs/papers/covariance-standardization-owner-contract.md b/docs/papers/covariance-standardization-owner-contract.md
new file mode 100644
index 000000000..501b6f494
--- /dev/null
+++ b/docs/papers/covariance-standardization-owner-contract.md
@@ -0,0 +1,43 @@
+# Covariance standardization owner contract
+
+## Scope
+
+`mlsirm-core` owns the reusable static numerical map from a covariance matrix to its correlation matrix:
+
+\[
+R = D^{-1/2}\Sigma D^{-1/2},
+\]
+
+where `D` is the diagonal of strictly positive marginal variances. The scalar self-standardization is `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`.
+
+This contract is intentionally domain-neutral. It does not decide event time, valid time, knowledge cutoff, state evolution, temporal identification, or product-specific model activation. Those policies belong to consuming bounded contexts such as TEPP Longitudinal Modeling.
+
+## Triggering reuse case
+
+TEPP PR #475 recovered the ctsem `TIPREDVARstd` scalar map inside `psychometric_core::event_time`. Review of that implementation showed that the EventTime admission rule is TEPP-owned, while the arithmetic itself is reusable covariance standardization and therefore belongs in fast-mlsirm. The owner implementation is tracked by issue #1720.
+
+The TEPP evidence is preserved rather than copied as a product name here: Driver, Oud, and Voelkle describe the continuous-time model and the time-independent predictor covariance family; the 2017-era ctsem summary implementation forms standardized covariance quantities using inverse marginal standard deviations. fast-mlsirm exposes the generic normalization only. TEPP can later bind `TIPREDVARstd` through an ACL after a released version exists.
+
+## Numerical contract
+
+`fast_mlsirm.covariance_standardization@1.0.0` requires:
+
+- finite covariance cells;
+- a non-empty square row-major matrix;
+- strictly positive diagonal variances;
+- symmetry within an explicit binary64 tolerance;
+- pairwise covariance magnitude consistent with `|r| <= 1` up to floating-point tolerance.
+
+The implementation divides sequentially by each marginal standard deviation rather than first multiplying two square roots. This avoids overflowing a representable correlation solely because `sqrt(v_i) * sqrt(v_j)` exceeds the finite binary64 range. Near-boundary correlations within rounding tolerance are clamped to `[-1, 1]`; materially invalid pairwise covariance fails closed.
+
+The contract does not claim to prove full positive semidefiniteness. Model-specific PSD admission remains a separate invariant.
+
+## Recovery and parity
+
+Unit/property-style fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, asymmetry, and pairwise covariance beyond the correlation bound.
+
+Before TEPP removes its local duplicate, the TEPP adapter must prove parity with its preserved `TIPREDVARstd` fixtures while continuing to enforce EventTime semantics outside this kernel.
+
+## Reference
+
+Driver, C. C., Oud, J. H. L., & Voelkle, M. C. (2017). Continuous time structural equation modeling with R package ctsem. *Journal of Statistical Software, 77*(5), 1–35. https://doi.org/10.18637/jss.v077.i05

From 69ae2c40c4a459f9e946f9dc0850019dde06b25a Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 04:50:17 +0900
Subject: [PATCH 04/19] test(core): reject heuristic covariance admission

---
 .../covariance_standardization_contract.rs    | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 crates/mlsirm-core/tests/covariance_standardization_contract.rs

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
new file mode 100644
index 000000000..2d46df93c
--- /dev/null
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -0,0 +1,25 @@
+use mlsirm_core::covariance_standardization::{
+    standardize_covariance_matrix, CovarianceStandardizationError,
+};
+
+#[test]
+fn covariance_matrix_requires_exact_symmetry() {
+    let epsilon_offset = 8.0 * f64::EPSILON;
+    let covariance = [1.0, 0.25, 0.25 + epsilon_offset, 1.0];
+
+    assert_eq!(
+        standardize_covariance_matrix(&covariance, 2),
+        Err(CovarianceStandardizationError::NonSymmetricCovariance),
+    );
+}
+
+#[test]
+fn covariance_matrix_never_clamps_an_out_of_range_correlation() {
+    let above_one = 1.0 + 64.0 * f64::EPSILON;
+    let covariance = [1.0, above_one, above_one, 1.0];
+
+    assert_eq!(
+        standardize_covariance_matrix(&covariance, 2),
+        Err(CovarianceStandardizationError::InvalidPairwiseCovariance),
+    );
+}

From 6ed74a7d37492976133bd5450efcf0c1d6ecaa80 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 04:52:55 +0900
Subject: [PATCH 05/19] fix(core): remove heuristic covariance admission

---
 .../src/covariance_standardization.rs         | 38 +++++++++----------
 1 file changed, 17 insertions(+), 21 deletions(-)

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
index de4e458bd..883955364 100644
--- a/crates/mlsirm-core/src/covariance_standardization.rs
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -8,8 +8,9 @@
 //! For a covariance matrix `Σ` with strictly positive diagonal `D`, the
 //! standardized matrix is `R = D^{-1/2} Σ D^{-1/2}`. The scalar specialization
 //! is therefore `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`.
-//! The implementation evaluates divisions sequentially to avoid overflowing a
-//! representable correlation by first forming `sqrt(v_i) * sqrt(v_j)`.
+//! Matrix entries are divided by each marginal standard deviation sequentially,
+//! so the implementation does not form `sqrt(v_i) * sqrt(v_j)`, whose product
+//! could overflow even when the standardized correlation is representable.
 //!
 //! The TEPP migration that motivated this owner contract concerns ctsem's
 //! `TIPREDVARstd`, but ctsem names, clocks, state equations, and event semantics
@@ -41,9 +42,9 @@ pub enum CovarianceStandardizationError {
     NonFiniteInput,
     /// A variance on the diagonal is zero or negative and cannot be standardized.
     NonPositiveVariance,
-    /// Mirrored covariance cells disagree beyond floating-point tolerance.
+    /// Mirrored covariance cells are not exactly equal in binary64.
     NonSymmetricCovariance,
-    /// A pairwise covariance implies an absolute correlation materially above one.
+    /// A pairwise covariance implies an absolute correlation above one.
     InvalidPairwiseCovariance,
     /// A finite input produced a non-finite standardized result.
     NonFiniteResult,
@@ -100,11 +101,11 @@ pub fn standardize_variance(
 /// Convert a finite symmetric covariance matrix to a correlation matrix.
 ///
 /// `covariance` is row-major with shape `dimension × dimension`. Every
-/// diagonal variance must be strictly positive. Symmetry is checked with a
-/// scale-aware binary64 tolerance; off-diagonal pairs that only differ by
-/// rounding are averaged safely (`a/2 + b/2`) before standardization so the
-/// returned matrix is exactly symmetric. Pairwise correlations whose absolute
-/// value exceeds one by more than floating-point tolerance fail closed.
+/// diagonal variance must be strictly positive. Mirrored off-diagonal cells
+/// must be exactly equal in binary64. This contract does not invent a
+/// floating-point tolerance or clamp an out-of-range correlation into the
+/// admissible interval; callers that need approximate-symmetry preprocessing
+/// must perform and document that operation before calling this kernel.
 ///
 /// This routine validates the pairwise covariance bounds but does not claim a
 /// full positive-semidefinite proof. A caller that requires PSD admission must
@@ -113,8 +114,8 @@ pub fn standardize_variance(
 /// # Errors
 ///
 /// Returns a typed error for invalid shape, non-finite input, non-positive
-/// diagonal variance, material asymmetry, an impossible pairwise covariance,
-/// or non-finite output arithmetic.
+/// diagonal variance, asymmetric mirrored cells, an impossible pairwise
+/// covariance, or non-finite output arithmetic.
 pub fn standardize_covariance_matrix(
     covariance: &[f64],
     dimension: usize,
@@ -149,25 +150,20 @@ pub fn standardize_covariance_matrix(
         for column in (row + 1)..dimension {
             let upper = covariance[row * dimension + column];
             let lower = covariance[column * dimension + row];
-            let symmetry_scale = upper.abs().max(lower.abs()).max(1.0);
-            let symmetry_tolerance = 64.0 * f64::EPSILON * symmetry_scale;
-            if (upper - lower).abs() > symmetry_tolerance {
+            if upper != lower {
                 return Err(CovarianceStandardizationError::NonSymmetricCovariance);
             }
 
-            let symmetric_covariance = upper * 0.5 + lower * 0.5;
             let standardized =
-                (symmetric_covariance / standard_deviations[row]) / standard_deviations[column];
+                (upper / standard_deviations[row]) / standard_deviations[column];
             if !standardized.is_finite() {
                 return Err(CovarianceStandardizationError::NonFiniteResult);
             }
-            let bound_tolerance = 128.0 * f64::EPSILON;
-            if standardized.abs() > 1.0 + bound_tolerance {
+            if standardized.abs() > 1.0 {
                 return Err(CovarianceStandardizationError::InvalidPairwiseCovariance);
             }
-            let bounded = standardized.clamp(-1.0, 1.0);
-            correlation[row * dimension + column] = bounded;
-            correlation[column * dimension + row] = bounded;
+            correlation[row * dimension + column] = standardized;
+            correlation[column * dimension + row] = standardized;
         }
     }
 

From feaf80fe1246c0ad125f5770635cfb780a416aaf Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:03:49 +0900
Subject: [PATCH 06/19] test(core): pin covariance contract identity

---
 .../tests/covariance_standardization_contract.rs         | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index 2d46df93c..74a5ef8b9 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -1,7 +1,16 @@
 use mlsirm_core::covariance_standardization::{
     standardize_covariance_matrix, CovarianceStandardizationError,
+    COVARIANCE_STANDARDIZATION_CONTRACT_VERSION,
 };
 
+#[test]
+fn covariance_standardization_contract_version_is_stable() {
+    assert_eq!(
+        COVARIANCE_STANDARDIZATION_CONTRACT_VERSION,
+        "fast_mlsirm.covariance_standardization@1.0.0",
+    );
+}
+
 #[test]
 fn covariance_matrix_requires_exact_symmetry() {
     let epsilon_offset = 8.0 * f64::EPSILON;

From 13d2159ae15903b26b5f002c9caa583db5b110b1 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:04:38 +0900
Subject: [PATCH 07/19] test(core): preserve public covariance contract
 coverage

---
 .../covariance_standardization_contract.rs    | 91 ++++++++++++++++++-
 1 file changed, 89 insertions(+), 2 deletions(-)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index 74a5ef8b9..e0c649d0a 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -1,14 +1,101 @@
 use mlsirm_core::covariance_standardization::{
-    standardize_covariance_matrix, CovarianceStandardizationError,
+    standardize_covariance_matrix, standardize_variance, CovarianceStandardizationError,
     COVARIANCE_STANDARDIZATION_CONTRACT_VERSION,
 };
 
+fn assert_close(actual: f64, expected: f64, tolerance: f64) {
+    assert!(
+        (actual - expected).abs() <= tolerance,
+        "actual={actual:?} expected={expected:?} tolerance={tolerance:?}",
+    );
+}
+
 #[test]
-fn covariance_standardization_contract_version_is_stable() {
+fn public_covariance_standardization_contract_is_versioned_and_scale_invariant() {
     assert_eq!(
         COVARIANCE_STANDARDIZATION_CONTRACT_VERSION,
         "fast_mlsirm.covariance_standardization@1.0.0",
     );
+
+    for variance in [
+        f64::MIN_POSITIVE,
+        1.0e-200,
+        0.25,
+        1.0,
+        6.4,
+        1.0e200,
+        f64::MAX,
+    ] {
+        assert_close(
+            standardize_variance(variance).expect("strictly positive finite variance"),
+            1.0,
+            8.0e-15,
+        );
+    }
+}
+
+#[test]
+fn public_scalar_contract_fails_closed_for_unstandardizable_variance() {
+    assert_eq!(
+        standardize_variance(0.0),
+        Err(CovarianceStandardizationError::NonPositiveVariance),
+    );
+    assert_eq!(
+        standardize_variance(-1.0),
+        Err(CovarianceStandardizationError::NonPositiveVariance),
+    );
+    for variance in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
+        assert_eq!(
+            standardize_variance(variance),
+            Err(CovarianceStandardizationError::NonFiniteInput),
+        );
+    }
+}
+
+#[test]
+fn public_matrix_contract_recovers_known_correlation_and_scale_invariance() {
+    let base = [4.0, -3.0, -3.0, 9.0];
+    let scaled = [4.0e100, -3.0e100, -3.0e100, 9.0e100];
+    let base_result = standardize_covariance_matrix(&base, 2).expect("base covariance");
+    let scaled_result = standardize_covariance_matrix(&scaled, 2).expect("scaled covariance");
+    let expected = [1.0, -0.5, -0.5, 1.0];
+
+    for ((base_value, scaled_value), expected_value) in base_result
+        .iter()
+        .zip(scaled_result.iter())
+        .zip(expected.iter())
+    {
+        assert_close(*base_value, *expected_value, 8.0e-15);
+        assert_close(*scaled_value, *expected_value, 8.0e-15);
+    }
+}
+
+#[test]
+fn public_matrix_contract_fails_closed_for_invalid_covariance() {
+    assert_eq!(
+        standardize_covariance_matrix(&[], 0),
+        Err(CovarianceStandardizationError::InvalidShape),
+    );
+    assert_eq!(
+        standardize_covariance_matrix(&[1.0, 0.0, 0.0], 2),
+        Err(CovarianceStandardizationError::InvalidShape),
+    );
+    assert_eq!(
+        standardize_covariance_matrix(&[1.0, f64::NAN, f64::NAN, 1.0], 2),
+        Err(CovarianceStandardizationError::NonFiniteInput),
+    );
+    assert_eq!(
+        standardize_covariance_matrix(&[0.0], 1),
+        Err(CovarianceStandardizationError::NonPositiveVariance),
+    );
+    assert_eq!(
+        standardize_covariance_matrix(&[1.0, 0.2, 0.3, 1.0], 2),
+        Err(CovarianceStandardizationError::NonSymmetricCovariance),
+    );
+    assert_eq!(
+        standardize_covariance_matrix(&[1.0, 1.1, 1.1, 1.0], 2),
+        Err(CovarianceStandardizationError::InvalidPairwiseCovariance),
+    );
 }
 
 #[test]

From 82d4362082ac50e8987f96b731126b8e68896d8e Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:07:59 +0900
Subject: [PATCH 08/19] test(core): reproduce valid boundary covariance
 roundoff

---
 .../tests/covariance_standardization_contract.rs  | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index e0c649d0a..96c6d3b70 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -98,6 +98,21 @@ fn public_matrix_contract_fails_closed_for_invalid_covariance() {
     );
 }
 
+#[test]
+fn valid_boundary_covariance_survives_one_ulp_standardization_roundoff() {
+    // These finite binary64 values satisfy c² <= v1*v2 exactly as represented,
+    // while sequential f64 division rounds the raw ratio to next_up(1.0).
+    let variance_one = 2.0533691813403163e-253;
+    let variance_two = 3.3250793271932294e47;
+    let covariance = 2.612970611386659e-103;
+    let matrix = [variance_one, covariance, covariance, variance_two];
+
+    let correlation = standardize_covariance_matrix(&matrix, 2)
+        .expect("an exactly admissible binary64 covariance must not be rejected");
+    assert_eq!(correlation[1], 1.0);
+    assert_eq!(correlation[2], 1.0);
+}
+
 #[test]
 fn covariance_matrix_requires_exact_symmetry() {
     let epsilon_offset = 8.0 * f64::EPSILON;

From e7a9ebf1a35162af60b30c836965fa5fc0afa6c5 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:11:22 +0900
Subject: [PATCH 09/19] fix(core): admit exact-valid covariance boundaries

---
 .../src/covariance_standardization.rs         | 137 ++++++++++++++++--
 1 file changed, 122 insertions(+), 15 deletions(-)

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
index 883955364..640427c5d 100644
--- a/crates/mlsirm-core/src/covariance_standardization.rs
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -26,6 +26,7 @@
 //! use case; the matrix identity implemented here is the ordinary definition
 //! converting a covariance matrix to its correlation matrix.
 
+use std::cmp::Ordering;
 use std::error::Error;
 use std::fmt::{Display, Formatter};
 
@@ -44,7 +45,7 @@ pub enum CovarianceStandardizationError {
     NonPositiveVariance,
     /// Mirrored covariance cells are not exactly equal in binary64.
     NonSymmetricCovariance,
-    /// A pairwise covariance implies an absolute correlation above one.
+    /// A represented covariance pair violates `c² <= v_i * v_j` exactly.
     InvalidPairwiseCovariance,
     /// A finite input produced a non-finite standardized result.
     NonFiniteResult,
@@ -68,6 +69,78 @@ impl Display for CovarianceStandardizationError {
 
 impl Error for CovarianceStandardizationError {}
 
+/// Return an exact integer-significand representation of one finite binary64 value.
+///
+/// The returned pair `(significand, exponent)` satisfies
+/// `abs(value) = significand * 2^exponent` without floating-point rounding.
+fn binary64_components(value: f64) -> (u64, i32) {
+    let bits = value.abs().to_bits();
+    let exponent_bits = ((bits >> 52) & 0x7ff) as i32;
+    let fraction = bits & ((1_u64 << 52) - 1);
+    if exponent_bits == 0 {
+        (fraction, -1074)
+    } else {
+        ((1_u64 << 52) | fraction, exponent_bits - 1023 - 52)
+    }
+}
+
+/// Compare two non-negative exact integers multiplied by powers of two.
+fn scaled_integer_le(
+    left_significand: u128,
+    left_exponent: i32,
+    right_significand: u128,
+    right_exponent: i32,
+) -> bool {
+    if left_significand == 0 {
+        return true;
+    }
+    if right_significand == 0 {
+        return false;
+    }
+
+    let left_bits = 128_i32 - left_significand.leading_zeros() as i32;
+    let right_bits = 128_i32 - right_significand.leading_zeros() as i32;
+    let left_top_bit = left_exponent + left_bits - 1;
+    let right_top_bit = right_exponent + right_bits - 1;
+
+    match left_top_bit.cmp(&right_top_bit) {
+        Ordering::Less => true,
+        Ordering::Greater => false,
+        Ordering::Equal => {
+            if left_exponent >= right_exponent {
+                let shift = (left_exponent - right_exponent) as u32;
+                left_significand
+                    .checked_shl(shift)
+                    .is_some_and(|aligned| aligned <= right_significand)
+            } else {
+                let shift = (right_exponent - left_exponent) as u32;
+                right_significand
+                    .checked_shl(shift)
+                    .is_some_and(|aligned| left_significand <= aligned)
+            }
+        }
+    }
+}
+
+/// Check the covariance Cauchy-Schwarz bound exactly for represented binary64 inputs.
+fn pairwise_covariance_is_admissible(covariance: f64, variance_one: f64, variance_two: f64) -> bool {
+    let (covariance_significand, covariance_exponent) = binary64_components(covariance);
+    let (variance_one_significand, variance_one_exponent) = binary64_components(variance_one);
+    let (variance_two_significand, variance_two_exponent) = binary64_components(variance_two);
+
+    let covariance_square =
+        u128::from(covariance_significand) * u128::from(covariance_significand);
+    let variance_product =
+        u128::from(variance_one_significand) * u128::from(variance_two_significand);
+
+    scaled_integer_le(
+        covariance_square,
+        covariance_exponent * 2,
+        variance_product,
+        variance_one_exponent + variance_two_exponent,
+    )
+}
+
 /// Standardize one finite, strictly positive variance against itself.
 ///
 /// The arithmetic is evaluated rather than replaced with a hard-coded `1.0` so
@@ -102,10 +175,17 @@ pub fn standardize_variance(
 ///
 /// `covariance` is row-major with shape `dimension × dimension`. Every
 /// diagonal variance must be strictly positive. Mirrored off-diagonal cells
-/// must be exactly equal in binary64. This contract does not invent a
-/// floating-point tolerance or clamp an out-of-range correlation into the
-/// admissible interval; callers that need approximate-symmetry preprocessing
-/// must perform and document that operation before calling this kernel.
+/// must be exactly equal in binary64. Pairwise admissibility is decided from
+/// the exact represented binary64 integers using `c² <= v_i * v_j`, so an
+/// invalid covariance is never accepted merely because floating-point division
+/// rounded its correlation back into range.
+///
+/// After exact admission, sequential division can round a mathematically valid
+/// boundary correlation just outside `[-1, 1]`. Only then is the numerical
+/// result projected back to that mathematically certified interval. This is a
+/// consequence of the exact bound, not an empirical epsilon or tolerance.
+/// Callers that need approximate-symmetry preprocessing must define and validate
+/// that policy explicitly before calling this kernel.
 ///
 /// This routine validates the pairwise covariance bounds but does not claim a
 /// full positive-semidefinite proof. A caller that requires PSD admission must
@@ -141,9 +221,8 @@ pub fn standardize_covariance_matrix(
 
     let mut correlation = vec![0.0; expected_len];
     for index in 0..dimension {
-        correlation[index * dimension + index] = standardize_variance(
-            covariance[index * dimension + index],
-        )?;
+        correlation[index * dimension + index] =
+            standardize_variance(covariance[index * dimension + index])?;
     }
 
     for row in 0..dimension {
@@ -154,16 +233,20 @@ pub fn standardize_covariance_matrix(
                 return Err(CovarianceStandardizationError::NonSymmetricCovariance);
             }
 
+            let variance_row = covariance[row * dimension + row];
+            let variance_column = covariance[column * dimension + column];
+            if !pairwise_covariance_is_admissible(upper, variance_row, variance_column) {
+                return Err(CovarianceStandardizationError::InvalidPairwiseCovariance);
+            }
+
             let standardized =
                 (upper / standard_deviations[row]) / standard_deviations[column];
             if !standardized.is_finite() {
                 return Err(CovarianceStandardizationError::NonFiniteResult);
             }
-            if standardized.abs() > 1.0 {
-                return Err(CovarianceStandardizationError::InvalidPairwiseCovariance);
-            }
-            correlation[row * dimension + column] = standardized;
-            correlation[column * dimension + row] = standardized;
+            let bounded = standardized.clamp(-1.0, 1.0);
+            correlation[row * dimension + column] = bounded;
+            correlation[column * dimension + row] = bounded;
         }
     }
 
@@ -173,7 +256,8 @@ pub fn standardize_covariance_matrix(
 #[cfg(test)]
 mod tests {
     use super::{
-        CovarianceStandardizationError, standardize_covariance_matrix, standardize_variance,
+        scaled_integer_le, CovarianceStandardizationError, standardize_covariance_matrix,
+        standardize_variance,
     };
 
     fn assert_close(actual: f64, expected: f64, tolerance: f64) {
@@ -183,6 +267,16 @@ mod tests {
         );
     }
 
+    #[test]
+    fn exact_scaled_integer_comparison_handles_zero_magnitude_and_alignment() {
+        assert!(scaled_integer_le(0, -100, 1, -1000));
+        assert!(!scaled_integer_le(1, 0, 0, 0));
+        assert!(scaled_integer_le(1, 0, 1, 1));
+        assert!(!scaled_integer_le(1, 1, 1, 0));
+        assert!(scaled_integer_le(1, 1, 2, 0));
+        assert!(scaled_integer_le(2, 0, 1, 1));
+    }
+
     #[test]
     fn scalar_reference_recovers_one_across_positive_scales() {
         for variance in [
@@ -194,7 +288,11 @@ mod tests {
             1.0e200,
             f64::MAX,
         ] {
-            assert_close(standardize_variance(variance).expect("positive variance"), 1.0, 8.0e-15);
+            assert_close(
+                standardize_variance(variance).expect("positive variance"),
+                1.0,
+                8.0e-15,
+            );
         }
     }
 
@@ -230,6 +328,15 @@ mod tests {
         }
     }
 
+    #[test]
+    fn matrix_standardization_accepts_zero_covariance_with_subnormal_variance() {
+        let smallest_subnormal = f64::from_bits(1);
+        let covariance = [smallest_subnormal, 0.0, 0.0, 1.0];
+        let correlation = standardize_covariance_matrix(&covariance, 2).expect("valid covariance");
+        assert_eq!(correlation[1], 0.0);
+        assert_eq!(correlation[2], 0.0);
+    }
+
     #[test]
     fn matrix_standardization_fails_closed_for_shape_and_numeric_defects() {
         assert_eq!(

From e8278ac7925c38b1b7921f686b4da12b7937a539 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:12:03 +0900
Subject: [PATCH 10/19] docs(core): specify exact covariance boundary policy

---
 .../covariance-standardization-owner-contract.md     | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/docs/papers/covariance-standardization-owner-contract.md b/docs/papers/covariance-standardization-owner-contract.md
index 501b6f494..913892f59 100644
--- a/docs/papers/covariance-standardization-owner-contract.md
+++ b/docs/papers/covariance-standardization-owner-contract.md
@@ -25,16 +25,18 @@ The TEPP evidence is preserved rather than copied as a product name here: Driver
 - finite covariance cells;
 - a non-empty square row-major matrix;
 - strictly positive diagonal variances;
-- symmetry within an explicit binary64 tolerance;
-- pairwise covariance magnitude consistent with `|r| <= 1` up to floating-point tolerance.
+- mirrored covariance cells that are exactly equal in the supplied binary64 representation;
+- exact pairwise admissibility of the represented values under `c² <= v_i * v_j`.
 
-The implementation divides sequentially by each marginal standard deviation rather than first multiplying two square roots. This avoids overflowing a representable correlation solely because `sqrt(v_i) * sqrt(v_j)` exceeds the finite binary64 range. Near-boundary correlations within rounding tolerance are clamped to `[-1, 1]`; materially invalid pairwise covariance fails closed.
+Pairwise admission does not use an empirical epsilon. Each finite binary64 value is decomposed into its exact integer significand and power-of-two exponent, and the two products in `c² <= v_i * v_j` are compared exactly with `u128` significand products plus exponent alignment. Consequently, a genuinely invalid represented covariance fails closed even if later floating-point division would round its correlation back into range.
 
-The contract does not claim to prove full positive semidefiniteness. Model-specific PSD admission remains a separate invariant.
+The correlation itself is evaluated by sequentially dividing the covariance by each marginal standard deviation rather than first multiplying the two square roots. This avoids overflowing a representable result solely because `sqrt(v_i) * sqrt(v_j)` exceeds the finite binary64 range. Sequential division can still round a mathematically admissible boundary correlation one representable value outside `[-1, 1]`. After the exact admissibility proof, and only after that proof, the computed value is projected back to `[-1, 1]`. This projection is therefore a consequence of the exact represented-input bound rather than a tolerance for invalid covariance.
+
+The contract does not claim to prove full positive semidefiniteness. Model-specific PSD admission remains a separate invariant. Callers that wish to accept approximately symmetric observations must define and validate that preprocessing policy before calling this exact low-level kernel.
 
 ## Recovery and parity
 
-Unit/property-style fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, asymmetry, and pairwise covariance beyond the correlation bound.
+Unit and public-contract fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, exact-symmetry rejection, pairwise covariance beyond the exact correlation bound, zero covariance with a subnormal variance, and a finite exact-valid boundary case whose sequential binary64 divisions round to `next_up(1.0)` before the bound-certified projection.
 
 Before TEPP removes its local duplicate, the TEPP adapter must prove parity with its preserved `TIPREDVARstd` fixtures while continuing to enforce EventTime semantics outside this kernel.
 

From a15fba524c74ac336827068afe610f3818d8a3aa Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:16:59 +0900
Subject: [PATCH 11/19] test(core): reproduce variance-order underflow

---
 .../tests/covariance_standardization_contract.rs | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index 96c6d3b70..c9386ab33 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -113,6 +113,22 @@ fn valid_boundary_covariance_survives_one_ulp_standardization_roundoff() {
     assert_eq!(correlation[2], 1.0);
 }
 
+#[test]
+fn extreme_variance_ordering_preserves_nonzero_correlation() {
+    let covariance = 1.0e-200;
+    let large_first = [f64::MAX, covariance, covariance, f64::from_bits(1)];
+    let small_first = [f64::from_bits(1), covariance, covariance, f64::MAX];
+
+    let large_first_result =
+        standardize_covariance_matrix(&large_first, 2).expect("valid covariance");
+    let small_first_result =
+        standardize_covariance_matrix(&small_first, 2).expect("permuted valid covariance");
+
+    assert!(large_first_result[1] > 0.0);
+    assert_eq!(large_first_result[1], small_first_result[1]);
+    assert_eq!(large_first_result[2], small_first_result[2]);
+}
+
 #[test]
 fn covariance_matrix_requires_exact_symmetry() {
     let epsilon_offset = 8.0 * f64::EPSILON;

From 6bda7010f2c9796037544446e26b591f2d4980fe Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:19:30 +0900
Subject: [PATCH 12/19] fix(core): normalize by smaller marginal scale first

---
 .../src/covariance_standardization.rs         | 43 ++++++++++++-------
 1 file changed, 28 insertions(+), 15 deletions(-)

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
index 640427c5d..019166886 100644
--- a/crates/mlsirm-core/src/covariance_standardization.rs
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -8,9 +8,10 @@
 //! For a covariance matrix `Σ` with strictly positive diagonal `D`, the
 //! standardized matrix is `R = D^{-1/2} Σ D^{-1/2}`. The scalar specialization
 //! is therefore `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`.
-//! Matrix entries are divided by each marginal standard deviation sequentially,
-//! so the implementation does not form `sqrt(v_i) * sqrt(v_j)`, whose product
-//! could overflow even when the standardized correlation is representable.
+//! Matrix entries are divided first by the smaller marginal standard deviation
+//! and then by the larger one. This avoids forming `sqrt(v_i) * sqrt(v_j)`,
+//! whose product can overflow, while also avoiding avoidable intermediate
+//! underflow when the two marginal scales differ by many orders of magnitude.
 //!
 //! The TEPP migration that motivated this owner contract concerns ctsem's
 //! `TIPREDVARstd`, but ctsem names, clocks, state equations, and event semantics
@@ -123,7 +124,11 @@ fn scaled_integer_le(
 }
 
 /// Check the covariance Cauchy-Schwarz bound exactly for represented binary64 inputs.
-fn pairwise_covariance_is_admissible(covariance: f64, variance_one: f64, variance_two: f64) -> bool {
+fn pairwise_covariance_is_admissible(
+    covariance: f64,
+    variance_one: f64,
+    variance_two: f64,
+) -> bool {
     let (covariance_significand, covariance_exponent) = binary64_components(covariance);
     let (variance_one_significand, variance_one_exponent) = binary64_components(variance_one);
     let (variance_two_significand, variance_two_exponent) = binary64_components(variance_two);
@@ -180,16 +185,19 @@ pub fn standardize_variance(
 /// invalid covariance is never accepted merely because floating-point division
 /// rounded its correlation back into range.
 ///
-/// After exact admission, sequential division can round a mathematically valid
-/// boundary correlation just outside `[-1, 1]`. Only then is the numerical
-/// result projected back to that mathematically certified interval. This is a
-/// consequence of the exact bound, not an empirical epsilon or tolerance.
-/// Callers that need approximate-symmetry preprocessing must define and validate
-/// that policy explicitly before calling this kernel.
+/// After exact admission, the covariance is divided by the smaller marginal
+/// standard deviation before the larger one. For an admissible pair the first
+/// quotient is bounded in magnitude by the larger standard deviation, avoiding
+/// overflow while minimizing avoidable intermediate underflow. The final
+/// division can still round a mathematically valid boundary correlation just
+/// outside `[-1, 1]`; only after the exact admission proof is that numerical
+/// result projected back to the certified interval. This is not an empirical
+/// epsilon or tolerance.
 ///
-/// This routine validates the pairwise covariance bounds but does not claim a
-/// full positive-semidefinite proof. A caller that requires PSD admission must
-/// apply that model-specific invariant separately.
+/// Callers that need approximate-symmetry preprocessing must define and validate
+/// that policy explicitly before calling this kernel. This routine validates
+/// pairwise covariance bounds but does not claim a full positive-semidefinite
+/// proof; consumers requiring PSD admission retain that model-specific invariant.
 ///
 /// # Errors
 ///
@@ -239,8 +247,13 @@ pub fn standardize_covariance_matrix(
                 return Err(CovarianceStandardizationError::InvalidPairwiseCovariance);
             }
 
-            let standardized =
-                (upper / standard_deviations[row]) / standard_deviations[column];
+            let (first_sd, second_sd) = if standard_deviations[row] <= standard_deviations[column]
+            {
+                (standard_deviations[row], standard_deviations[column])
+            } else {
+                (standard_deviations[column], standard_deviations[row])
+            };
+            let standardized = (upper / first_sd) / second_sd;
             if !standardized.is_finite() {
                 return Err(CovarianceStandardizationError::NonFiniteResult);
             }

From 4f043d104b4de964c2c36aca3f539b3a8b08f81c Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:20:27 +0900
Subject: [PATCH 13/19] docs(core): document scale-order invariant

---
 docs/papers/covariance-standardization-owner-contract.md | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/docs/papers/covariance-standardization-owner-contract.md b/docs/papers/covariance-standardization-owner-contract.md
index 913892f59..fdeb4421c 100644
--- a/docs/papers/covariance-standardization-owner-contract.md
+++ b/docs/papers/covariance-standardization-owner-contract.md
@@ -30,13 +30,15 @@ The TEPP evidence is preserved rather than copied as a product name here: Driver
 
 Pairwise admission does not use an empirical epsilon. Each finite binary64 value is decomposed into its exact integer significand and power-of-two exponent, and the two products in `c² <= v_i * v_j` are compared exactly with `u128` significand products plus exponent alignment. Consequently, a genuinely invalid represented covariance fails closed even if later floating-point division would round its correlation back into range.
 
-The correlation itself is evaluated by sequentially dividing the covariance by each marginal standard deviation rather than first multiplying the two square roots. This avoids overflowing a representable result solely because `sqrt(v_i) * sqrt(v_j)` exceeds the finite binary64 range. Sequential division can still round a mathematically admissible boundary correlation one representable value outside `[-1, 1]`. After the exact admissibility proof, and only after that proof, the computed value is projected back to `[-1, 1]`. This projection is therefore a consequence of the exact represented-input bound rather than a tolerance for invalid covariance.
+The correlation itself is evaluated without first multiplying the two marginal standard deviations. After exact pairwise admission, the covariance is divided by the smaller standard deviation first and the larger standard deviation second. For an admissible pair the first quotient is bounded by the larger standard deviation, which avoids overflow and also prevents an avoidable zero caused by dividing a very small covariance by the larger scale first. This ordering makes the numerical result invariant to exchanging the two variables even when their variances differ by hundreds of orders of magnitude.
+
+The final division can still round a mathematically admissible boundary correlation one representable value outside `[-1, 1]`. After the exact admissibility proof, and only after that proof, the computed value is projected back to `[-1, 1]`. This projection is therefore a consequence of the exact represented-input bound rather than a tolerance for invalid covariance.
 
 The contract does not claim to prove full positive semidefiniteness. Model-specific PSD admission remains a separate invariant. Callers that wish to accept approximately symmetric observations must define and validate that preprocessing policy before calling this exact low-level kernel.
 
 ## Recovery and parity
 
-Unit and public-contract fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, exact-symmetry rejection, pairwise covariance beyond the exact correlation bound, zero covariance with a subnormal variance, and a finite exact-valid boundary case whose sequential binary64 divisions round to `next_up(1.0)` before the bound-certified projection.
+Unit and public-contract fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, exact-symmetry rejection, pairwise covariance beyond the exact correlation bound, zero covariance with a subnormal variance, a finite exact-valid boundary case whose sequential binary64 divisions round to `next_up(1.0)` before the bound-certified projection, and permutation invariance for an extreme-scale covariance whose correlation would underflow to zero if the larger marginal standard deviation were divided first.
 
 Before TEPP removes its local duplicate, the TEPP adapter must prove parity with its preserved `TIPREDVARstd` fixtures while continuing to enforce EventTime semantics outside this kernel.
 

From 759c9dbddcf78f5dc4a6877df0ffa457c290eae3 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:28:04 +0900
Subject: [PATCH 14/19] test(core): reject signed-zero covariance asymmetry

---
 .../tests/covariance_standardization_contract.rs       | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index c9386ab33..799662932 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -140,6 +140,16 @@ fn covariance_matrix_requires_exact_symmetry() {
     );
 }
 
+#[test]
+fn covariance_matrix_rejects_signed_zero_mirror_mismatch() {
+    let covariance = [1.0, 0.0, -0.0, 1.0];
+
+    assert_eq!(
+        standardize_covariance_matrix(&covariance, 2),
+        Err(CovarianceStandardizationError::NonSymmetricCovariance),
+    );
+}
+
 #[test]
 fn covariance_matrix_never_clamps_an_out_of_range_correlation() {
     let above_one = 1.0 + 64.0 * f64::EPSILON;

From 73a930e76b084ea7eae08465bb0e72646fbdac6a Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 05:30:35 +0900
Subject: [PATCH 15/19] fix(core): preserve exact binary64 covariance symmetry

---
 crates/mlsirm-core/src/covariance_standardization.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
index 019166886..a5832bcc4 100644
--- a/crates/mlsirm-core/src/covariance_standardization.rs
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -237,7 +237,7 @@ pub fn standardize_covariance_matrix(
         for column in (row + 1)..dimension {
             let upper = covariance[row * dimension + column];
             let lower = covariance[column * dimension + row];
-            if upper != lower {
+            if upper.to_bits() != lower.to_bits() {
                 return Err(CovarianceStandardizationError::NonSymmetricCovariance);
             }
 

From eb4664b456248bf06e47466a515f3c819e45a4cb Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 06:03:27 +0900
Subject: [PATCH 16/19] test(core): require exact scalar covariance
 standardization

---
 .../tests/covariance_standardization_contract.rs       | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/crates/mlsirm-core/tests/covariance_standardization_contract.rs b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
index 799662932..e9e39d9bb 100644
--- a/crates/mlsirm-core/tests/covariance_standardization_contract.rs
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -18,19 +18,19 @@ fn public_covariance_standardization_contract_is_versioned_and_scale_invariant()
     );
 
     for variance in [
+        f64::from_bits(1),
         f64::MIN_POSITIVE,
         1.0e-200,
         0.25,
         1.0,
+        3.0,
         6.4,
         1.0e200,
         f64::MAX,
     ] {
-        assert_close(
-            standardize_variance(variance).expect("strictly positive finite variance"),
-            1.0,
-            8.0e-15,
-        );
+        let standardized =
+            standardize_variance(variance).expect("strictly positive finite variance");
+        assert_eq!(standardized.to_bits(), 1.0_f64.to_bits());
     }
 }
 

From d2c0a81505177810dc1f4ea9dc73eeeb93e3d6b8 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 06:05:16 +0900
Subject: [PATCH 17/19] fix(core): preserve exact scalar standardization
 identity

---
 .../src/covariance_standardization.rs         | 32 ++++++++-----------
 1 file changed, 14 insertions(+), 18 deletions(-)

diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
index a5832bcc4..f84dba2a2 100644
--- a/crates/mlsirm-core/src/covariance_standardization.rs
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -148,16 +148,18 @@ fn pairwise_covariance_is_admissible(
 
 /// Standardize one finite, strictly positive variance against itself.
 ///
-/// The arithmetic is evaluated rather than replaced with a hard-coded `1.0` so
-/// this scalar reference exercises the same normalization contract consumed by
-/// matrix standardization and downstream parity tests.
+/// A one-dimensional covariance standardizes exactly to unit correlation. After
+/// validating `variance > 0`, the algebraically equivalent ratio `variance /
+/// variance` keeps an actual binary64 arithmetic path while avoiding the
+/// avoidable rounding introduced by separately evaluating inverse square-root
+/// factors. Every admitted finite positive binary64 value therefore returns the
+/// exact binary64 representation of `1.0`.
 ///
 /// # Errors
 ///
 /// Returns [`CovarianceStandardizationError::NonFiniteInput`] for NaN or
-/// infinity, [`CovarianceStandardizationError::NonPositiveVariance`] for zero
-/// or a negative value, and [`CovarianceStandardizationError::NonFiniteResult`]
-/// if the arithmetic cannot produce a finite value.
+/// infinity and [`CovarianceStandardizationError::NonPositiveVariance`] for zero
+/// or a negative value.
 pub fn standardize_variance(
     variance: f64,
 ) -> Result<f64, CovarianceStandardizationError> {
@@ -168,12 +170,7 @@ pub fn standardize_variance(
         return Err(CovarianceStandardizationError::NonPositiveVariance);
     }
 
-    let inverse_sd = 1.0 / variance.sqrt();
-    let standardized = (variance * inverse_sd) * inverse_sd;
-    if !standardized.is_finite() {
-        return Err(CovarianceStandardizationError::NonFiniteResult);
-    }
-    Ok(standardized)
+    Ok(variance / variance)
 }
 
 /// Convert a finite symmetric covariance matrix to a correlation matrix.
@@ -291,21 +288,20 @@ mod tests {
     }
 
     #[test]
-    fn scalar_reference_recovers_one_across_positive_scales() {
+    fn scalar_reference_recovers_exact_one_across_positive_scales() {
         for variance in [
+            f64::from_bits(1),
             f64::MIN_POSITIVE,
             1.0e-200,
             0.25,
             1.0,
+            3.0,
             6.4,
             1.0e200,
             f64::MAX,
         ] {
-            assert_close(
-                standardize_variance(variance).expect("positive variance"),
-                1.0,
-                8.0e-15,
-            );
+            let standardized = standardize_variance(variance).expect("positive variance");
+            assert_eq!(standardized.to_bits(), 1.0_f64.to_bits());
         }
     }
 

From e1847c07fd7ef8331dcebd0dd588b1381cc1231d Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Wed, 2 Sep 2026 06:06:21 +0900
Subject: [PATCH 18/19] docs(core): record exact scalar covariance identity

---
 docs/papers/covariance-standardization-owner-contract.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/papers/covariance-standardization-owner-contract.md b/docs/papers/covariance-standardization-owner-contract.md
index fdeb4421c..2df806493 100644
--- a/docs/papers/covariance-standardization-owner-contract.md
+++ b/docs/papers/covariance-standardization-owner-contract.md
@@ -8,7 +8,7 @@
 R = D^{-1/2}\Sigma D^{-1/2},
 \]
 
-where `D` is the diagonal of strictly positive marginal variances. The scalar self-standardization is `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`.
+where `D` is the diagonal of strictly positive marginal variances. The scalar self-standardization is `(1 / sqrt(v)) * v * (1 / sqrt(v)) = 1` for finite `v > 0`. After validating the scalar variance, the implementation evaluates the algebraically equivalent ratio `v / v`, so every admitted finite positive binary64 value—including the smallest positive subnormal and `f64::MAX`—returns the exact binary64 representation of `1.0` rather than accumulating avoidable square-root rounding.
 
 This contract is intentionally domain-neutral. It does not decide event time, valid time, knowledge cutoff, state evolution, temporal identification, or product-specific model activation. Those policies belong to consuming bounded contexts such as TEPP Longitudinal Modeling.
 
@@ -38,7 +38,7 @@ The contract does not claim to prove full positive semidefiniteness. Model-speci
 
 ## Recovery and parity
 
-Unit and public-contract fixtures cover positive scalar magnitudes from `f64::MIN_POSITIVE` through `f64::MAX`, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, exact-symmetry rejection, pairwise covariance beyond the exact correlation bound, zero covariance with a subnormal variance, a finite exact-valid boundary case whose sequential binary64 divisions round to `next_up(1.0)` before the bound-certified projection, and permutation invariance for an extreme-scale covariance whose correlation would underflow to zero if the larger marginal standard deviation were divided first.
+Unit and public-contract fixtures cover the smallest positive subnormal scalar variance, `f64::MIN_POSITIVE`, ordinary values including `3.0`, very large finite values through `f64::MAX`, exact binary64 unit recovery, malformed scalar inputs, a known 2x2 covariance/correlation pair, multiplicative scale invariance, shape errors, non-finite cells, non-positive diagonals, exact-symmetry rejection, pairwise covariance beyond the exact correlation bound, zero covariance with a subnormal variance, a finite exact-valid boundary case whose sequential binary64 divisions round to `next_up(1.0)` before the bound-certified projection, and permutation invariance for an extreme-scale covariance whose correlation would underflow to zero if the larger marginal standard deviation were divided first.
 
 Before TEPP removes its local duplicate, the TEPP adapter must prove parity with its preserved `TIPREDVARstd` fixtures while continuing to enforce EventTime semantics outside this kernel.
 

From 28b0305595107fd0ba21d7b27c1ac5db68ae8bf1 Mon Sep 17 00:00:00 2001
From: Seongho Bae <me@seonghobae.me>
Date: Mon, 7 Sep 2026 11:36:30 +0900
Subject: [PATCH 19/19] docs(changelog): record covariance standardization

---
 docs/changelog.d/1722-covariance-standardization.md | 6 ++++++
 1 file changed, 6 insertions(+)
 create mode 100644 docs/changelog.d/1722-covariance-standardization.md

diff --git a/docs/changelog.d/1722-covariance-standardization.md b/docs/changelog.d/1722-covariance-standardization.md
new file mode 100644
index 000000000..2b3ad6be7
--- /dev/null
+++ b/docs/changelog.d/1722-covariance-standardization.md
@@ -0,0 +1,6 @@
+# Covariance standardization owner contract
+
+## Added
+
+- Add the Rust-owned `fast_mlsirm.covariance_standardization@1.0.0` public
+  contract for exact, overflow-safe covariance-to-correlation standardization.