diff --git a/crates/mlsirm-core/src/covariance_standardization.rs b/crates/mlsirm-core/src/covariance_standardization.rs
new file mode 100644
index 000000000..f84dba2a2
--- /dev/null
+++ b/crates/mlsirm-core/src/covariance_standardization.rs
@@ -0,0 +1,376 @@
+//! 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`.
+//! 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
+//! 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::cmp::Ordering;
+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 are not exactly equal in binary64.
+    NonSymmetricCovariance,
+    /// A represented covariance pair violates `c² <= v_i * v_j` exactly.
+    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 {}
+
+/// 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.
+///
+/// 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 and [`CovarianceStandardizationError::NonPositiveVariance`] for zero
+/// or a negative 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);
+    }
+
+    Ok(variance / 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. Mirrored off-diagonal cells
+/// 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, 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.
+///
+/// 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
+///
+/// Returns a typed error for invalid shape, non-finite input, non-positive
+/// diagonal variance, asymmetric mirrored cells, 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];
+            if upper.to_bits() != lower.to_bits() {
+                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 (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);
+            }
+            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::{
+        scaled_integer_le, 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 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_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,
+        ] {
+            let standardized = standardize_variance(variance).expect("positive variance");
+            assert_eq!(standardized.to_bits(), 1.0_f64.to_bits());
+        }
+    }
+
+    #[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_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!(
+            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)
+        );
+    }
+}
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;
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..e9e39d9bb
--- /dev/null
+++ b/crates/mlsirm-core/tests/covariance_standardization_contract.rs
@@ -0,0 +1,162 @@
+use mlsirm_core::covariance_standardization::{
+    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 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::from_bits(1),
+        f64::MIN_POSITIVE,
+        1.0e-200,
+        0.25,
+        1.0,
+        3.0,
+        6.4,
+        1.0e200,
+        f64::MAX,
+    ] {
+        let standardized =
+            standardize_variance(variance).expect("strictly positive finite variance");
+        assert_eq!(standardized.to_bits(), 1.0_f64.to_bits());
+    }
+}
+
+#[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]
+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 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;
+    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_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;
+    let covariance = [1.0, above_one, above_one, 1.0];
+
+    assert_eq!(
+        standardize_covariance_matrix(&covariance, 2),
+        Err(CovarianceStandardizationError::InvalidPairwiseCovariance),
+    );
+}
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.
diff --git a/docs/papers/covariance-standardization-owner-contract.md b/docs/papers/covariance-standardization-owner-contract.md
new file mode 100644
index 000000000..2df806493
--- /dev/null
+++ b/docs/papers/covariance-standardization-owner-contract.md
@@ -0,0 +1,47 @@
+# 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`. 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.
+
+## 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;
+- 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`.
+
+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 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 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.
+
+## 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