test: preserve Dev verification contracts through Refactor APIs

docs/cpp-object-oriented-modular-refactoring/merge-manifest.json

@@ -32530,5 +32530,119 @@
         "cause": "Temporary Python cmd argument quoting passed full command as quoted executable; corrected temporary runner only; first log retained as msbuild.log."
       }
     }
+  },
+  "t3_test_migration": {
+    "base_commit": "26009169bf51a05ce2516f73e75487f593d1bd25",
+    "scope": "검사와 읽기 전용 comparator helper만 변경; 제품·reference 수정 및 최종 삭제 없음",
+    "migration_table": "docs/cpp-object-oriented-modular-refactoring/test-migration.tsv",
+    "table_summary": {
+      "common_files": 130,
+      "byte_identical": 72,
+      "byte_different": 58,
+      "dev_only": [
+        "tests/test_cpp_policy_contract.py",
+        "tests/unit/math/high_low_arithmetic_original.h"
+      ],
+      "refactor_only": [
+        "tests/test_compare_hex8.py",
+        "tests/test_compare_hht_dynamic.py",
+        "tests/unit/elements/element_capabilities_test.cpp"
+      ],
+      "rows": 2876,
+      "missing_identities": [],
+      "deletions": []
+    },
+    "raw_evidence": ".build/merge-verification/t3-20261018-01",
+    "builds": {
+      "final-cmake-03": {
+        "command": "cmake --build .build/cmake-debug --target fesa_tests --parallel 2",
+        "exit": 0,
+        "seconds": 3.9731388092041016,
+        "log": ".build\\merge-verification\\t3-20261018-01\\final-cmake-03.log"
+      },
+      "final-native": {
+        "command": "msbuild FESA.sln /m:2 /p:Configuration=Debug /p:Platform=x64 /t:fesa_unit_tests,fesa_integration_tests,fesa_reference_tests",
+        "exit": 0,
+        "seconds": 4.966400384902954,
+        "log": ".build\\merge-verification\\t3-20261018-01\\final-native.log"
+      },
+      "ctest-discovery": {
+        "command": "ctest --test-dir .build/cmake-debug --show-only=json-v1",
+        "exit": 0,
+        "seconds": 2.725306987762451,
+        "log": ".build\\merge-verification\\t3-20261018-01\\ctest-discovery.log"
+      }
+    },
+    "registration": {
+      "fesa_unit_tests": {
+        "native_count": 90,
+        "cmake_count": 90,
+        "native_only": [],
+        "cmake_only": []
+      },
+      "fesa_integration_tests": {
+        "native_count": 11,
+        "cmake_count": 11,
+        "native_only": [],
+        "cmake_only": []
+      },
+      "fesa_reference_tests": {
+        "native_count": 17,
+        "cmake_count": 17,
+        "native_only": [],
+        "cmake_only": []
+      }
+    },
+    "discovery_counts": {
+      "fesa_unit_tests": 849,
+      "fesa_integration_tests": 72,
+      "fesa_reference_tests": 43
+    },
+    "ctest_discovered": 964,
+    "focused_results": {
+      "final-cpp-focused": {
+        "tests": 176,
+        "failures": 0,
+        "errors": 0,
+        "skipped": 0
+      },
+      "final-native-focused": {
+        "tests": 15,
+        "failures": 0,
+        "errors": 0,
+        "skipped": 0
+      },
+      "final-reference-focused": {
+        "tests": 28,
+        "failures": 0,
+        "errors": 0,
+        "skipped": 0
+      },
+      "python-focused": {
+        "tests": 145,
+        "failures": 0,
+        "errors": 0,
+        "skipped": 0
+      }
+    },
+    "source_snapshot": ".build/merge-verification/t3-20261018-01/final-source-hashes.json",
+    "product_final": {
+      "count": 274,
+      "mismatches": []
+    },
+    "reference_final": {
+      "dev_tracked_files": 44,
+      "byte_mismatches": [],
+      "hex8_frozen_hashes_matched": 14
+    },
+    "pending": "canonical shellNL-meshfine 세 변형 자료는 승인 대기 상태로 보존하며 복원 의존 검증은 미실행. 전체 CTest/Python/reference 검증은 T5 범위.",
+    "structural_migrations": [
+      "AddPair/MultiplyPair 및 ScalarPacketAdd bool 경계는 공통 checked 산술로 이관하고 두 기존 consumer의 독립 frozen oracle를 유지한다.",
+      "Cholesky inverse out-parameter·reason enum은 DenseFactor lower action과 typed Result 실패로 이관하며 이전 owned candidate 보존을 검사한다.",
+      "rotation inverse matrix API는 모든 basis solve로 이관하고 positive unit null의 동등한 analytic invariant를 유지한다."
+    ],
+    "removed_tests": [],
+    "policy_failure_reproduction": "native ExceptionHandling을 메모리에서 false로 바꾸면 현재 policy 검사가 AssertionError로 거부하고 원문 설정에서는 통과한다. 파일은 수정하지 않는다.",
+    "policy_failure_evidence": ".build/merge-verification/t3-20261018-01/policy-mutation.log"
   }
 }

scripts/AGENTS.md

@@ -0,0 +1,11 @@
+# 비교 보조 도구
+
+`verification/compare_hex8.py`와 `compare_hht_dynamic.py`는 입력·reference·결과를
+읽고 요청한 JSON report만 기록한다. Solver를 실행하지 않는다.
+
+- 실행 전에 argv, 의존성, report 경로와 덮어쓰기 방지 검사를 확인한다.
+- 원본 입력·CSV·TXT·결과 bytes와 frozen hash를 변경하지 않는다. 자료 부재나
+  hash 충돌은 실패 또는 미실행으로 보고하며 원본 보정으로 해결하지 않는다.
+- identity·schema·finite 사전 gate와 승인된 family별 판정 영향을 유지한다.
+  Synthetic HDF5 fixture의 comparator 통과를 실제 solver의 물리 검증으로 보고하지 않는다.
+- 실행 근거는 고유한 작업 directory에 남기고 기존 report와 로그를 덮어쓰지 않는다.

scripts/verification/compare_hex8.py

@@ -0,0 +1,807 @@
+"""Compare original Hex8 CSV/TXT results with Solid HDF5; never launch a solver.
+
+This entry point is intentionally restricted to reference/hexa's fixed original
+models. Identity/schema checks precede finite checks and the common tolerance.
+Rounded HDF5 observables are compared; internal high+low solver gates are not
+recomputed from them. Missing original modal Value data remain NOT RUN.
+"""
+
+import argparse
+import csv
+from decimal import Decimal, InvalidOperation
+import hashlib
+import json
+import math
+import os
+from pathlib import Path
+import re
+from typing import Any
+
+import h5py
+import numpy as np
+
+
+ROOT = Path(__file__).resolve().parents[2]
+REFERENCE_DIRECTORY = ROOT / "reference/hexa"
+STEP = "Step-1"
+INSTANCE = "PART-1-1"
+NODE_COUNT = 442
+ELEMENT_COUNT = 300
+POLICY = {"relative_tolerance": 0.05, "relative_rms_tolerance": 0.01,
+          "near_zero_ratio": 0.01, "zero_scale": "exact-zero"}
+CASES = {
+    "hexa_linearstatic": "linear-static",
+    "hexa_nonlinearstatic": "nonlinear-static-full-newton",
+    "hexa_dynlinear": "hht-linear-dynamics",
+    "hexa_dynnonlinear": "hht-nonlinear-dynamics",
+    "hexa_freq": "linear-frequency",
+}
+# Working-tree bytes (Windows CRLF), not Git's normalized LF blobs.
+SOURCE_SHA256 = {
+    "hexa_dynlinear.inp": "f8db81b8ea56b689ab0ad892e6a84b7c6a0162c9d3ff95328ea99722f3de078e",
+    "hexa_dynlinear_U.csv": "623c20489e8cc18d23fa1728563ce5a26059f321861514ef8569e3d204ef6cb5",
+    "hexa_dynnonlinear.inp": "1685f86f844e8020bf46e1747b949699ac9aff0d7570bd21ec95109399365ff9",
+    "hexa_dynnonlinear_U.csv": "f5832a24577f54838bb49dec86a0dc9668cf43e78953f718a9d3635e0f574226",
+    "hexa_freq.inp": "42cbcec2a8b738f8b6d2c1e13114b8840dccdde6fc3e6f8a43a3272e85d82347",
+    "hexa_freq.txt": "3d2dacbee13cd2cc8d3f7c1a38210611a72de41d5e31d169e636bf40e00d4370",
+    "hexa_linearstatic.inp": "67183909c5dc962a4d4b65ffb81f7bd65b72e89a5d015c8369a8e0226347a647",
+    "hexa_linearstatic_RF.csv": "3e04e71ec075ffef098cf4870d77f9c400557fccb375369da8bf67cf67c1e36b",
+    "hexa_linearstatic_S.csv": "8840a502b9d28a91e5685349301a56fed37d2ac62327b71304716b714db577e8",
+    "hexa_linearstatic_U.csv": "1ff49aea3ec601cfd796d6a2fb6a4f91be3fc3de9426cdf3744d981216d4bd5f",
+    "hexa_nonlinearstatic.inp": "7858254025b96c6bfe3fb6498605ecf9c473425c2bc3593f5fbdc6902fd7c6cb",
+    "hexa_nonlinearstatic_RF.csv": "3ba67a4fb75b2ce3edee13d35b2ea66ac7d52fc16e612c086ed1122a834fee7b",
+    "hexa_nonlinearstatic_S.csv": "a8866a730f0d65b1e3921d504d255b72158b4a2b9ee30fc691712d94982b07ba",
+    "hexa_nonlinearstatic_U.csv": "bd6c2fe2b2bb2ffb0bc2ea52f67513c3967a396b5707bbebc9533fde7be25ac6",
+}
+BASE_COLUMNS = ["ODB Name", "Step", "Frame", "Part Instance Name"]
+TAIL_COLUMNS = ["X", "Y", "Z", "Section Name", "Material Name", "Section Point"]
+COMPONENTS = {"U": ["U1", "U2", "U3"], "RF": ["RF1", "RF2", "RF3"],
+              "S": ["S11", "S22", "S33", "S12", "S13", "S23"]}
+STRESS_DERIVED = ["Mises", "Max. Principal", "Max. Principal (Abs)",
+                  "Mid. Principal", "Min. Principal", "Tresca", "Pressure", "Third Invariant"]
+NATURAL_NODES = np.array([[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1],
+                          [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]])
+NATURAL_POINTS = np.array([[-1, -1, -1], [1, -1, -1], [-1, 1, -1], [1, 1, -1],
+                           [-1, -1, 1], [1, -1, 1], [-1, 1, 1], [1, 1, 1]]) / math.sqrt(3.0)
+FRAME_PATTERN = re.compile(
+    r"Increment\s+(\d+):\s*Step Time\s*=\s*([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)"
+)
+
+
+class ComparisonError(ValueError):
+    def __init__(self, code: str, message: str):
+        super().__init__(message)
+        self.code = code
+
+
+def require(condition: bool, code: str, message: str):
+    if not condition:
+        raise ComparisonError(code, message)
+
+
+def canonical(value: str) -> str:
+    """Abaqus ASCII case folding only; distinct instances never coalesce."""
+    return value.translate(str.maketrans("abcdefghijklmnopqrstuvwxyz", "ABCDEFGHIJKLMNOPQRSTUVWXYZ"))
+
+
+def label(value: str) -> int:
+    require(re.fullmatch(r"[0-9]+", value) is not None, "invalid-source-label", value)
+    number = int(value)
+    require(0 < number <= 2**63 - 1, "invalid-source-label", value)
+    return number
+
+
+def lexical_absolute(path: Path) -> str:
+    """Match C++ absolute.lexically_normal without resolving physical aliases."""
+    return Path(os.path.abspath(path)).as_posix()
+
+
+def sha256(path: Path) -> str:
+    with path.open("rb") as stream:
+        return hashlib.file_digest(stream, "sha256").hexdigest()
+
+
+def fnv1a64(raw: bytes) -> str:
+    value = 14695981039346656037
+    for byte in raw:
+        value = ((value ^ byte) * 1099511628211) & (2**64 - 1)
+    return f"fnv1a64:{value:016x}"
+
+
+def input_inventory(path: Path, case: str) -> dict[str, Any]:
+    """Read only the frozen single-part identity-instance original subset."""
+    require(case in CASES, "unknown-source-case", case)
+    require(sha256(path) == SOURCE_SHA256[f"{case}.inp"], "source-input-hash-mismatch", str(path))
+    nodes, elements, sets, boundaries = {}, {}, {}, []
+    keyword, options, in_assembly = "", {}, False
+    model = {"case": case, "nodes": nodes, "elements": elements, "node_sets": sets,
+             "section_line": None, "source_path": lexical_absolute(path),
+             "content_identity": fnv1a64(path.read_bytes())}
+    for line_number, raw in enumerate(path.read_text(encoding="utf-8-sig").splitlines(), 1):
+        line = raw.strip()
+        if not line or line.startswith("**"):
+            continue
+        if line.startswith("*"):
+            tokens = [token.strip() for token in line[1:].split(",")]
+            keyword = canonical(tokens[0])
+            options = {canonical(token.split("=", 1)[0]): token.split("=", 1)[1].strip()
+                       for token in tokens[1:] if "=" in token}
+            if keyword == "ASSEMBLY":
+                in_assembly = True
+            elif keyword == "END ASSEMBLY":
+                in_assembly = False
+            elif keyword == "INSTANCE":
+                require(canonical(options["NAME"]) == INSTANCE and len(tokens) == 3,
+                        "unsupported-source-instance", line)
+            elif keyword == "SOLID SECTION":
+                model["section_line"] = line_number
+                model["section_elset"] = options["ELSET"]
+                model["material"] = options["MATERIAL"]
+            elif keyword == "STEP":
+                require(options["NAME"] == STEP, "wrong-source-step", line)
+            continue
+        tokens = [token.strip() for token in line.split(",") if token.strip()]
+        if keyword == "NODE":
+            node = label(tokens[0])
+            require(node not in nodes and len(tokens) == 4, "source-node-inventory", line)
+            nodes[node] = tuple(float(value) for value in tokens[1:])
+        elif keyword == "ELEMENT":
+            element = label(tokens[0])
+            require(element not in elements and len(tokens) == 9 and options["TYPE"] == "C3D8",
+                    "source-element-inventory", line)
+            elements[element] = tuple(label(value) for value in tokens[1:])
+        elif keyword == "NSET" and in_assembly:
+            name = canonical(options["NSET"])
+            require(canonical(options["INSTANCE"]) == INSTANCE, "wrong-source-instance", line)
+            sets.setdefault(name, []).extend(label(value) for value in tokens)
+        elif keyword == "BOUNDARY":
+            first, last = int(tokens[1]), int(tokens[2])
+            require(len(tokens) == 3 or float(tokens[3]) == 0.0, "unsupported-source-boundary", line)
+            boundaries.append((canonical(tokens[0]), first, last))
+        elif keyword == "DENSITY":
+            model["density"] = float(tokens[0])
+        elif keyword == "ELASTIC":
+            model["elastic"] = tuple(float(value) for value in tokens)
+    require(set(nodes) == set(range(1, NODE_COUNT + 1))
+            and set(elements) == set(range(1, ELEMENT_COUNT + 1)),
+            "source-model-inventory", "Expected original 442 nodes and 300 elements.")
+    constrained = set()
+    for name, first, last in boundaries:
+        require(name in sets and 1 <= first <= last <= 6, "source-boundary-identity", name)
+        constrained.update((node, component) for node in sets[name]
+                           for component in range(first, min(last, 3) + 1))
+    model["constrained"] = constrained
+    return model
+
+
+def point_coordinates(model, element: int, point: int):
+    xi = NATURAL_POINTS[point - 1]
+    shape = np.prod(1.0 + NATURAL_NODES * xi, axis=1) / 8.0
+    coordinates = np.array([model["nodes"][node] for node in model["elements"][element]])
+    terms = shape[:, None] * coordinates
+    return np.sum(terms, axis=0), np.sum(np.abs(terms), axis=0)
+
+
+def read_reference(path: Path, case: str, quantity: str, model):
+    """Validate all row identities before converting any physical values."""
+    require(quantity in COMPONENTS, "unknown-quantity", quantity)
+    dynamic = case.startswith("hexa_dyn")
+    require(not dynamic or quantity == "U", "wrong-reference-quantity", quantity)
+    extra = (["Element Label", "IntPt"] if quantity == "S" else ["Node Label"])
+    suffix = ([f"S-{name}" for name in STRESS_DERIVED] if quantity == "S"
+              else [f"{quantity}-Magnitude"])
+    headers = BASE_COLUMNS + extra + TAIL_COLUMNS + suffix + [f"{quantity}-{c}" for c in COMPONENTS[quantity]]
+    records = {}
+    with path.open(newline="", encoding="utf-8-sig") as stream:
+        reader = csv.reader(stream)
+        fields = [field.strip() for field in next(reader, [])]
+        require(len(fields) == len(set(fields)), "duplicate-reference-header", str(path))
+        require(set(fields) == set(headers), "wrong-reference-columns", str(fields))
+        for row_number, values in enumerate(reader, 2):
+            require(len(values) == len(fields), "wrong-reference-row-width", str(row_number))
+            row = dict(zip(fields, [value.strip() for value in values]))
+            require(row["ODB Name"].replace("\\", "/").rsplit("/", 1)[-1] == f"{case}.odb",
+                    "wrong-reference-case", str(row_number))
+            require(row["Step"] == STEP, "wrong-reference-step", str(row_number))
+            require(canonical(row["Part Instance Name"]) == INSTANCE,
+                    "wrong-reference-instance", str(row_number))
+            match = FRAME_PATTERN.fullmatch(row["Frame"])
+            require(match is not None, "invalid-reference-frame", str(row_number))
+            increment = int(match[1])
+            time = Decimal(match[2])
+            require(time.is_finite() and increment in (range(101) if dynamic else (10,))
+                    and time == (Decimal(increment) / 10 if dynamic else Decimal(1)),
+                    "wrong-reference-frame-time", str(row_number))
+            entity = label(row["Element Label" if quantity == "S" else "Node Label"])
+            point = label(row["IntPt"]) if quantity == "S" else 0
+            require(entity in model["elements" if quantity == "S" else "nodes"],
+                    "extra-reference-entity", str(row_number))
+            require(quantity != "S" or 1 <= point <= 8, "wrong-reference-point", str(row_number))
+            require(row["Section Point"] == ""
+                    and canonical(row["Material Name"]) == ("MATERIAL-1" if quantity == "S" else "NA")
+                    and canonical(row["Section Name"]) == ("ASSEMBLY_PART-1-1_SET-1" if quantity == "S" else "NA"),
+                    "wrong-reference-section", str(row_number))
+            key = (increment, entity, point)
+            require(key not in records, "duplicate-reference-identity", str(key))
+            records[key] = row
+    entities = model["elements" if quantity == "S" else "nodes"]
+    expected = {(n, entity, p) for n in (range(101) if dynamic else (10,))
+                for entity in entities for p in (range(1, 9) if quantity == "S" else (0,))}
+    require(records.keys() == expected, "reference-inventory", "Missing reference identities.")
+    return records
+
+
+def reference_values(records, quantity, model):
+    result = {}
+    for key in sorted(records):
+        row = records[key]
+        values = [float(row[f"{quantity}-{component}"]) for component in COMPONENTS[quantity]]
+        coordinates = np.array([float(row[axis]) for axis in ("X", "Y", "Z")])
+        # Every exported numeric column is checked, including unused derived values.
+        derived = STRESS_DERIVED if quantity == "S" else ["Magnitude"]
+        other = [float(row[f"{quantity}-{name}"]) for name in derived]
+        require(all(math.isfinite(value) for value in values + other) and np.isfinite(coordinates).all(),
+                "nonfinite-reference-value", str(key))
+        if quantity == "S":
+            expected, scale = point_coordinates(model, key[1], key[2])
+        else:
+            expected = np.array(model["nodes"][key[1]])
+            scale = np.abs(expected)
+        budget = 16.0 * np.finfo(np.float32).eps * scale
+        require((np.abs(coordinates - expected) <= budget).all(),
+                "reference-coordinate-convention", str(key))
+        result[key] = values
+    return result
+
+
+def evaluate_family(identity, result, reference):
+    require(len(result) == len(reference) == len(identity["rows"]) and len(result) > 0,
+            "scalar-inventory", identity["quantity"])
+    require(all(math.isfinite(value) for value in (*result, *reference)),
+            "nonfinite-comparison-value", "Comparison values must be finite.")
+    scale = max(abs(value) for value in reference)
+    band = POLICY["near_zero_ratio"] * scale
+    rows, error_norm = [], 0.0
+    for key, actual, expected in zip(identity["rows"], result, reference):
+        error = abs(actual - expected)
+        require(math.isfinite(error), "nonfinite-comparison-value", str(key))
+        error_norm = math.hypot(error_norm, error)
+        require(math.isfinite(error_norm), "nonfinite-comparison-value", str(key))
+        relative = None
+        if scale == 0.0:
+            branch, threshold, passed = "zero-scale-exact", 0.0, error == 0.0
+        elif abs(expected) <= band:
+            branch, threshold, passed = "near-zero", band, error <= band
+        else:
+            branch, threshold = "relative", POLICY["relative_tolerance"] * abs(expected)
+            relative = error / abs(expected)
+            require(math.isfinite(relative), "nonfinite-comparison-value", str(key))
+            passed = relative <= POLICY["relative_tolerance"]
+        rows.append({**key, "fesa_value": actual, "reference_value": expected,
+                     "absolute_error": error, "relative_error": relative, "branch": branch,
+                     "threshold": threshold, "passed": passed})
+    rows_passed = all(row["passed"] for row in rows)
+    absolute_rms = error_norm / math.sqrt(len(rows))
+    relative_rms = absolute_rms / scale if scale else None
+    require(relative_rms is None or math.isfinite(relative_rms), "nonfinite-comparison-value", "RMS")
+    rms_passed = rows_passed if not scale else relative_rms <= POLICY["relative_rms_tolerance"]
+    return {**{k: v for k, v in identity.items() if k != "rows"}, "impact": "blocking",
+            "row_count": len(rows), "reference_scale": scale, "near_zero_band": band,
+            "maximum_absolute_error": max(row["absolute_error"] for row in rows),
+            "absolute_rms": absolute_rms, "relative_rms": relative_rms,
+            "policy": dict(POLICY), "rows": rows, "rows_passed": rows_passed,
+            "rms_passed": rms_passed, "passed": rows_passed and rms_passed,
+            "diagnostic_code": "zero-reference-scale-nonzero-error" if not scale and not rows_passed else None}
+
+
+def read_modal_reference(path: Path):
+    """Read the supplied hexa_freq.txt Value=lambda, Freq=cycles/time format.
+
+    This artifact has no step field: its fixed case is bound to the sole Step-1
+    of the hash-checked hexa_freq.inp. Freq is finite auxiliary data; comparison
+    uses the printed Value directly, without frequency conversion or rounding.
+    """
+    mapping = {}
+    for line in path.read_text(encoding="utf-8-sig").splitlines():
+        if not line.strip():
+            continue
+        match = re.fullmatch(
+            r"\s*Mode\s+([0-9]+):\s+Value\s*=\s*(\S+)\s+Freq\s*=\s*(\S+)\s+\(cycles/time\)\s*", line)
+        require(match is not None, "wrong-modal-reference-format", line)
+        mode = int(match[1])
+        require(mode in range(1, 11), "wrong-modal-mode", str(mode))
+        require(mode not in mapping, "duplicate-modal-mode", str(mode))
+        mapping[mode] = (match[2], match[3])
+    require(set(mapping) == set(range(1, 11)), "modal-mode-inventory", "Lowest ten modes required.")
+    rows = []
+    for mode in range(1, 11):
+        value, frequency = (float(token) for token in mapping[mode])
+        require(math.isfinite(value) and math.isfinite(frequency), "nonfinite-modal-reference", str(mode))
+        require(value >= 0 and frequency >= 0, "invalid-modal-reference-value", str(mode))
+        rows.append({"step": STEP, "mode": mode, "quantity": "eigenvalue",
+                     "unit_dimensions": "1/time^2", "value": value})
+    require(all(a["value"] <= b["value"] for a, b in zip(rows, rows[1:])),
+            "wrong-modal-order", str(path))
+    return rows
+
+
+def compare_modal(reference, result, quantity: str, unit: str):
+    """Compare explicitly typed modal rows using the common tolerance.
+
+    Each row explicitly supplies step, mode, quantity, unit_dimensions and value.
+    Callers must obtain those semantics from a genuine supplied artifact first.
+    Eigenvalues and frequencies are never silently converted into one another.
+    """
+    units = {"eigenvalue": "1/time^2", "frequency": "1/time"}
+    require(quantity in units and unit == units[quantity], "wrong-modal-unit", unit)
+    inventories = []
+    for rows in (reference, result):
+        mapping = {}
+        for row in rows:
+            require(row["quantity"] == quantity, "wrong-modal-quantity", str(row["mode"]))
+            require(row["unit_dimensions"] == unit, "wrong-modal-unit", str(row["mode"]))
+            require(row["step"] == STEP, "wrong-modal-step", str(row["mode"]))
+            mode = row["mode"]
+            require(type(mode) is int and mode in range(1, 11), "wrong-modal-mode", str(mode))
+            require(mode not in mapping, "duplicate-modal-mode", str(mode))
+            mapping[mode] = row["value"]
+        require(set(mapping) == set(range(1, 11)), "modal-mode-inventory", "Lowest ten modes required.")
+        inventories.append(mapping)
+    families = [evaluate_family(
+        {"identity": f"hexa_freq:{STEP}:mode-{mode}:{quantity}",
+         "case": "hexa_freq", "step": STEP, "mode": mode, "quantity": quantity,
+         "unit_dimensions": unit, "rows": [{"mode": mode, "component": quantity}]},
+        [inventories[1][mode]], [inventories[0][mode]]) for mode in range(1, 11)]
+    return {"families": families, "passed": all(family["passed"] for family in families)}
+
+
+def text(value) -> str:
+    return value.decode("utf-8", errors="strict") if isinstance(value, bytes) else str(value)
+
+
+def string_attribute(obj, name, expected=None):
+    require(name in obj.attrs, "missing-result-attribute", f"{obj.name}:{name}")
+    attr = obj.attrs.get_id(name)
+    string = h5py.check_string_dtype(attr.dtype)
+    require(attr.shape == () and string is not None and string.length is None
+            and string.encoding == "utf-8", "wrong-result-attribute-type", f"{obj.name}:{name}")
+    value = text(obj.attrs[name])
+    require(bool(value) and (expected is None or value == expected),
+            "wrong-result-attribute", f"{obj.name}:{name}={value}")
+    return value
+
+
+def numeric_attribute(obj, name, dtype, expected):
+    require(name in obj.attrs, "missing-result-attribute", f"{obj.name}:{name}")
+    attr = obj.attrs.get_id(name)
+    require(attr.shape == () and attr.dtype == np.dtype(dtype),
+            "wrong-result-attribute-type", f"{obj.name}:{name}")
+    require(obj.attrs[name] == expected, "wrong-result-attribute", f"{obj.name}:{name}")
+
+
+def dataset(file, path, shape):
+    require(path in file and isinstance(file[path], h5py.Dataset), "missing-result-dataset", path)
+    value = file[path]
+    require(value.shape == shape and value.chunks is None and value.compression is None,
+            "wrong-result-shape", path)
+    return value
+
+
+def table(file, path, count, fields):
+    value = dataset(file, path, (count,))
+    require(value.dtype.names == tuple(name for name, _ in fields), "wrong-result-table-schema", path)
+    for name, dtype in fields:
+        actual = value.dtype[name]
+        if dtype == "string":
+            info = h5py.check_string_dtype(actual)
+            valid = info is not None and info.length is None and info.encoding == "utf-8"
+        else:
+            valid = actual == np.dtype(dtype)
+        require(valid, "wrong-result-table-schema", f"{path}:{name}")
+    return value[:]
+
+
+def read_results(path: Path, model):
+    """Read complete official schema/identities; defer physical finite checks."""
+    case = model["case"]
+    dynamic, nonlinear = case.startswith("hexa_dyn"), "nonlinear" in case
+    frequency = case == "hexa_freq"
+    arrays, frames, modal = [], {}, {}
+    with h5py.File(path, "r") as file:
+        require("metadata" in file and isinstance(file["metadata"], h5py.Group),
+                "missing-result-metadata", str(path))
+        metadata = file["metadata"]
+        numeric_attribute(metadata, "schema_version", "<u8", 1)
+        for name, expected in {
+            "feature_id": "hex8-structural-analysis", "analysis_procedure": CASES[case],
+            "element_formulation": "FESA-Hex8", "integration_rule": "gauss-2x2x2",
+            "mass_formulation": "lumped-reference", "source_step_name": STEP,
+            "unit_system_label": "user-consistent-unspecified",
+            "coordinate_convention": "global-cartesian; reference-strain; current-cauchy-stress",
+            "component_order": "node-major; UX,UY,UZ",
+            "source_input_identity": f"path={model['source_path']};content_identity={model['content_identity']}",
+        }.items():
+            string_attribute(metadata, name, expected)
+        string_attribute(metadata, "solver_version")
+        nodes = table(file, "model/nodes", NODE_COUNT, [
+            ("internal_node_id", "<u8"), ("instance_name", "string"),
+            ("source_label", "string"), ("coordinates", ("<f8", (3,)))])
+        node_map = {}
+        for i, row in enumerate(nodes):
+            node = label(text(row["source_label"]))
+            require(int(row["internal_node_id"]) == i and canonical(text(row["instance_name"])) == INSTANCE,
+                    "wrong-result-node-identity", str(i))
+            require(node not in node_map and node in model["nodes"], "result-node-inventory", str(node))
+            node_map[node] = i
+        require(node_map.keys() == model["nodes"].keys(), "result-node-inventory", str(path))
+        elements = table(file, "model/elements", ELEMENT_COUNT, [
+            ("internal_element_id", "<u8"), ("instance_name", "string"), ("source_label", "string"),
+            ("source_element_type", "string"), ("internal_formulation", "string"),
+            ("node_internal_ids", ("<u8", (8,))), ("solid_section_internal_id", "<u8"),
+            ("material_internal_id", "<u8")])
+        element_map = {}
+        for i, row in enumerate(elements):
+            element = label(text(row["source_label"]))
+            require(int(row["internal_element_id"]) == i and canonical(text(row["instance_name"])) == INSTANCE,
+                    "wrong-result-element-identity", str(i))
+            require(element not in element_map and element in model["elements"],
+                    "result-element-inventory", str(element))
+            require(text(row["source_element_type"]) == "C3D8"
+                    and text(row["internal_formulation"]) == "FESA-Hex8"
+                    and tuple(row["node_internal_ids"]) == tuple(node_map[n] for n in model["elements"][element])
+                    and row["solid_section_internal_id"] == 0 and row["material_internal_id"] == 0,
+                    "wrong-result-element-connectivity", str(element))
+            element_map[element] = i
+        materials = table(file, "model/solid/materials", 1, [
+            ("internal_material_id", "<u8"), ("name", "string"), ("youngs_modulus", "<f8"),
+            ("poisson_ratio", "<f8"), ("density_present", "u1"), ("density", "<f8")])
+        material = materials[0]
+        require(material["internal_material_id"] == 0 and canonical(text(material["name"])) == "MATERIAL-1"
+                and material["density_present"] == 1, "wrong-result-material-identity", str(path))
+        sections = table(file, "model/solid/sections", 1, [
+            ("internal_section_id", "<u8"), ("source_file", "string"), ("source_line", "<u8"),
+            ("source_elset", "string"), ("material_internal_id", "<u8")])
+        section = sections[0]
+        require(section["internal_section_id"] == 0 and section["material_internal_id"] == 0
+                and text(section["source_file"]) == model["source_path"]
+                and section["source_line"] == model["section_line"]
+                and canonical(text(section["source_elset"])) == canonical(model["section_elset"]),
+                "wrong-result-section-identity", str(path))
+        points = table(file, "model/solid_integration_points", ELEMENT_COUNT * 8, [
+            ("internal_element_id", "<u8"), ("point_id", "u1"),
+            ("natural_coordinates", ("<f8", (3,))), ("reference_coordinates", ("<f8", (3,)))])
+        for i, row in enumerate(points):
+            require(row["internal_element_id"] == i // 8 and row["point_id"] == i % 8 + 1,
+                    "wrong-result-point-identity", str(i))
+        mask = dataset(file, "model/nodal_constraint_mask", (NODE_COUNT, 3))
+        require(mask.dtype == np.dtype("u1"), "wrong-result-mask-type", str(path))
+        string_attribute(mask, "component_names", "UX,UY,UZ")
+        mask = mask[:]
+        expected_mask = np.array([[(node, c) in model["constrained"] for c in (1, 2, 3)]
+                                  for node in node_map], dtype="u1")
+        require(np.array_equal(mask, expected_mask), "wrong-result-constraint-identity", str(path))
+        require("steps" in file and set(file["steps"]) == {STEP}, "wrong-result-step", str(path))
+        step = file[f"steps/{STEP}"]
+        expected_step_groups = {"eigen"} if frequency else ({"frames", "accepted_frames"} if dynamic or nonlinear else {"frames"})
+        require(set(step) == expected_step_groups, "wrong-result-step-inventory", str(path))
+        if frequency:
+            modal = read_modal_tables(file, arrays, int((mask == 0).sum()))
+        else:
+            count = 101 if dynamic else (10 if nonlinear else 1)
+            accepted = None
+            if dynamic or nonlinear:
+                coordinate = "time" if dynamic else "load_factor"
+                flag = "is_initial" if dynamic else "is_factor_one"
+                accepted = table(file, f"steps/{STEP}/accepted_frames", count, [
+                    ("frame_index", "<u8"), ("increment_index", "<u8"), (coordinate, "<f8"),
+                    ("iteration_count", "<u8"), (flag, "u1")])
+                for frame, row in enumerate(accepted):
+                    increment = frame if dynamic else frame + 1
+                    value = increment * 0.1
+                    expected_flag = int(frame == (0 if dynamic else count - 1))
+                    require(row["frame_index"] == frame and row["increment_index"] == increment
+                            and row[coordinate] == value and row[flag] == expected_flag
+                            and (row["iteration_count"] == 0 if dynamic and frame == 0
+                                 else 1 <= row["iteration_count"] <= 100),
+                            "wrong-result-frame-schedule", str(frame))
+            frame_root = f"steps/{STEP}/frames"
+            require(set(file[frame_root]) == {str(i) for i in range(count)},
+                    "result-frame-inventory", str(path))
+            for frame in range(count):
+                increment = frame if dynamic else (frame + 1 if nonlinear else 10)
+                value = increment * 0.1 if dynamic or nonlinear else 1.0
+                frames[increment] = read_frame(file, frame, increment, value, dynamic, nonlinear, arrays)
+        # Identity and schema checks above deliberately finish before physical checks.
+        arrays.extend([("model/nodes/coordinates", nodes["coordinates"]),
+                       ("model/points/natural", points["natural_coordinates"]),
+                       ("model/points/reference", points["reference_coordinates"]),
+                       ("model/material", np.array([material[n] for n in ("youngs_modulus", "poisson_ratio", "density")]))])
+        evidence = {"source_input_identity": string_attribute(metadata, "source_input_identity"),
+                    "source_node_count": NODE_COUNT, "source_element_count": ELEMENT_COUNT,
+                    "accepted_frame_count": len(frames), "constrained_scalar_count": int(mask.sum()),
+                    "result_node_map": node_map, "result_element_map": element_map}
+    return {"frames": frames, "modal": modal, "nodes": nodes, "points": points,
+            "material": material, "arrays": arrays, "mask": mask, "evidence": evidence}
+
+
+def physical_dataset(file, path, shape, names, units, location, frame=None,
+                     increment=None, coordinate=None, dynamic=False, extra=None):
+    value = dataset(file, path, shape)
+    require(value.dtype == np.dtype("<f8"), "wrong-result-value-type", path)
+    attributes = {"component_names": names, "unit_dimensions": units,
+                  "coordinate_system": "global-cartesian", "location": location}
+    if frame is not None:
+        attributes["step_name"] = STEP
+    attributes.update(extra or {})
+    for name, expected in attributes.items():
+        string_attribute(value, name, expected)
+    if frame is not None:
+        numeric_attribute(value, "frame_index", "<u8", frame)
+        numeric_attribute(value, "increment_index", "<u8", increment)
+        numeric_attribute(value, "time" if dynamic else "load_factor", "<f8", coordinate)
+    return value[()]
+
+
+def read_frame(file, frame, increment, coordinate, dynamic, nonlinear, arrays):
+    root = f"steps/{STEP}/frames/{frame}"
+    require(set(file[root]) == {"nodal", "solid", "global"}, "wrong-result-frame-inventory", root)
+    definitions = {
+        "nodal/displacement": ((NODE_COUNT, 3), "U1,U2,U3", "length", "nodal", {}),
+        "nodal/reaction": ((NODE_COUNT, 3), "RF1,RF2,RF3", "force", "nodal", {}),
+        "nodal/residual": ((NODE_COUNT, 3), "R1,R2,R3", "force", "nodal", {}),
+        "solid/stress": ((ELEMENT_COUNT, 8, 6), "S11,S22,S33,S12,S13,S23", "force/length^2",
+                         "integration-point", {"stress_measure": "cauchy"}),
+        "solid/strain": ((ELEMENT_COUNT, 8, 6), "E11,E22,E33,E12,E13,E23", "dimensionless",
+                         "integration-point", {"strain_measure": "green-lagrange" if nonlinear else "infinitesimal",
+                                               "shear_convention": "tensor"}),
+        "global/strain_energy": ((), "scalar", "force*length", "global",
+                                 {"definition": "sum(W*detJ0)" if nonlinear else "0.5*u^T*K0*u"}),
+    }
+    if dynamic:
+        definitions.update({
+            "nodal/velocity": ((NODE_COUNT, 3), "V1,V2,V3", "length/time", "nodal", {}),
+            "nodal/acceleration": ((NODE_COUNT, 3), "A1,A2,A3", "length/time^2", "nodal", {"acceleration_kind": "physical"}),
+            "nodal/algorithmic_acceleration": ((NODE_COUNT, 3), "A1,A2,A3", "length/time^2", "nodal", {"acceleration_kind": "algorithmic"}),
+            "global/kinetic_energy": ((), "scalar", "force*length", "global", {"definition": "0.5*v^T*M*v"}),
+            "global/external_work": ((), "scalar", "force*length", "global", {"definition": "sum(0.5*(F_old+F_new)*delta_u)"}),
+            "global/energy_balance_error": ((), "scalar", "force*length", "global", {"definition": "T+U-T0-U0-Wext"}),
+            "global/linear_momentum": ((3,), "P1,P2,P3", "mass*length/time", "global", {}),
+            "global/angular_momentum": ((3,), "L1,L2,L3", "mass*length^2/time", "global", {"origin": "global0"}),
+        })
+    for group in ("nodal", "solid", "global"):
+        expected = {key.split("/", 1)[1] for key in definitions if key.startswith(group + "/")}
+        require(set(file[f"{root}/{group}"]) == expected, "wrong-result-quantity-inventory", root + "/" + group)
+    values = {}
+    for quantity, (shape, names, units, location, extra) in definitions.items():
+        path = root + "/" + quantity
+        values[quantity] = physical_dataset(file, path, shape, names, units, location,
+                                            frame, increment, coordinate, dynamic, extra)
+        arrays.append((path, values[quantity]))
+    return values
+
+
+def read_modal_tables(file, arrays, free_count):
+    root = f"steps/{STEP}/eigen"
+    require(set(file[root]) == {"modes", "frequencies", "vectors", "verification"},
+            "wrong-modal-result-inventory", root)
+    modes = table(file, root + "/modes", 10, [("mode_id", "<u8"), ("eigenvalue", "<f8"),
+                                                  ("classification", "u1"), ("cluster_id", "<u8")])
+    for name, component, units in (("modes", "eigenvalue", "1/time^2"),
+                                    ("frequencies", "frequency", "1/time")):
+        obj = file[root + "/" + name]
+        for attr, expected in {"component_names": component, "unit_dimensions": units,
+                               "coordinate_system": "global-cartesian", "location": "global"}.items():
+            string_attribute(obj, attr, expected)
+    frequencies = table(file, root + "/frequencies", sum(row["eigenvalue"] >= 0 for row in modes),
+                        [("mode_id", "<u8"), ("frequency", "<f8")])
+    require(list(modes["mode_id"]) == list(range(1, 11)) and (modes["classification"] <= 1).all()
+            and modes["cluster_id"][0] == 1
+            and all(b in (a, a + 1) for a, b in zip(modes["cluster_id"], modes["cluster_id"][1:])),
+            "wrong-modal-mode", root)
+    require(list(frequencies["mode_id"]) == [row["mode_id"] for row in modes if row["eigenvalue"] >= 0],
+            "wrong-modal-frequency-inventory", root)
+    vectors = physical_dataset(file, root + "/vectors", (10, NODE_COUNT, 3),
+                                "UX,UY,UZ", "mass-normalized", "nodal",
+                                extra={"normalization": "phi^T*M*phi=1"})
+    verification = root + "/verification"
+    require(set(file[verification]) == {"pairs", "eta_gram", "max_eta_orth"},
+            "wrong-modal-verification-inventory", verification)
+    pair_object = file[verification + "/pairs"]
+    require(isinstance(pair_object, h5py.Dataset) and len(pair_object.shape) == 1,
+            "wrong-modal-pair-shape", verification + "/pairs")
+    pair_count = pair_object.shape[0]
+    require(10 <= pair_count <= free_count, "wrong-modal-pair-count", str(pair_count))
+    pairs = table(file, verification + "/pairs", pair_count, [
+        ("sorted_rank", "<u8"), ("solver_column", "<u8"), ("eigenvalue", "<f8"),
+        ("modal_mass", "<f8"), ("rayleigh_quotient", "<f8"), ("eta_r", "<f8"),
+        ("eta_norm", "<f8"), ("eta_energy", "<f8")])
+    require(list(pairs["sorted_rank"]) == list(range(1, pair_count + 1))
+            and len(set(pairs["solver_column"])) == pair_count
+            and (pairs["solver_column"] < pair_count).all(), "wrong-modal-pair-identity", verification)
+    for name in ("eta_gram", "max_eta_orth"):
+        value = dataset(file, verification + "/" + name, ())
+        require(value.dtype == np.dtype("<f8"), "wrong-result-value-type", name)
+        arrays.append((verification + "/" + name, value[()]))
+    arrays.extend([(root + "/vectors", vectors), (root + "/lambda", modes["eigenvalue"]),
+                   (root + "/frequency", frequencies["frequency"])])
+    arrays.extend([(verification + "/" + name, pairs[name]) for name in pairs.dtype.names[2:]])
+    return {"modes": modes, "frequencies": frequencies, "pairs": pairs}
+
+
+def check_result_values(result, model):
+    for path, values in result["arrays"]:
+        require(np.isfinite(values).all(), "nonfinite-result-value", path)
+    for row in result["nodes"]:
+        node = label(text(row["source_label"]))
+        require(tuple(row["coordinates"]) == model["nodes"][node], "wrong-result-coordinate", str(node))
+    material = result["material"]
+    require(tuple(material[name] for name in ("youngs_modulus", "poisson_ratio")) == model["elastic"]
+            and material["density"] == model["density"], "wrong-result-material-value", model["case"])
+    for i, row in enumerate(result["points"]):
+        element = list(result["evidence"]["result_element_map"])[i // 8]
+        expected, scale = point_coordinates(model, element, i % 8 + 1)
+        require(np.array_equal(row["natural_coordinates"], NATURAL_POINTS[i % 8])
+                and (np.abs(row["reference_coordinates"] - expected) <= 16 * np.finfo(float).eps * scale).all(),
+                "result-coordinate-convention", str(i))
+    for increment, values in result["frames"].items():
+        require((values["nodal/reaction"][result["mask"] == 0] == 0.0).all(),
+                "nonzero-free-result-reaction", str(increment))
+        if model["case"].startswith("hexa_dyn") and increment == 0:
+            require((values["nodal/displacement"] == 0.0).all() and (values["nodal/velocity"] == 0.0).all(),
+                    "nonzero-initial-result-state", str(increment))
+    if result["modal"]:
+        modes, frequencies, pairs = (result["modal"][key] for key in ("modes", "frequencies", "pairs"))
+        require((np.diff(modes["eigenvalue"]) >= 0).all() and (np.diff(pairs["eigenvalue"]) >= 0).all(),
+                "wrong-modal-order", model["case"])
+        for row in frequencies:
+            expected = math.sqrt(modes[row["mode_id"] - 1]["eigenvalue"]) / (2 * math.pi)
+            require(row["frequency"] == expected, "wrong-modal-frequency", str(row["mode_id"]))
+
+
+def compare(case: str, input_path: Path, results_path: Path, reference_directory: Path = REFERENCE_DIRECTORY,
+            execution_record: Path | None = None):
+    report = {"schema": "hex8-reference-comparison-v1", "case": case, "policy": dict(POLICY),
+              "identities": {}, "precheck": {"passed": False}, "families": [],
+              "blocking_failures": [], "blocking_passed": False, "verdict": "failed"}
+    try:
+        require(case in CASES, "unknown-source-case", case)
+        quantities = ["U"] if case.startswith("hexa_dyn") else (["U", "RF", "S"] if case != "hexa_freq" else [])
+        paths = {"input": input_path, "results": results_path}
+        paths.update({quantity: reference_directory / f"{case}_{quantity}.csv" for quantity in quantities})
+        modal_path = reference_directory / "hexa_freq.txt"
+        if case == "hexa_freq" and modal_path.exists():
+            paths["Value"] = modal_path
+        for name, path in paths.items():
+            report["identities"][name] = {"path": lexical_absolute(path), "sha256": sha256(path)}
+        model = input_inventory(input_path, case)
+        report["input_copy_matches_approved_bytes"] = True
+        report["approved_input_sha256"] = SOURCE_SHA256[f"{case}.inp"]
+        report["execution"] = {"verified": False, "scope": "file comparison only; execution record not supplied"}
+        if execution_record is not None:
+            execution = json.loads(execution_record.read_text(encoding="utf-8-sig"))
+            command = execution.get("command")
+            require(execution.get("case") == case and execution.get("native_exit") == 0
+                    and type(execution.get("native_exit")) is int
+                    and execution.get("original_sha256") == report["identities"]["input"]["sha256"]
+                    and execution.get("scratch_sha256") == report["identities"]["input"]["sha256"]
+                    and execution.get("result_sha256") == report["identities"]["results"]["sha256"]
+                    and execution.get("source_unchanged") is True
+                    and execution.get("scratch_unchanged") is True
+                    and execution.get("result_exists") is True
+                    and isinstance(command, list) and len(command) == 2
+                    and all(isinstance(item, str) for item in command)
+                    and Path(lexical_absolute(Path(command[1]))) == Path(lexical_absolute(input_path)),
+                    "invalid-execution-evidence", str(execution_record))
+            report["execution"] = {"verified": True, "native_exit": 0, "command": command,
+                                   "record_path": lexical_absolute(execution_record),
+                                   "record_sha256": sha256(execution_record)}
+        records = {}
+        modal_reference = None
+        if "Value" in paths:
+            require(report["identities"]["Value"]["sha256"] == SOURCE_SHA256["hexa_freq.txt"],
+                    "reference-hash-mismatch", str(modal_path))
+            modal_reference = read_modal_reference(modal_path)
+        for quantity in quantities:
+            require(report["identities"][quantity]["sha256"] == SOURCE_SHA256[f"{case}_{quantity}.csv"],
+                    "reference-hash-mismatch", str(paths[quantity]))
+            records[quantity] = read_reference(paths[quantity], case, quantity, model)
+        result = read_results(results_path, model)
+        references = {quantity: reference_values(rows, quantity, model) for quantity, rows in records.items()}
+        check_result_values(result, model)
+        report["precheck"] = {"passed": True, **result["evidence"],
+                              "reference_row_counts": {q: len(rows) for q, rows in records.items()},
+                              "coordinate_check": "auxiliary; fixed identity; 16*float32-eps uncancelled source-coordinate budget"}
+        # Internal maps are useful to validation but unnecessary in the public report.
+        del report["precheck"]["result_node_map"], report["precheck"]["result_element_map"]
+        if case == "hexa_freq":
+            if modal_reference is None:
+                report["verdict"] = "not-run"
+                report["reason"] = "Original modal Value reference is absent; no format/quantity/unit assumption made."
+                return report
+            actual = [{"step": STEP, "mode": int(row["mode_id"]), "quantity": "eigenvalue",
+                       "unit_dimensions": "1/time^2", "value": float(row["eigenvalue"])}
+                      for row in result["modal"]["modes"]]
+            modal = compare_modal(modal_reference, actual, "eigenvalue", "1/time^2")
+            report["precheck"]["reference_row_counts"]["Value"] = len(modal_reference)
+            report["modal_reference"] = {"format": "Mode n: Value = lambda Freq = f (cycles/time)",
+                                         "step_binding": "sole Step-1 of approved hexa_freq.inp",
+                                         "quantity": "eigenvalue", "unit_dimensions": "1/time^2",
+                                         "frequency_role": "finite auxiliary; no conversion"}
+            report["families"] = modal["families"]
+            for family in report["families"]:
+                family["reference_sha256"] = report["identities"]["Value"]["sha256"]
+            report["blocking_failures"] = [f["identity"] for f in modal["families"] if not f["passed"]]
+            report["blocking_passed"] = modal["passed"]
+            report["verdict"] = "passed" if modal["passed"] else "failed"
+            return report
+        families = []
+        for increment in sorted(result["frames"]):
+            if not case.startswith("hexa_dyn") and increment != 10:
+                continue
+            values = result["frames"][increment]
+            for quantity in quantities:
+                actual, expected, identities = [], [], []
+                for key in sorted(references[quantity]):
+                    if key[0] != increment:
+                        continue
+                    _, entity, point = key
+                    index = result["evidence"]["result_element_map" if quantity == "S" else "result_node_map"][entity]
+                    row = values[{"U": "nodal/displacement", "RF": "nodal/reaction", "S": "solid/stress"}[quantity]][index]
+                    if quantity == "S":
+                        row = row[point - 1]
+                    for c, component in enumerate(COMPONENTS[quantity]):
+                        if quantity == "RF" and (entity, c + 1) not in model["constrained"]:
+                            continue
+                        identities.append({"instance": INSTANCE, "component": component,
+                                           "element_label" if quantity == "S" else "node_label": entity,
+                                           **({"integration_point": point} if quantity == "S" else {})})
+                        actual.append(float(row[c]))
+                        expected.append(references[quantity][key][c])
+                identity = {"identity": f"{case}:{STEP}:increment-{increment}:{quantity}-global",
+                            "case": case, "step": STEP, "increment": increment, "quantity": quantity,
+                            "unit_dimensions": {"U": "length", "RF": "force", "S": "force/length^2"}[quantity],
+                            "coordinate_system": "global-cartesian", "rows": identities,
+                            "reference_sha256": report["identities"][quantity]["sha256"]}
+                if case.startswith("hexa_dyn"):
+                    identity.update({"hdf5_time": increment * 0.1,
+                                     "reference_decimal_time": str(Decimal(increment) / 10)})
+                families.append(evaluate_family(identity, actual, expected))
+        report["families"] = families
+        report["blocking_failures"] = [f["identity"] for f in families if not f["passed"]]
+        report["blocking_passed"] = not report["blocking_failures"]
+        report["verdict"] = "passed" if report["blocking_passed"] else "failed"
+    except (ComparisonError, OSError, ValueError, InvalidOperation, UnicodeError, KeyError, TypeError) as error:
+        report["precheck"] = {"passed": False, "code": getattr(error, "code", "invalid-comparison-input"),
+                              "message": str(error)}
+    return report
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--case", choices=CASES, required=True)
+    parser.add_argument("--input", type=Path, required=True, help="Actual byte-identical scratch input read by the solver.")
+    parser.add_argument("--results", type=Path, required=True)
+    parser.add_argument("--references", type=Path, default=REFERENCE_DIRECTORY)
+    parser.add_argument("--execution-record", type=Path, help="Closed native CLI record with input/result hashes; never executed here.")
+    parser.add_argument("--report", type=Path, required=True)
+    args = parser.parse_args()
+    protected = {args.input.resolve(), args.results.resolve()}
+    if args.execution_record is not None:
+        protected.add(args.execution_record.resolve())
+    protected.update(path.resolve() for path in args.references.glob("*"))
+    protected.update(path.resolve() for path in REFERENCE_DIRECTORY.glob("*"))
+    report_path = args.report.resolve()
+    hardlink_alias = args.report.exists() and any(
+        path.exists() and args.report.samefile(path) for path in protected)
+    if report_path in protected or hardlink_alias:
+        parser.error("The report path must not overwrite an input, reference or result.")
+    report = compare(args.case, args.input, args.results, args.references, args.execution_record)
+    args.report.parent.mkdir(parents=True, exist_ok=True)
+    args.report.write_text(json.dumps(report, sort_keys=True, indent=2, allow_nan=False) + "\n", encoding="utf-8")
+    print(json.dumps({"case": args.case, "verdict": report["verdict"], "precheck": report["precheck"],
+                      "blocking_failure_count": len(report["blocking_failures"]), "report": str(args.report)}, allow_nan=False))
+    return 2 if not report["precheck"]["passed"] or report["verdict"] == "not-run" else (0 if report["blocking_passed"] else 1)
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

scripts/verification/compare_hht_dynamic.py

@@ -0,0 +1,436 @@
+"""Compare the approved 101-frame HHT reference without modifying artifacts.
+
+Run with ``uv run --with h5py python -B scripts/verification/compare_hht_dynamic.py``.
+Only the documented displacement/rotation reference quantities are compared;
+this is not an independent validation of inertia, energy, or convergence.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+from decimal import Decimal, InvalidOperation
+import hashlib
+import json
+import math
+from pathlib import Path
+import re
+from typing import Any
+
+import h5py
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[2]
+CASE = "shellNL-dyn-meshfine"
+REFERENCE_DIRECTORY = ROOT / "reference" / CASE
+NODE_LABELS = (2, 21, 22, 69, 70, 165, 309, 442, 586)
+INSTANCE = "PART-1-1"
+STEP = "Step-1"
+NODE_COUNT = 729
+FRAME_COUNT = 101
+INCREMENT = 0.1
+COMPONENTS = ("U-U1", "U-U2", "U-U3", "UR-UR1", "UR-UR2", "UR-UR3")
+CSV_COLUMNS = ("Frame", "Part Instance Name", "Node Label", *COMPONENTS)
+FRAME_PATTERN = re.compile(
+    r"Increment\s+(\d+):\s*Step Time\s*=\s*"
+    r"([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)"
+)
+POLICY = {
+    "policy_id": "common-reference-tolerance",
+    "near_zero_ratio": 0.01,
+    "relative_tolerance": 0.05,
+    "relative_rms_tolerance": 0.01,
+}
+
+
+class ComparisonError(ValueError):
+    """An identity/schema/finite failure that must precede tolerance evaluation."""
+
+    def __init__(self, code: str, message: str):
+        super().__init__(message)
+        self.code = code
+
+
+def require(condition: bool, code: str, message: str) -> None:
+    if not condition:
+        raise ComparisonError(code, message)
+
+
+def positive_label(text: str) -> int:
+    require(bool(re.fullmatch(r"\d+", text)), "invalid-node-label", repr(text))
+    label = int(text)
+    require(0 < label <= 2**63 - 1, "invalid-node-label", repr(text))
+    return label
+
+
+def canonical_instance(text: str) -> str:
+    return text.translate(str.maketrans("abcdefghijklmnopqrstuvwxyz", "ABCDEFGHIJKLMNOPQRSTUVWXYZ"))
+
+
+def input_inventory(path: Path) -> tuple[dict[int, tuple[float, ...]], str]:
+    """Read identity and schedule authority from this specific approved INP."""
+    raw = path.read_bytes()
+    nodes: dict[int, tuple[float, ...]] = {}
+    parts: list[str] = []
+    instances: list[tuple[str, str]] = []
+    steps: list[dict[str, str]] = []
+    dynamic_rows: list[tuple[float, float]] = []
+    keyword = ""
+    in_part = False
+    for line_number, line in enumerate(raw.decode("utf-8-sig").splitlines(), 1):
+        text = line.strip()
+        if not text or text.startswith("**"):
+            continue
+        fields = [field.strip() for field in text.split(",")]
+        if text.startswith("*"):
+            keyword = fields[0].upper()
+            parameters = dict(
+                (name.upper(), value.strip())
+                for field in fields[1:] if "=" in field
+                for name, value in [field.split("=", 1)]
+            )
+            if keyword == "*PART":
+                in_part = True
+                parts.append(canonical_instance(parameters.get("NAME", "")))
+            elif keyword == "*END PART":
+                in_part = False
+            elif keyword == "*INSTANCE":
+                instances.append((canonical_instance(parameters.get("NAME", "")),
+                                  canonical_instance(parameters.get("PART", ""))))
+            elif keyword == "*STEP":
+                steps.append(parameters)
+            elif keyword == "*DYNAMIC":
+                require([field.upper() for field in fields] == ["*DYNAMIC", "DIRECT"],
+                        "wrong-input-procedure", f"line {line_number}")
+            continue
+        if keyword == "*NODE":
+            require(in_part and len(fields) == 4, "invalid-input-node", f"line {line_number}")
+            label = positive_label(fields[0])
+            require(label not in nodes, "duplicate-input-node", str(label))
+            coordinates = tuple(float(value) for value in fields[1:])
+            require(all(math.isfinite(value) for value in coordinates),
+                    "nonfinite-input-node", str(label))
+            nodes[label] = coordinates
+        elif keyword == "*DYNAMIC":
+            if fields[-1] == "":
+                fields.pop()
+            require(len(fields) == 2, "wrong-input-schedule", f"line {line_number}")
+            dynamic_rows.append((float(fields[0]), float(fields[1])))
+    require(parts == ["PART-1"] and instances == [(INSTANCE, "PART-1")],
+            "wrong-input-instance", "Expected the original single identity instance.")
+    require(len(nodes) == NODE_COUNT and set(NODE_LABELS) <= nodes.keys(),
+            "wrong-input-node-inventory", f"Expected {NODE_COUNT} source nodes.")
+    require(len(steps) == 1 and steps[0].get("NAME") == STEP
+            and steps[0].get("NLGEOM", "").upper() == "YES"
+            and steps[0].get("INC") == "1000", "wrong-input-step", str(steps))
+    require(dynamic_rows == [(INCREMENT, 10.0)], "wrong-input-schedule", str(dynamic_rows))
+    value = 14695981039346656037
+    for byte in raw:
+        value = ((value ^ byte) * 1099511628211) & (2**64 - 1)
+    return nodes, f"fnv1a64:{value:016x}"
+
+
+def read_reference(path: Path) -> dict[tuple[int, int], tuple[float, ...]]:
+    rows: dict[tuple[int, int], tuple[float, ...]] = {}
+    with path.open(encoding="utf-8-sig", newline="") as stream:
+        reader = csv.reader(stream)
+        header = [field.strip() for field in next(reader, [])]
+        require(len(header) == len(set(header)), "duplicate-reference-column", str(header))
+        require(set(CSV_COLUMNS) <= set(header), "missing-reference-column", str(header))
+        for line_number, fields in enumerate(reader, 2):
+            if not fields or all(not field.strip() for field in fields):
+                continue
+            require(len(fields) == len(header), "invalid-reference-row", f"line {line_number}")
+            row = dict(zip(header, (field.strip() for field in fields)))
+            match = FRAME_PATTERN.fullmatch(row["Frame"])
+            require(match is not None, "invalid-reference-frame", f"line {line_number}")
+            increment = int(match.group(1))
+            time = Decimal(match.group(2))
+            require(time.is_finite() and time == Decimal(increment) / Decimal(10),
+                    "wrong-reference-time", f"line {line_number}: {row['Frame']}")
+            require(canonical_instance(row["Part Instance Name"]) == INSTANCE,
+                    "wrong-reference-instance", f"line {line_number}")
+            label = positive_label(row["Node Label"])
+            key = (increment, label)
+            require(key not in rows, "duplicate-reference-row", str(key))
+            values = tuple(float(row[component]) for component in COMPONENTS)
+            require(all(math.isfinite(value) for value in values),
+                    "nonfinite-reference-value", str(key))
+            rows[key] = values
+    expected = {(frame, label) for frame in range(FRAME_COUNT) for label in NODE_LABELS}
+    require(rows.keys() == expected, "reference-row-inventory",
+            f"missing={sorted(expected - rows.keys())}; extra={sorted(rows.keys() - expected)}")
+    return rows
+
+
+def text_value(value: Any) -> str:
+    return value.decode("utf-8") if isinstance(value, bytes) else str(value)
+
+
+def string_attribute(obj: Any, name: str, expected: str | None = None) -> str:
+    require(name in obj.attrs, "missing-hdf5-attribute", f"{obj.name}@{name}")
+    attribute = obj.attrs.get_id(name)
+    string = h5py.check_string_dtype(attribute.dtype)
+    require(attribute.shape == () and string is not None and string.length is None
+            and string.encoding == "utf-8", "wrong-hdf5-attribute-type", f"{obj.name}@{name}")
+    value = text_value(obj.attrs[name])
+    require(bool(value) and (expected is None or value == expected),
+            "wrong-hdf5-attribute-value", f"{obj.name}@{name}: {value!r}")
+    return value
+
+
+def numeric_attribute(obj: Any, name: str, dtype: str, expected: int | float) -> None:
+    require(name in obj.attrs, "missing-hdf5-attribute", f"{obj.name}@{name}")
+    attribute = obj.attrs.get_id(name)
+    require(attribute.shape == () and attribute.dtype == np.dtype(dtype),
+            "wrong-hdf5-attribute-type", f"{obj.name}@{name}")
+    value = obj.attrs[name]
+    require(math.isfinite(value) and value == expected, "wrong-hdf5-attribute-value",
+            f"{obj.name}@{name}: {value!r}")
+
+
+def dataset(file: h5py.File, path: str, shape: tuple[int, ...]) -> h5py.Dataset:
+    require(path in file and isinstance(file[path], h5py.Dataset), "missing-hdf5-dataset", path)
+    result = file[path]
+    require(result.shape == shape and result.chunks is None and result.compression is None,
+            "wrong-hdf5-dataset-layout", path)
+    return result
+
+
+def read_results(path: Path, nodes: dict[int, tuple[float, ...]],
+                 content_identity: str) -> tuple[list[np.ndarray], dict[str, Any]]:
+    """Precheck every result node/frame/component, then project by source key."""
+    projected: list[np.ndarray] = []
+    with h5py.File(path, "r") as file:
+        require("metadata" in file and isinstance(file["metadata"], h5py.Group),
+                "missing-hdf5-metadata", "/metadata")
+        metadata = file["metadata"]
+        for name, expected in {
+            "feature_id": "hht-implicit-dynamics", "analysis_procedure": "hht-implicit-dynamics",
+            "mass_formulation": "nodal-lumped-kinetic", "internal_formulation": "FESA-MITC4-total-lagrangian",
+            "initial_state": "zero-q-v-step-load-0+",
+            "acceleration_convention": "physical-endpoint-and-auxiliary-separated",
+            "source_step_name": STEP, "step_name": STEP, "unit_system_label": "consistent-input-units",
+        }.items():
+            string_attribute(metadata, name, expected)
+        string_attribute(metadata, "solver_version")
+        source_identity = string_attribute(metadata, "source_input_identity")
+        require(source_identity.endswith(";content_identity=" + content_identity),
+                "wrong-input-content-identity", source_identity)
+        require("schema_version" in metadata.attrs, "missing-hdf5-attribute", "/metadata@schema_version")
+        version_attribute = metadata.attrs.get_id("schema_version")
+        require(version_attribute.shape == () and version_attribute.dtype == np.dtype("<u8"),
+                "wrong-hdf5-attribute-type", "/metadata@schema_version")
+        schema_version = int(metadata.attrs["schema_version"])
+        require(schema_version in (1, 2), "wrong-hdf5-attribute-value", "/metadata@schema_version")
+        if schema_version == 2:
+            string_attribute(metadata, "time_integration_formulation",
+                             "reduced-coordinate-HHT-with-physical-rate-reconstruction")
+        elif "time_integration_formulation" in metadata.attrs:
+            string_attribute(metadata, "time_integration_formulation", "full-coordinate-HHT")
+        for name, dtype, value in (
+            ("schema_version", "<u8", schema_version), ("increment_count", "<u8", 100),
+            ("maximum_increments", "<u8", 1000), ("alpha", "<f8", -0.05),
+            ("beta", "<f8", 0.275625), ("gamma", "<f8", 0.55),
+            ("initial_increment", "<f8", INCREMENT), ("time_period", "<f8", 10.0),
+        ):
+            numeric_attribute(metadata, name, dtype, value)
+        table = dataset(file, "model/nodes", (NODE_COUNT,))
+        require(table.dtype.names == ("internal_node_id", "instance_name", "source_label", "coordinates"),
+                "wrong-node-schema", str(table.dtype))
+        require(table.dtype["internal_node_id"] == np.dtype("<u8")
+                and table.dtype["coordinates"] == np.dtype(("<f8", (3,))),
+                "wrong-node-schema", str(table.dtype))
+        for name in ("instance_name", "source_label"):
+            string = h5py.check_string_dtype(table.dtype[name])
+            require(string is not None and string.length is None and string.encoding == "utf-8",
+                    "wrong-node-schema", name)
+        node_indices: dict[int, int] = {}
+        for index, node in enumerate(table[:]):
+            require(int(node["internal_node_id"]) == index, "wrong-internal-node-identity", str(index))
+            instance = canonical_instance(text_value(node["instance_name"]))
+            label = positive_label(text_value(node["source_label"]))
+            require(instance == INSTANCE, "wrong-result-instance", instance)
+            require(label not in node_indices, "duplicate-result-node", str(label))
+            require(label in nodes, "extra-result-node", str(label))
+            require(np.isfinite(node["coordinates"]).all(), "nonfinite-result-node", str(label))
+            require(tuple(node["coordinates"]) == nodes[label], "wrong-result-coordinate", str(label))
+            node_indices[label] = index
+        require(node_indices.keys() == nodes.keys(), "result-node-inventory", "Source node inventories differ.")
+        require("steps" in file and isinstance(file["steps"], h5py.Group)
+                and set(file["steps"]) == {STEP}, "wrong-result-step", "Expected one accepted step.")
+        accepted = dataset(file, f"steps/{STEP}/accepted_frames", (FRAME_COUNT,))
+        expected_fields = (("frame_index", "<u8"), ("increment_index", "<u8"),
+                           ("time", "<f8"), ("iteration_count", "<u8"), ("is_initial", "u1"))
+        require(accepted.dtype.names == tuple(name for name, _ in expected_fields)
+                and all(accepted.dtype[name] == np.dtype(dtype) for name, dtype in expected_fields),
+                "wrong-frame-table-schema", str(accepted.dtype))
+        frame_path = f"steps/{STEP}/frames"
+        require(frame_path in file and isinstance(file[frame_path], h5py.Group)
+                and set(file[frame_path]) == {str(i) for i in range(FRAME_COUNT)},
+                "result-frame-inventory", "Missing, extra, or noncanonical frame group.")
+        selection = [node_indices[label] for label in NODE_LABELS]
+        for frame, row in enumerate(accepted[:]):
+            scheduled_time = float(frame) * INCREMENT
+            require(int(row["frame_index"]) == frame and int(row["increment_index"]) == frame
+                    and math.isfinite(row["time"]) and float(row["time"]) == scheduled_time
+                    and int(row["is_initial"]) == int(frame == 0)
+                    and 0 <= int(row["iteration_count"]) <= 100
+                    and (frame != 0 or (int(row["iteration_count"]) == 0 and not np.signbit(row["time"]))),
+                    "wrong-result-frame-time", str(frame))
+            displacement = dataset(file, f"{frame_path}/{frame}/nodal/displacement", (NODE_COUNT, 6))
+            require(displacement.dtype == np.dtype("<f8"), "wrong-displacement-type", str(frame))
+            for name, value in {
+                "component_names": "UX,UY,UZ,URX,URY,URZ",
+                "component_unit_dimensions": "length,length,length,1,1,1",
+                "coordinate_system": "global-cartesian", "location": "nodal",
+                "step_name": STEP, "source_step_name": STEP,
+                "rotation_representation": "total-rotation-vector-chart",
+                "quantity_semantics": "global-translation-and-total-rotation-vector",
+            }.items():
+                string_attribute(displacement, name, value)
+            for name, dtype, value in (
+                ("schema_version", "<u8", schema_version), ("frame_index", "<u8", frame),
+                ("increment_index", "<u8", frame), ("time", "<f8", scheduled_time),
+                ("is_initial", "u1", int(frame == 0)),
+            ):
+                numeric_attribute(displacement, name, dtype, value)
+            if schema_version == 2:
+                auxiliary = dataset(file, f"{frame_path}/{frame}/nodal/algorithmic_acceleration", (NODE_COUNT, 6))
+                require(auxiliary.dtype == np.dtype("<f8"), "wrong-auxiliary-type", str(frame))
+                for name, value in {
+                    "component_names": "AUX1,AUX2,AUX3,AUX4,AUX5,AUX6",
+                    "component_unit_dimensions": "length/time^2,length/time^2,length/time^2,1/time^2,1/time^2,1/time^2",
+                    "coordinate_system": "fixed-initial-global-chart", "location": "nodal",
+                    "step_name": STEP, "source_step_name": STEP,
+                    "rotation_representation": "total-rotation-vector-chart",
+                    "quantity_semantics": "full-chart-HHT-filter-combination-not-predictor-driver",
+                }.items():
+                    string_attribute(auxiliary, name, value)
+                for name, dtype, value in (
+                    ("schema_version", "<u8", schema_version), ("frame_index", "<u8", frame),
+                    ("increment_index", "<u8", frame), ("time", "<f8", scheduled_time),
+                    ("is_initial", "u1", int(frame == 0)),
+                ):
+                    numeric_attribute(auxiliary, name, dtype, value)
+                require(np.isfinite(auxiliary[:]).all(), "nonfinite-auxiliary-value", str(frame))
+            values = displacement[:]
+            require(np.isfinite(values).all(), "nonfinite-result-value", str(frame))
+            if frame == 0:
+                require((values == 0.0).all(), "nonzero-initial-displacement", "Frame zero must have q0=0.")
+            projected.append(values[selection])
+    return projected, {"source_input_identity": source_identity, "source_node_count": NODE_COUNT,
+                       "accepted_frame_count": FRAME_COUNT, "checked_result_scalars": NODE_COUNT * FRAME_COUNT * 6,
+                       "projected_node_labels": list(NODE_LABELS)}
+
+
+def evaluate_family(increment: int, quantity: str, result: list[float],
+                    reference: list[float]) -> dict[str, Any]:
+    """Apply the common policy exactly, including zero-scale and family RMS."""
+    require(len(result) == len(reference) and len(result) > 0, "scalar-inventory", quantity)
+    require(all(math.isfinite(value) for value in (*result, *reference)),
+            "nonfinite-comparison-value", quantity)
+    scale = max(abs(value) for value in reference)
+    band = POLICY["near_zero_ratio"] * scale
+    rows = []
+    error_norm = 0.0
+    for index, (actual, expected) in enumerate(zip(result, reference)):
+        error = abs(actual - expected)
+        require(math.isfinite(error), "nonfinite-comparison-value", quantity)
+        error_norm = math.hypot(error_norm, error)
+        require(math.isfinite(error_norm), "nonfinite-comparison-value", quantity)
+        relative_error = None
+        if scale == 0.0:
+            branch, threshold, passed = "zero-scale-exact", 0.0, error == 0.0
+        elif abs(expected) <= band:
+            branch, threshold, passed = "near-zero", band, error <= band
+        else:
+            branch = "relative"
+            threshold = POLICY["relative_tolerance"] * abs(expected)
+            relative_error = error / abs(expected)
+            require(math.isfinite(relative_error), "nonfinite-comparison-value", quantity)
+            passed = relative_error <= POLICY["relative_tolerance"]
+        rows.append({"instance": INSTANCE, "node_label": NODE_LABELS[index // 3],
+                     "component": f"{quantity}{index % 3 + 1}", "fesa_value": actual,
+                     "reference_value": expected, "absolute_error": error,
+                     "relative_error": relative_error, "branch": branch,
+                     "threshold": threshold, "passed": passed})
+    rows_passed = all(row["passed"] for row in rows)
+    rms = error_norm / math.sqrt(len(rows)) / scale if scale else None
+    require(rms is None or math.isfinite(rms), "nonfinite-comparison-value", quantity)
+    rms_passed = rows_passed if scale == 0.0 else rms <= POLICY["relative_rms_tolerance"]
+    return {"identity": f"{quantity}-global:{STEP}:increment-{increment}", "case": CASE,
+            "step": STEP, "increment": increment, "hdf5_time": float(increment) * INCREMENT,
+            "reference_decimal_time": str(Decimal(increment) / Decimal(10)),
+            "quantity": quantity, "components": [f"{quantity}{i}" for i in (1, 2, 3)],
+            "unit_dimension": "length" if quantity == "U" else "rad",
+            "coordinate_system": "global-cartesian" if quantity == "U" else "fixed-initial-global-total-rotation-vector",
+            "impact": "blocking" if quantity == "U" else "warning-only", "row_count": len(rows),
+            "reference_scale": scale, "near_zero_band": band, "rows": rows,
+            "rows_passed": rows_passed, "relative_rms": rms, "rms_passed": rms_passed,
+            "diagnostic_code": "zero-reference-scale-nonzero-error" if scale == 0.0 and not rows_passed else None,
+            "passed": rows_passed and rms_passed}
+
+
+def compare(input_path: Path, reference_path: Path, results_path: Path) -> dict[str, Any]:
+    report: dict[str, Any] = {
+        "schema": "hht-reference-comparison-v1", "case": CASE, "policy": dict(POLICY),
+        "identities": {}, "precheck": {"passed": False}, "families": [],
+        "blocking_failures": [], "warnings": [], "blocking_passed": False, "verdict": "failed",
+    }
+    try:
+        for name, path in (("input", input_path), ("reference", reference_path), ("results", results_path)):
+            report["identities"][name] = {"path": str(path.absolute())}
+            with path.open("rb") as stream:
+                report["identities"][name]["sha256"] = hashlib.file_digest(stream, "sha256").hexdigest()
+        nodes, content_identity = input_inventory(input_path)
+        reference = read_reference(reference_path)
+        result, evidence = read_results(results_path, nodes, content_identity)
+        families = []
+        for frame in range(FRAME_COUNT):
+            for quantity, first in (("U", 0), ("UR", 3)):
+                actual = [float(value) for row in result[frame] for value in row[first:first + 3]]
+                expected = [value for label in NODE_LABELS for value in reference[(frame, label)][first:first + 3]]
+                families.append(evaluate_family(frame, quantity, actual, expected))
+        report["precheck"] = {"passed": True, "reference_row_count": len(reference),
+                              "reference_frame_count": FRAME_COUNT, **evidence}
+        report["families"] = families
+        for family in families:
+            if family["passed"]:
+                continue
+            if family["impact"] == "blocking":
+                report["blocking_failures"].append(family["identity"])
+            else:
+                if not family["rows_passed"]:
+                    report["warnings"].append({"family": family["identity"], "code": "row-tolerance-failure"})
+                if not family["rms_passed"]:
+                    report["warnings"].append({"family": family["identity"], "code": "relative-rms-failure"})
+        report["blocking_passed"] = not report["blocking_failures"]
+        report["verdict"] = "passed" if report["blocking_passed"] else "failed"
+    except (ComparisonError, OSError, ValueError, InvalidOperation, UnicodeError) as error:
+        report["precheck"] = {"passed": False, "code": getattr(error, "code", "invalid-comparison-input"),
+                              "message": str(error)}
+    return report
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--input", type=Path, default=REFERENCE_DIRECTORY / f"{CASE}.inp")
+    parser.add_argument("--reference", type=Path, default=REFERENCE_DIRECTORY / f"{CASE} displacements.csv")
+    parser.add_argument("--results", type=Path)
+    parser.add_argument("--report", type=Path, required=True)
+    args = parser.parse_args()
+    results = args.results or args.input.with_name(args.input.stem + "_results.h5")
+    if args.report.resolve() in {args.input.resolve(), args.reference.resolve(), results.resolve()}:
+        parser.error("The report path must not overwrite an input, reference, or solver result.")
+    report = compare(args.input, args.reference, results)
+    args.report.parent.mkdir(parents=True, exist_ok=True)
+    args.report.write_text(json.dumps(report, indent=2, sort_keys=True, allow_nan=False) + "\n", encoding="utf-8")
+    print(json.dumps({"verdict": report["verdict"], "precheck": report["precheck"],
+                      "blocking_failure_count": len(report["blocking_failures"]),
+                      "warning_count": len(report["warnings"]), "report": str(args.report)}, allow_nan=False))
+    return 2 if not report["precheck"]["passed"] else (0 if report["blocking_passed"] else 1)
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tests/test_cpp_policy_contract.py

@@ -96,3 +96,19 @@ def test_p_style_002_msvc_targets_enable_standard_exception_unwinding():
         assert "/EHsc" in compile_options.group("body"), (
             f"{target} must enable standard C++ exception unwinding under MSVC"
         )
+
+
+    # The maintained solution uses the shared native property sheet; both paths
+    # must keep standard exception unwinding. Doxygen remains optional above.
+    import xml.etree.ElementTree as ET
+
+    ns = {"m": "http://schemas.microsoft.com/developer/msbuild/2003"}
+    props = ET.parse(ROOT / "msvc" / "Fesa.Common.props")
+    assert props.findtext(".//m:ExceptionHandling", namespaces=ns) == "Sync"
+    for target in ("fesa_solver", "fesa_cli", "fesa_unit_tests",
+                   "fesa_integration_tests", "fesa_reference_tests"):
+        project = ET.parse(ROOT / "msvc" / f"{target}.vcxproj")
+        imports = [item.attrib.get("Project") for item in project.findall("m:Import", ns)]
+        assert "Fesa.Common.props" in imports, target
+        overrides = project.findall(".//m:ExceptionHandling", ns)
+        assert all(item.text == "Sync" for item in overrides), target

tests/unit/analysis/eigen_reduction_test.cpp

@@ -1,3 +1,4 @@
+#include <limits>
 #include "analysis/eigen_reduction.h"
 
 #include <gtest/gtest.h>
@@ -46,6 +47,37 @@ fesa::ShellNodeInitialFrame TiltedFrame() {
   return {0U, {0.6, 0.0, 0.8}, {0.0, 1.0, 0.0}, {-0.8, 0.0, 0.6}};
 }
 
+TEST(EigenReduction, RejectsNonfiniteUsedBasisAndPreservesPriorCandidate) {
+  const fesa::ShellNodeInitialFrame frame{0U, {0, 0, 1}, {1, 0, 0}, {0, 1, 0}};
+  const auto mass = NodalMass(frame.director);
+  fesa::Matrix stiffness(6U, 6U);
+  for (std::size_t i = 0; i < 5U; ++i) stiffness(i, i) = mass(i, i) * 7.0;
+  stiffness(5U, 5U) = .25;
+  const auto matrix = Sparse(stiffness);
+  const std::vector<std::uint8_t> mask(6U, 0U);
+  const auto prior = fesa::eigen_internal::BuildEigenReduction(
+      matrix, {mass}, {frame}, mask, {.25}, 2.0);
+  ASSERT_TRUE(prior.HasValue());
+  const auto prior_values = prior.Value().whitened_operator.Values();
+  for (const double invalid : {std::numeric_limits<double>::quiet_NaN(),
+                                std::numeric_limits<double>::infinity()}) {
+    for (const bool first_tangent : {false, true}) {
+      auto bad_frame = frame;
+      (first_tangent ? bad_frame.tangent_a : bad_frame.tangent_b)[0] = invalid;
+      const auto rejected = fesa::eigen_internal::BuildEigenReduction(
+          matrix, {mass}, {bad_frame}, mask, {.25}, 2.0);
+      ASSERT_FALSE(rejected.HasValue());
+      EXPECT_EQ(rejected.GetStatus().Category(), fesa::FailureCategory::kSolver);
+      EXPECT_EQ(rejected.GetStatus().Diagnostics()[0].code,
+                "invalid-eigen-mass-space");
+      EXPECT_EQ(prior.Value().whitened_operator.Values(), prior_values);
+    }
+  }
+  const auto rerun = fesa::eigen_internal::BuildEigenReduction(
+      matrix, {mass}, {frame}, mask, {.25}, 2.0);
+  ASSERT_TRUE(rerun.HasValue());
+  EXPECT_EQ(rerun.Value().whitened_operator.Values(), prior_values);
+}
 }  // namespace
 
 TEST(EigenReduction, KeepsPositivePartialGlobalRotationAndOriginalPencil) {

tests/unit/analysis/hht_reduced_shell_test.cpp

@@ -5,6 +5,8 @@
 #include <cmath>
 #include <cfenv>
 #include <cstddef>
+#include <cstring>
+#include <functional>
 #include <iostream>
 #include <memory>
 #include <limits>
@@ -19,6 +21,7 @@
 #endif
 
 #include "analysis/hht_reduced_internal.h"
+#include "math/dense_blas_internal.h"
 #include "fesa/analysis/analysis_model.h"
 #include "fesa/analysis/hht_time_integration.h"
 #include "fesa/assembly/parallel_for.h"
@@ -26,6 +29,71 @@
 
 namespace {
 
+class RestoreTrialEnvironment {
+ public:
+  RestoreTrialEnvironment()
+      : rounding_(std::fegetround()),
+        dynamic_(mkl_get_dynamic()),
+        tls_(mkl_set_num_threads_local(0)) {
+#if defined(_M_X64)
+    mxcsr_ = _mm_getcsr();
+    _mm_setcsr((_mm_getcsr() | 0x1F80U) & ~0x8040U);
+#endif
+    static_cast<void>(std::fesetround(FE_TONEAREST));
+    mkl_set_dynamic(0);
+  }
+  ~RestoreTrialEnvironment() {
+    mkl_set_dynamic(dynamic_);
+    static_cast<void>(mkl_set_num_threads_local(tls_));
+    static_cast<void>(std::fesetround(rounding_));
+#if defined(_M_X64)
+    _mm_setcsr(mxcsr_);
+#endif
+  }
+
+ private:
+  int rounding_;
+  int dynamic_;
+  int tls_;
+#if defined(_M_X64)
+  unsigned int mxcsr_;
+#endif
+};
+
+void ExpectSameBits(const double left, const double right) {
+  EXPECT_EQ(std::memcmp(&left, &right, sizeof(double)), 0);
+}
+
+void ExpectSameOperators(
+    const fesa::hht_reduced_internal::ConstraintOperators& left,
+    const fesa::hht_reduced_internal::ConstraintOperators& right) {
+  ASSERT_EQ(left.rows.size(), right.rows.size());
+  for (std::size_t row = 0U; row < left.rows.size(); ++row) {
+    const auto& a = left.rows[row];
+    const auto& b = right.rows[row];
+    ExpectSameBits(a.value, b.value);
+    ExpectSameBits(a.value_scale, b.value_scale);
+    ExpectSameBits(a.time_first, b.time_first);
+    ExpectSameBits(a.time_second, b.time_second);
+    for (const auto& pair : {std::make_pair(&a.gradient, &b.gradient),
+                             std::make_pair(&a.time_gradient, &b.time_gradient),
+                             std::make_pair(&a.third_directional_gradient,
+                                            &b.third_directional_gradient)}) {
+      ASSERT_EQ(pair.first->size(), pair.second->size());
+      for (std::size_t i = 0U; i < pair.first->size(); ++i) {
+        EXPECT_EQ((*pair.first)[i].dof, (*pair.second)[i].dof);
+        ExpectSameBits((*pair.first)[i].value, (*pair.second)[i].value);
+      }
+    }
+    ASSERT_EQ(a.hessian.size(), b.hessian.size());
+    for (std::size_t i = 0U; i < a.hessian.size(); ++i) {
+      EXPECT_EQ(a.hessian[i].row, b.hessian[i].row);
+      EXPECT_EQ(a.hessian[i].column, b.hessian[i].column);
+      ExpectSameBits(a.hessian[i].value, b.hessian[i].value);
+    }
+  }
+}
+
 struct Fixture {
   std::unique_ptr<fesa::Domain> domain;
   fesa::DofManager dofs;
@@ -394,6 +462,479 @@ TEST(HhtReducedChart, InvertibleGraphDoesNotCertifySelectorTransversality) {
       assembly, chart.Value(), solver).HasValue());
 }
 
+TEST(HhtReducedOperators,
+     CompleteTrialConstraintDataMatchesFreshAndOwnsCopies) {
+  RestoreTrialEnvironment restore;
+  const auto fixture = MakeFixture(false, true);
+  const auto shells =
+      fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  auto q = zero;
+  auto v = zero;
+  auto force = zero;
+  auto rate = zero;
+  auto acceleration = zero;
+  for (std::size_t i = 0U; i < zero.Size(); ++i) {
+    q[i] = 0.001 * static_cast<double>((i % 5U) + 1U);
+    v[i] = 0.002 * static_cast<double>((i % 7U) + 1U);
+    if (i % 6U >= 3U) {
+      force[i] = 0.03;
+      rate[i] = -0.02;
+      acceleration[i] = 0.01;
+    }
+  }
+  auto assembled = Assemble(fixture, shells.Value(), q, v, zero, force, true);
+  ASSERT_TRUE(assembled.HasValue());
+  ASSERT_TRUE(assembled.Value().trial_constraints.has_value());
+  ASSERT_EQ(assembled.Value().trial_constraints->shells.size(),
+            shells.Value().size());
+  const auto chart = fesa::hht_reduced_internal::BuildFixedChart(
+      fixture.dofs, assembled.Value().rotations);
+  ASSERT_TRUE(chart.HasValue());
+  for (std::size_t e = 0U; e < shells.Value().size(); ++e) {
+    const auto& saved = assembled.Value().trial_constraints->shells[e];
+    ASSERT_TRUE(saved.directional_internal);
+    fesa::NonlinearMitc4TrialDofs trial{};
+    fesa::Vector direction{24U};
+    for (std::size_t node = 0U; node < 4U; ++node) {
+      const auto first = 6U * node;
+      trial.displacement[node] =
+          fesa::Vector3{saved.q[first], saved.q[first + 1U], saved.q[first + 2U]};
+      trial.rotation_vector[node] = fesa::Vector3{
+          saved.q[first + 3U], saved.q[first + 4U], saved.q[first + 5U]};
+    }
+    for (std::size_t i = 0U; i < 24U; ++i) direction[i] = saved.v[i];
+    const auto response = shells.Value()[e].internal->EvaluateTrial(trial, &direction);
+    ASSERT_TRUE(response.HasValue());
+    for (std::size_t node = 0U; node < 4U; ++node) {
+      ExpectSameBits(saved.drilling_constraint[node],
+                     response.Value().drilling_constraint[node]);
+      for (std::size_t i = 0U; i < 24U; ++i) {
+        ExpectSameBits(saved.drilling_constraint_gradient(node, i),
+                       response.Value().drilling_constraint_gradient(node, i));
+        ExpectSameBits(
+            saved.drilling_constraint_third_directional_gradient[node][i],
+            response.Value()
+                .drilling_constraint_third_directional_gradient[node][i]);
+        for (std::size_t j = 0U; j < 24U; ++j)
+          ExpectSameBits(
+              saved.drilling_constraint_hessian[node](i, j),
+              response.Value().drilling_constraint_hessian[node](i, j));
+      }
+    }
+  }
+  const auto evaluate = [&](const fesa::hht_internal::DynamicAssembly& value) {
+    return fesa::hht_reduced_internal::EvaluateConstraintOperators(
+        *fixture.domain, fixture.dofs, shells.Value(), value, chart.Value(),
+        rate, acceleration);
+  };
+  auto fresh = assembled.Value();
+  fresh.trial_constraints.reset();
+  const auto reference = evaluate(fresh);
+  const auto reused = evaluate(assembled.Value());
+  ASSERT_TRUE(reference.HasValue());
+  ASSERT_TRUE(reused.HasValue());
+  ExpectSameOperators(reused.Value(), reference.Value());
+  auto copy = assembled.Value();
+  copy.trial_constraints->shells[0U].drilling_constraint[0U] += 7.0;
+  // Test-only finite marker proves the selected route actually reads its owner.
+  const auto marker = evaluate(copy);
+  ASSERT_TRUE(marker.HasValue());
+  EXPECT_NE(marker.Value().rows[0U].value, reused.Value().rows[0U].value);
+  EXPECT_NE(
+      copy.trial_constraints->shells[0U].drilling_constraint[0U],
+      assembled.Value().trial_constraints->shells[0U].drilling_constraint[0U]);
+  auto moved = std::move(assembled.Value());
+  const auto moved_result = evaluate(moved);
+  ASSERT_TRUE(moved_result.HasValue());
+  ExpectSameOperators(moved_result.Value(), reference.Value());
+  fesa::DynamicAcceptedFrameState frame{};
+  frame.generalized_coordinates = q;
+  frame.velocity = v;
+  frame.acceleration = zero;
+  frame.algorithmic_acceleration = zero;
+  frame.external_force = force;
+  const auto parameters = fesa::MakeHhtParameters(-0.05);
+  ASSERT_TRUE(parameters.HasValue());
+  fesa::TbbParallelFor parallel;
+  const auto tbb = fesa::hht_internal::EvaluateAssembly(
+      *fixture.domain, fixture.dofs, shells.Value(), frame, zero, 0.1,
+      parameters.Value(), parallel, true);
+  ASSERT_TRUE(tbb.HasValue());
+  const auto tbb_result = evaluate(tbb.Value());
+  ASSERT_TRUE(tbb_result.HasValue());
+  ExpectSameOperators(tbb_result.Value(), reference.Value());
+}
+
+TEST(HhtReducedOperators,
+     ChangedTrialConstraintContextFallsBackInOriginalPosition) {
+  RestoreTrialEnvironment restore;
+  const auto fixture = MakeFixture(false);
+  const auto shells =
+      fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  const auto assembled =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  ASSERT_TRUE(assembled.HasValue());
+  ASSERT_TRUE(assembled.Value().trial_constraints.has_value());
+  ASSERT_TRUE(
+      assembled.Value().trial_constraints->shells[0U].directional_internal);
+  const auto chart = fesa::hht_reduced_internal::BuildFixedChart(
+      fixture.dofs, assembled.Value().rotations);
+  ASSERT_TRUE(chart.HasValue());
+  const auto evaluate = [&](const fesa::hht_internal::DynamicAssembly& value) {
+    return fesa::hht_reduced_internal::EvaluateConstraintOperators(
+        *fixture.domain, fixture.dofs, shells.Value(), value, chart.Value(),
+        zero, zero);
+  };
+  using Assembly = fesa::hht_internal::DynamicAssembly;
+  std::vector<std::function<void(Assembly&)>> changes{
+      [](Assembly& a) { a.frame.generalized_coordinates[2U] = 0.001; },
+      [](Assembly& a) { a.frame.velocity[3U] = 0.003; },
+      [](Assembly& a) { a.frame.time = 0.1; },
+      [](Assembly& a) { ++a.frame.frame_index; },
+      [](Assembly& a) { ++a.frame.increment_index; },
+      [](Assembly& a) { ++a.frame.iteration_count; },
+      [](Assembly& a) { a.frame.is_initial = true; },
+      [](Assembly& a) { a.trial_constraints->source_path = "different.inp"; },
+      [](Assembly& a) { ++a.trial_constraints->shells[0U].element_index; },
+      [](Assembly& a) { ++a.trial_constraints->shells[0U].scatter[0U]; },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].directional_internal = false;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].node_source_ids[0U].instance_name +=
+            "x";
+      },
+      [](Assembly& a) {
+        ++a.trial_constraints->shells[0U].node_source_ids[0U].source_label;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].node_source_ids[0U].source_label_text +=
+            "x";
+      },
+      [](Assembly& a) {
+        ++a.trial_constraints->shells[0U].configuration.source.line;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].configuration.source.file =
+            "different.inp";
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].configuration.identity += "x";
+      },
+      [](Assembly& a) { a.trial_constraints->shells[0U].q[0U] = -0.0; },
+      [](Assembly& a) { a.trial_constraints->shells[0U].v[0U] = -0.0; },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].drilling_constraint_gradient =
+            fesa::Matrix{0U, 0U};
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].drilling_constraint_hessian[2U] =
+            fesa::Matrix{0U, 0U};
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U]
+            .drilling_constraint_third_directional_gradient[1U] =
+            fesa::Vector{0U};
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].drilling_constraint_hessian[0U](
+            0U, 0U) = std::numeric_limits<double>::infinity();
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U]
+            .evaluation_context.dense_policy.blas_max_threads = 2;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U]
+            .evaluation_context.dense_policy.dynamic = 1;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U]
+            .evaluation_context.dense_policy.dynamic = -1;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].evaluation_context.mxcsr_control ^=
+            0x8000U;
+      },
+      [](Assembly& a) {
+        a.trial_constraints->shells[0U].evaluation_context.supported = false;
+      }};
+  for (std::size_t i = 0U; i < 64U; ++i)
+    changes.push_back([i](Assembly& a) {
+      a.trial_constraints->shells[0U].configuration.values[i] += 0.125;
+    });
+  for (const double invalid : {std::numeric_limits<double>::quiet_NaN(),
+                                std::numeric_limits<double>::infinity(),
+                                -std::numeric_limits<double>::infinity()})
+    changes.push_back([invalid](Assembly& a) {
+      a.trial_constraints->shells[0U].configuration.values[60U] = invalid;
+    });
+  for (std::size_t i = 0U; i < changes.size(); ++i) {
+    SCOPED_TRACE(i);
+    auto candidate = assembled.Value();
+    candidate.trial_constraints->shells[0U].drilling_constraint[0U] = 17.0;
+    changes[i](candidate);
+    auto fresh = candidate;
+    fresh.trial_constraints.reset();
+    const auto reference = evaluate(fresh);
+    const auto actual = evaluate(candidate);
+    ASSERT_TRUE(reference.HasValue());
+    ASSERT_TRUE(actual.HasValue());
+    ExpectSameOperators(actual.Value(), reference.Value());
+  }
+  for (std::size_t kind = 0U; kind < 2U; ++kind) {
+    auto candidate = assembled.Value();
+    if (kind == 0U) candidate.trial_constraints->shells.clear();
+    if (kind == 1U) candidate.trial_constraints.reset();
+    auto fresh = assembled.Value();
+    fresh.trial_constraints.reset();
+    const auto reference = evaluate(fresh);
+    const auto actual = evaluate(candidate);
+    ASSERT_TRUE(reference.HasValue());
+    ASSERT_TRUE(actual.HasValue());
+    ExpectSameOperators(actual.Value(), reference.Value());
+  }
+  // Change the actual consumer environment, preserving the original fallback.
+  for (std::size_t mode = 0U; mode < 5U; ++mode) {
+    SCOPED_TRACE(mode);
+    auto candidate = assembled.Value();
+    candidate.trial_constraints->shells[0U].drilling_constraint[0U] = 17.0;
+    auto fresh = candidate;
+    fresh.trial_constraints.reset();
+#if defined(_M_X64)
+    const auto csr = _mm_getcsr();
+    if (mode == 0U) _mm_setcsr(csr ^ 0x8000U);
+    if (mode == 1U) _mm_setcsr(csr ^ 0x0040U);
+#endif
+    const int rounding = std::fegetround();
+    if (mode == 2U) static_cast<void>(std::fesetround(FE_DOWNWARD));
+    if (mode == 3U) static_cast<void>(mkl_set_num_threads_local(2));
+    if (mode == 4U) mkl_set_dynamic(1);
+    const auto reference = evaluate(fresh);
+    const auto actual = evaluate(candidate);
+    if (mode == 4U) mkl_set_dynamic(0);
+    if (mode == 3U) EXPECT_EQ(mkl_set_num_threads_local(0), 2);
+    static_cast<void>(std::fesetround(rounding));
+#if defined(_M_X64)
+    _mm_setcsr(csr);
+#endif
+    ASSERT_TRUE(reference.HasValue());
+    ASSERT_TRUE(actual.HasValue());
+    ExpectSameOperators(actual.Value(), reference.Value());
+  }
+  auto failed_trial = assembled.Value();
+  failed_trial.frame.generalized_coordinates[0U] =
+      std::numeric_limits<double>::infinity();
+  auto fresh_failure = failed_trial;
+  fresh_failure.trial_constraints.reset();
+  const auto failed = evaluate(failed_trial);
+  const auto original = evaluate(fresh_failure);
+  ASSERT_FALSE(failed.HasValue());
+  ASSERT_FALSE(original.HasValue());
+  ASSERT_EQ(failed.GetStatus().Diagnostics().size(),
+            original.GetStatus().Diagnostics().size());
+  EXPECT_EQ(failed.GetStatus().Diagnostics()[0U].code,
+            original.GetStatus().Diagnostics()[0U].code);
+  EXPECT_EQ(failed.GetStatus().Diagnostics()[0U].message,
+            original.GetStatus().Diagnostics()[0U].message);
+  const auto rerun = evaluate(assembled.Value());
+  auto fresh = assembled.Value();
+  fresh.trial_constraints.reset();
+  const auto reference = evaluate(fresh);
+  ASSERT_TRUE(rerun.HasValue());
+  ASSERT_TRUE(reference.HasValue());
+  ExpectSameOperators(rerun.Value(), reference.Value());
+}
+
+TEST(HhtReducedOperators, TrialContextControlsAndInitialFalseFailureLifetime) {
+  RestoreTrialEnvironment restore;
+  const auto context = fesa::hht_internal::ReadTrialEvaluationContext();
+  ASSERT_TRUE(fesa::hht_internal::ReusableTrialContext(context, context));
+  for (const unsigned int bit : {0x0040U, 0x8000U, 0x2000U, 0x0080U}) {
+    auto changed = context;
+    changed.mxcsr_control ^= bit;
+    EXPECT_FALSE(fesa::hht_internal::ReusableTrialContext(context, changed));
+  }
+  auto unmasked = context;
+  unmasked.mxcsr_control &= ~0x0080U;
+  EXPECT_FALSE(fesa::hht_internal::ReusableTrialContext(unmasked, unmasked));
+  static_cast<void>(mkl_set_num_threads_local(2));
+  const auto positive = fesa::hht_internal::ReadTrialEvaluationContext();
+  EXPECT_FALSE(fesa::hht_internal::ReusableTrialContext(positive, positive));
+  EXPECT_EQ(mkl_set_num_threads_local(0), 2);
+  mkl_set_dynamic(1);
+  const auto dynamic = fesa::hht_internal::ReadTrialEvaluationContext();
+  EXPECT_FALSE(fesa::hht_internal::ReusableTrialContext(dynamic, dynamic));
+  mkl_set_dynamic(0);
+  const auto fixture = MakeFixture(false);
+  const auto shells =
+      fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  const auto position =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, false);
+  ASSERT_TRUE(position.HasValue());
+  EXPECT_FALSE(position.Value().trial_constraints.has_value());
+  mkl_set_dynamic(1);
+  const auto dynamic_trial =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  mkl_set_dynamic(0);
+  ASSERT_TRUE(dynamic_trial.HasValue());
+  ASSERT_TRUE(dynamic_trial.Value().trial_constraints.has_value());
+  EXPECT_FALSE(
+      dynamic_trial.Value().trial_constraints->shells[0U].directional_internal);
+  static_cast<void>(mkl_set_num_threads_local(2));
+  const auto positive_trial =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 2);
+  ASSERT_TRUE(positive_trial.HasValue());
+  ASSERT_TRUE(positive_trial.Value().trial_constraints.has_value());
+  EXPECT_FALSE(positive_trial.Value()
+                   .trial_constraints->shells[0U]
+                   .directional_internal);
+  fesa::DynamicAcceptedFrameState frame{};
+  frame.is_initial = true;
+  for (auto* vector :
+       {&frame.generalized_coordinates, &frame.velocity, &frame.acceleration,
+        &frame.algorithmic_acceleration, &frame.external_force})
+    *vector = zero;
+  fesa::SerialParallelFor parallel;
+  const auto parameters = fesa::MakeHhtParameters(-0.05);
+  ASSERT_TRUE(parameters.HasValue());
+  const auto initial = fesa::hht_internal::EvaluateAssembly(
+      *fixture.domain, fixture.dofs, shells.Value(), frame, zero, 0.1,
+      parameters.Value(), parallel, true);
+  ASSERT_TRUE(initial.HasValue());
+  EXPECT_FALSE(initial.Value().trial_constraints.has_value());
+  auto invalid = zero;
+  invalid[0U] = std::numeric_limits<double>::quiet_NaN();
+  const auto failure =
+      Assemble(fixture, shells.Value(), invalid, zero, zero, zero, true);
+  EXPECT_FALSE(failure.HasValue());
+  const auto current =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  const auto next =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  ASSERT_TRUE(current.HasValue());
+  ASSERT_TRUE(next.HasValue());
+  ASSERT_EQ(current.Value().trial_constraints->shells.size(),
+            shells.Value().size());
+  ASSERT_EQ(next.Value().trial_constraints->shells.size(),
+            shells.Value().size());
+  EXPECT_NE(&current.Value().trial_constraints->shells[0U],
+            &next.Value().trial_constraints->shells[0U]);
+}
+
+TEST(HhtReducedOperators, NonfiniteMatchingTrialKeysPreserveOriginalFailure) {
+  RestoreTrialEnvironment restore;
+  const auto fixture = MakeFixture(false);
+  const auto shells =
+      fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  const auto assembled =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  ASSERT_TRUE(assembled.HasValue());
+  ASSERT_TRUE(assembled.Value().trial_constraints.has_value());
+  ASSERT_TRUE(
+      assembled.Value().trial_constraints->shells[0U].directional_internal);
+  const auto chart = fesa::hht_reduced_internal::BuildFixedChart(
+      fixture.dofs, assembled.Value().rotations);
+  ASSERT_TRUE(chart.HasValue());
+  const auto evaluate = [&](const fesa::hht_internal::DynamicAssembly& value) {
+    return fesa::hht_reduced_internal::EvaluateConstraintOperators(
+        *fixture.domain, fixture.dofs, shells.Value(), value, chart.Value(),
+        zero, zero);
+  };
+  const std::array<double, 3> invalid{std::numeric_limits<double>::quiet_NaN(),
+                                      std::numeric_limits<double>::infinity(),
+                                      -std::numeric_limits<double>::infinity()};
+  for (std::size_t channel = 0U; channel < 2U; ++channel)
+    for (std::size_t i = 0U; i < invalid.size(); ++i) {
+      SCOPED_TRACE(channel);
+      SCOPED_TRACE(i);
+      auto candidate = assembled.Value();
+      auto& saved = candidate.trial_constraints->shells[0U];
+      if (channel == 0U) {
+        candidate.frame.generalized_coordinates[0U] = invalid[i];
+        saved.q[0U] = invalid[i];
+      } else {
+        candidate.frame.velocity[0U] = invalid[i];
+        saved.v[0U] = invalid[i];
+      }
+      auto fresh = candidate;
+      fresh.trial_constraints.reset();
+      const auto original = evaluate(fresh);
+      const auto actual = evaluate(candidate);
+      ASSERT_FALSE(original.HasValue());
+      ASSERT_FALSE(actual.HasValue());
+      EXPECT_EQ(actual.GetStatus().Category(), original.GetStatus().Category());
+      const auto& a = actual.GetStatus().Diagnostics();
+      const auto& b = original.GetStatus().Diagnostics();
+      ASSERT_EQ(a.size(), b.size());
+      for (std::size_t j = 0U; j < a.size(); ++j) {
+        EXPECT_EQ(a[j].severity, b[j].severity);
+        EXPECT_EQ(a[j].code, b[j].code);
+        EXPECT_EQ(a[j].location.file, b[j].location.file);
+        EXPECT_EQ(a[j].location.line, b[j].location.line);
+        EXPECT_EQ(a[j].keyword, b[j].keyword);
+        EXPECT_EQ(a[j].entity_identity, b[j].entity_identity);
+        EXPECT_EQ(a[j].message, b[j].message);
+      }
+    }
+}
+
+TEST(HhtReducedOperators, NonfiniteMatchingTimeKeyMissPreservesFreshSuccess) {
+  RestoreTrialEnvironment restore;
+  const auto fixture = MakeFixture(false);
+  const auto shells =
+      fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  const auto assembled =
+      Assemble(fixture, shells.Value(), zero, zero, zero, zero, true);
+  ASSERT_TRUE(assembled.HasValue());
+  ASSERT_TRUE(assembled.Value().trial_constraints.has_value());
+  const auto chart = fesa::hht_reduced_internal::BuildFixedChart(
+      fixture.dofs, assembled.Value().rotations);
+  ASSERT_TRUE(chart.HasValue());
+  const auto evaluate = [&](const fesa::hht_internal::DynamicAssembly& value) {
+    return fesa::hht_reduced_internal::EvaluateConstraintOperators(
+        *fixture.domain, fixture.dofs, shells.Value(), value, chart.Value(),
+        zero, zero);
+  };
+  const auto parameters = fesa::MakeHhtParameters(-0.05);
+  ASSERT_TRUE(parameters.HasValue());
+  fesa::SerialParallelFor parallel;
+  for (const double bad_time : {std::numeric_limits<double>::quiet_NaN(),
+                                std::numeric_limits<double>::infinity(),
+                                -std::numeric_limits<double>::infinity()}) {
+    auto candidate = assembled.Value();
+    candidate.frame.time = bad_time;
+    candidate.trial_constraints->time = bad_time;
+    candidate.trial_constraints->shells[0U].drilling_constraint[0U] = 17.0;
+    auto fresh = candidate;
+    fresh.trial_constraints.reset();
+    const auto reference = evaluate(fresh);
+    const auto actual = evaluate(candidate);
+    ASSERT_TRUE(reference.HasValue());
+    ASSERT_TRUE(actual.HasValue());
+    ExpectSameOperators(actual.Value(), reference.Value());
+    const auto producer = fesa::hht_internal::EvaluateAssembly(
+        *fixture.domain, fixture.dofs, shells.Value(), candidate.frame, zero,
+        0.1, parameters.Value(), parallel, true);
+    ASSERT_TRUE(producer.HasValue());
+    EXPECT_FALSE(producer.Value().trial_constraints.has_value());
+    const auto producer_result = evaluate(producer.Value());
+    ASSERT_TRUE(producer_result.HasValue());
+    ExpectSameOperators(producer_result.Value(), reference.Value());
+  }
+}
+
 TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) {
   const auto fixture = MakeFixture(false);
   const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain,
@@ -1141,4 +1682,58 @@ TEST(HhtReducedChannelNewton, RestrictedSingularityIsNotFullInvertibility) {
   EXPECT_FALSE(solver.FactorizeFresh(restricted.Value().tangent).IsOk());
 }
 
+TEST(HhtReducedOperators, NullMetricKeepsFiniteLeftAssociatedProducts) {
+  const auto fixture = MakeFixture(false);
+  const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain,
+                                                           fixture.dofs);
+  ASSERT_TRUE(shells.HasValue());
+  const fesa::Vector zero(fixture.dofs.FullDofCount());
+  const auto assembled = Assemble(fixture, shells.Value(), zero, zero, zero,
+                                  zero, false);
+  ASSERT_TRUE(assembled.HasValue());
+  auto assembly = assembled.Value();
+  const double maximum = (std::numeric_limits<double>::max)();
+  const double epsilon = std::ldexp(1.0, -26);
+  const double t = std::sqrt(1.0 - epsilon * epsilon);
+  const std::array<double, 3U> n{epsilon, 0.6 * t, 0.8 * t};
+  const std::array<double, 9U> block{0.0, maximum, maximum,
+                                    maximum, 0.0, 0.0,
+                                    maximum, 0.0, 0.0};
+  double original_metric = 0.0;
+  for (std::size_t i = 0U; i < 3U; ++i)
+    for (std::size_t j = 0U; j < 3U; ++j)
+      original_metric += (n[i] * block[3U * i + j]) * n[j];
+  ASSERT_TRUE(std::isfinite(original_metric));
+  ASSERT_GT(original_metric, 0.0);
+  ASSERT_FALSE(std::isfinite(block[1U] * n[1U] + block[2U] * n[2U]));
+  std::vector<fesa::CooContribution> entries;
+  for (std::size_t node = 0U; node < assembly.rotations.size(); ++node) {
+    assembly.rotations[node].drill_null = fesa::Vector3(n);
+    for (std::size_t i = 0U; i < 3U; ++i)
+      for (std::size_t j = 0U; j < 3U; ++j)
+        entries.push_back({6U * node + 3U + i, 6U * node + 3U + j,
+                           block[3U * i + j], node, 3U * i + j});
+  }
+  const auto metric = fesa::hht_internal::MakeSparse(
+      fixture.dofs.FullDofCount(), fixture.dofs.FullDofCount(),
+      std::move(entries));
+  ASSERT_TRUE(metric.HasValue());
+  assembly.drilling_metric = metric.Value();
+  const auto chart = fesa::hht_reduced_internal::BuildFixedChart(
+      fixture.dofs, assembly.rotations);
+  ASSERT_TRUE(chart.HasValue());
+  fesa::hht_reduced_internal::ConstraintOperators operators;
+  operators.rows.resize(chart.Value().null_nodes.size());
+  for (std::size_t row = 0U; row < operators.rows.size(); ++row)
+    for (std::size_t component = 0U; component < 3U; ++component)
+      operators.rows[row].gradient.push_back(
+          {6U * chart.Value().null_nodes[row] + 3U + component,
+           original_metric * n[component]});
+  fesa::MklPardisoTangentSolver solver;
+  const auto checked = fesa::hht_reduced_internal::CheckGraphRegularity(
+      assembly, chart.Value(), operators, solver);
+  ASSERT_TRUE(checked.HasValue());
+  EXPECT_GE(checked.Value(), std::ldexp(1.0, -32));
+}
+
 }  // namespace

tests/unit/assembly/load_assembler_test.cpp

@@ -481,3 +481,47 @@ TEST(LoadAssembly, ZeroLoadsRemainZero) {
   ASSERT_TRUE(constrained_rhs.HasValue());
   EXPECT_EQ(constrained_rhs.Value().Size(), 0U);
 }
+
+TEST(LoadAssembly, PreservesOrderedCancellationInPrescribedCorrection) {
+  const std::filesystem::path source{"models/load-assembly.inp"};
+  auto fixture = MakeFixture(
+      1U, {},
+      {{"10", 1, 1, 1.0, {source, 21U}},
+       {"10", 3, 3, 1.0, {source, 22U}},
+       {"10", 5, 5, 1.0, {source, 23U}}},
+      {});
+  const double large = 72057594037927936.0;  // 2^56, exactly representable.
+  const auto kfc = MakeDenseSparse(
+      3U, 3U, {large, 1.0, -large, large, -large, 1.0, 0.0, 0.0, 0.0});
+
+  const auto rhs = fesa::LoadAssembler::EffectiveFreeRhs(
+      fesa::Vector{6U}, kfc, fixture.dofs->PrescribedValues(), *fixture.dofs);
+
+  ASSERT_TRUE(rhs.HasValue());
+  ASSERT_EQ(rhs.Value().Size(), 3U);
+  EXPECT_DOUBLE_EQ(rhs.Value()[0U], 0.0);
+  EXPECT_DOUBLE_EQ(rhs.Value()[1U], -1.0);
+  EXPECT_DOUBLE_EQ(rhs.Value()[2U], 0.0);
+}
+
+TEST(LoadAssembly, SelectsFirstOrderedOverflowDiagnostic) {
+  const std::filesystem::path source{"models/load-assembly.inp"};
+  auto fixture = MakeFixture(
+      1U, {},
+      {{"10", 1, 1, 1.0, {source, 21U}},
+       {"10", 3, 3, 1.0, {source, 22U}},
+       {"10", 5, 5, 2.0, {source, 23U}}},
+      {});
+  const double maximum = (std::numeric_limits<double>::max)();
+  const auto kfc = MakeDenseSparse(
+      3U, 3U, {maximum, maximum, 0.0, 0.0, 0.0, maximum, 0.0, 0.0, 0.0});
+
+  const auto rhs = fesa::LoadAssembler::EffectiveFreeRhs(
+      fesa::Vector{6U}, kfc, fixture.dofs->PrescribedValues(), *fixture.dofs);
+
+  ExpectFailureCode(rhs, "nonfinite-load-accumulation");
+  ASSERT_EQ(rhs.GetStatus().Diagnostics().size(), 1U);
+  EXPECT_EQ(rhs.GetStatus().Diagnostics()[0U].entity_identity, "0");
+  EXPECT_NE(rhs.GetStatus().Diagnostics()[0U].message.find("row sum"),
+            std::string::npos);
+}

tests/unit/elements/finite_rotation_test.cpp

@@ -196,6 +196,21 @@ TEST(FiniteRotation,
     const Matrix left_jacobian = LeftJacobianSO3(phi);
     EXPECT_TRUE(IsFinite(left_jacobian));
 
+    // Preserve the original complete inverse-action check through the solve API.
+    for (std::size_t column = 0; column < kRotationDimension; ++column) {
+      const Vector3 basis{column == 0 ? 1.0 : 0.0,
+                          column == 1 ? 1.0 : 0.0,
+                          column == 2 ? 1.0 : 0.0};
+      const auto solved = SpatialMomentSO3(phi, basis);
+      ASSERT_TRUE(solved.HasValue());
+      EXPECT_TRUE(std::isfinite(solved.Value().X()));
+      EXPECT_TRUE(std::isfinite(solved.Value().Y()));
+      EXPECT_TRUE(std::isfinite(solved.Value().Z()));
+      const auto identity_column = Multiply(Transpose(left_jacobian), solved.Value());
+      for (std::size_t row = 0; row < kRotationDimension; ++row)
+        EXPECT_NEAR(identity_column[row], row == column ? 1.0 : 0.0, 1.0e-11);
+    }
+
     const Vector3 load{0.7, -0.2, 0.4};
     const auto spatial = SpatialMomentSO3(phi, load);
     ASSERT_TRUE(spatial.HasValue());

tests/unit/elements/hex8_element_test.cpp

@@ -1180,6 +1180,35 @@ TEST(Hex8Element, MeanVolumeWarpedLowOnlyStateRetainsWeightedCubicEnergy) {
   EXPECT_NEAR(actual_cubic.high, cubic * c * c * c, 1e-11 * c * c * c);
   EXPECT_NE(actual_cubic.high, 0);
 }
+TEST(Hex8Element, BackendOperandRejectionRetainsOriginalStiffnessDiagnostic) {
+  Coordinates coordinates{};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i)
+      coordinates[a][i] = 0.01 * kSigns[a][i];
+  // Original native probe: finite moduli/volume, but mu*gradient overflows.
+  // The failed diagnostic candidate must not become a scalar success fallback.
+  const auto element = Make(coordinates, 1e308, -0.5);
+  ASSERT_FALSE(element.HasValue());
+  EXPECT_EQ(element.GetStatus().Diagnostics()[0].code, "invalid-hex8");
+  EXPECT_EQ(element.GetStatus().Diagnostics()[0].message,
+            "Hex8 reference stiffness is nonfinite.");
+}
+
+TEST(Hex8Element, TensorFailureRetainsPriorityOverEnergyAndTangentOverflow) {
+  Coordinates coordinates{};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i)
+      coordinates[a][i] = 0.01 * kSigns[a][i];
+  const auto element = Make(coordinates, 1e200, -0.5);
+  ASSERT_TRUE(element.HasValue());
+  const auto displacement = Affine(coordinates,
+      Tensor{{{1e100, 0, 0}, {0, 1e100, 0}, {0, 0, 1e100}}});
+  const auto response = element.Value().EvaluateNonlinear(displacement);
+  ASSERT_FALSE(response.HasValue());
+  EXPECT_EQ(response.GetStatus().Diagnostics()[0].code, "invalid-hex8");
+  EXPECT_EQ(response.GetStatus().Diagnostics()[0].message,
+            "Hex8 trial tensor is nonfinite.");
+}
 }  // namespace
 
 TEST(Hex8Element, RejectsNonfiniteCanonicalJacobianActionBeforeFiniteGradients) {

tests/unit/elements/hex8_high_low_test.cpp

@@ -2,16 +2,836 @@
 
 #include <gtest/gtest.h>
 
+#include <array>
+#include <atomic>
+#include <cfenv>
 #include <cmath>
+#include <cstdint>
+#include <cstring>
+#include <future>
+#include <limits>
+#include <random>
+#include <thread>
+#include <type_traits>
+#include <utility>
+#include <vector>
 
 #include "fesa/materials/isotropic_linear_elastic_material.h"
 #include "../../../src/fesa/math/high_low_arithmetic.h"
 #include "../../../src/fesa/elements/hex8_geometry.h"
+#include "../../../src/fesa/elements/hex8_pair_simd.h"
+#include "../math/high_low_arithmetic_original.h"
 
 namespace {
 using fesa::math_internal::HighLow;
 using fesa::math_internal::Subtract;
 
+std::uint64_t PairBits(double value) {
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  return bits;
+}
+
+void ExpectOriginalProduct(HighLow left, HighLow right) {
+  const auto actual = fesa::math_internal::Multiply(left, right);
+  const auto expected = fesa::high_low_original::Multiply(
+      {left.high, left.low}, {right.high, right.low});
+  EXPECT_EQ(PairBits(actual.high), PairBits(expected.high));
+  EXPECT_EQ(PairBits(actual.low), PairBits(expected.low));
+  const auto actual_sum = fesa::math_internal::Add(left, right);
+  const auto expected_sum = fesa::high_low_original::Add(
+      {left.high, left.low}, {right.high, right.low});
+  EXPECT_EQ(PairBits(actual_sum.high), PairBits(expected_sum.high));
+  EXPECT_EQ(PairBits(actual_sum.low), PairBits(expected_sum.low));
+}
+
+TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) {
+  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
+  const double minimum = std::numeric_limits<double>::denorm_min();
+  const double maximum = std::numeric_limits<double>::max();
+  const double infinity = std::numeric_limits<double>::infinity();
+  const double nan = std::numeric_limits<double>::quiet_NaN();
+  std::vector<HighLow> pairs{{}, {-0.0, 0.0}, {0.0, -0.0}, {-0.0, -0.0},
+      {1, -0.0}, {-1, -0.0}, {1, 1}, {0, minimum}, {infinity, 0},
+      {0, infinity}, {nan, 0}, {0, nan}, {maximum}, {-maximum},
+      {minimum}, {-minimum}, {std::numeric_limits<double>::min()},
+      {0.5}, {std::nextafter(0.5, 1.0)},
+      {1, std::ldexp(1.0, -53)},
+      {std::nextafter(1.0, 2.0), -std::ldexp(1.0, -53)}};
+  for (const auto left : pairs)
+    for (const auto right : pairs) ExpectOriginalProduct(left, right);
+  for (int exponent = -1074; exponent <= 1023; ++exponent) {
+    const double value = std::ldexp(1.0, exponent);
+    for (double sign : {-1.0, 1.0})
+      for (const auto right : pairs) {
+        ExpectOriginalProduct(HighLow{sign * value}, right);
+        ExpectOriginalProduct(right, HighLow{sign * value});
+      }
+  }
+  std::mt19937_64 random(0x48e81234U);
+  for (int sample = 0; sample < 12000; ++sample) {
+    double left = 0;
+    double right = 0;
+    const auto left_bits = random();
+    const auto right_bits = random();
+    std::memcpy(&left, &left_bits, sizeof(left));
+    std::memcpy(&right, &right_bits, sizeof(right));
+    const auto pair = fesa::math_internal::TwoSum(left, right);
+    ExpectOriginalProduct(pair, HighLow{});
+    ExpectOriginalProduct(HighLow{}, pair);
+    ExpectOriginalProduct(pair, HighLow{right});
+    ExpectOriginalProduct(HighLow{right}, pair);
+    ExpectOriginalProduct(pair, pair);
+  }
+}
+
+fesa::hex8_internal::PacketPair RawPacketPair(HighLow value) {
+  return {value.high, value.low};
+}
+
+fesa::hex8_internal::PairPacket FrozenMaterialPacket(
+    const fesa::hex8_internal::MaterialPacketInput& input) {
+  using Original = fesa::high_low_original::HighLow;
+  const auto original = [](fesa::hex8_internal::PacketPair value) {
+    return Original{value.high, value.low};
+  };
+  fesa::hex8_internal::PairPacket result;
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const auto first = fesa::high_low_original::Multiply(
+        original(input.deviatoric_left.values[lane]),
+        original(input.deviatoric_right.values[lane]));
+    const auto second = fesa::high_low_original::Multiply(
+        original(input.shear_left.values[lane]), original(input.shear_right.values[lane]));
+    const auto third = fesa::high_low_original::Multiply(
+        original(input.tensor_left.values[lane]), original(input.gradient_dot.values[lane]));
+    const auto value = fesa::high_low_original::Add(
+        fesa::high_low_original::Add(first, second), third);
+    const auto accumulated = fesa::high_low_original::Add(
+        original(input.accumulated.values[lane]),
+        fesa::high_low_original::Multiply(value, original(input.weight)));
+    result.values[lane] = {accumulated.high, accumulated.low};
+  }
+  return result;
+}
+
+bool ExpectFrozenMaterialPacket(
+    const fesa::hex8_internal::MaterialPacketInput& input) {
+  const auto expected = FrozenMaterialPacket(input);
+  fesa::hex8_internal::PairPacket scalar, selected;
+  fesa::hex8_internal::EvaluateMaterialPacketScalar(input, &scalar);
+  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  evaluator(input, &selected);
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    SCOPED_TRACE(lane);
+    EXPECT_EQ(PairBits(scalar.values[lane].high), PairBits(expected.values[lane].high));
+    EXPECT_EQ(PairBits(scalar.values[lane].low), PairBits(expected.values[lane].low));
+    EXPECT_EQ(PairBits(selected.values[lane].high), PairBits(expected.values[lane].high));
+    EXPECT_EQ(PairBits(selected.values[lane].low), PairBits(expected.values[lane].low));
+  }
+  return evaluator != fesa::hex8_internal::EvaluateMaterialPacketScalar;
+}
+
+TEST(Hex8HighLow, FourColumnMaterialPacketMatchesFrozenScalarBits) {
+  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
+  using fesa::hex8_internal::MaterialPacketInput;
+  RecordProperty("private_packet_backend",
+      fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+          fesa::hex8_internal::EvaluateMaterialPacketScalar ? "scalar" : "avx2");
+  std::size_t packet_checks = 0, avx2_evaluator_checks = 0;
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const double maximum = std::numeric_limits<double>::max();
+  const std::vector<HighLow> boundary{{}, {-0.0, 0}, {0, -0.0},
+      {1}, {-1}, {1, -0.0}, {1, std::ldexp(1.0, -54)},
+      {-1, std::ldexp(1.0, -54)}, {tiny}, {-tiny},
+      {std::numeric_limits<double>::min()}, {maximum}, {-maximum},
+      {1, tiny}, {1, -tiny}, {std::nextafter(1.0, 2.0)},
+      {1, 1}, {0, tiny}, {std::numeric_limits<double>::infinity(), 0},
+      {std::numeric_limits<double>::quiet_NaN(), 0},
+      {std::numeric_limits<double>::signaling_NaN(), 0},
+      {0, std::numeric_limits<double>::signaling_NaN()}};
+  for (std::size_t first = 0; first < boundary.size(); ++first)
+    for (std::size_t second = 0; second < boundary.size(); ++second) {
+      MaterialPacketInput input{};
+      input.weight = {0.125};
+      for (std::size_t lane = 0; lane < 4; ++lane) {
+        input.deviatoric_left.values[lane] = RawPacketPair(boundary[first]);
+        input.deviatoric_right.values[lane] = RawPacketPair(boundary[second]);
+        input.shear_left.values[lane] = {2, std::ldexp(1.0, -54)};
+        input.shear_right.values[lane] = {-0.5, std::ldexp(1.0, -56)};
+        input.tensor_left.values[lane] = {0.75};
+        input.gradient_dot.values[lane] = {-0.25};
+        input.accumulated.values[lane] = lane == 0 ? fesa::hex8_internal::PacketPair{0, 0}
+            : fesa::hex8_internal::PacketPair{0.125, 0};
+      }
+      avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
+      ++packet_checks;
+    }
+  for (int exponent = -1074; exponent <= 1023; ++exponent) {
+    MaterialPacketInput input{};
+    input.weight = {1};
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      input.deviatoric_left.values[lane] = {
+          std::ldexp(lane % 2U == 0 ? 1.0 : -1.0, exponent)};
+      input.deviatoric_right.values[lane] = {lane < 2 ? 0.5 : 2.0};
+    }
+    avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
+    ++packet_checks;
+  }
+  std::mt19937_64 random(0x4c0116U);
+  for (std::size_t sample = 0; sample < 2048; ++sample) {
+    MaterialPacketInput input{};
+    const auto pair = [&]() {
+      const double high = std::ldexp(
+          1.0 + static_cast<double>(random() % 256U) / 256.0,
+          static_cast<int>(random() % 201U) - 100);
+      const double signed_high = (random() & 1U) != 0 ? high : -high;
+      const double low = sample % 3U == 0 ? 0.0 : std::ldexp(high, -55);
+      return fesa::hex8_internal::PacketPair{
+          signed_high, (random() & 1U) != 0 ? low : -low};
+    };
+    input.weight = sample % 2U == 0
+        ? fesa::hex8_internal::PacketPair{0.125, 0} : pair();
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      input.deviatoric_left.values[lane] = pair();
+      input.deviatoric_right.values[lane] = pair();
+      input.shear_left.values[lane] = pair();
+      input.shear_right.values[lane] = pair();
+      input.tensor_left.values[lane] = pair();
+      input.gradient_dot.values[lane] = pair();
+      input.accumulated.values[lane] = sample % 4U == 0
+          ? fesa::hex8_internal::PacketPair{0, 0} : pair();
+    }
+    avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
+    ++packet_checks;
+  }
+  RecordProperty("packet_oracle_checks", static_cast<int>(packet_checks));
+  RecordProperty("avx2_evaluator_packet_checks",
+                 static_cast<int>(avx2_evaluator_checks));
+}
+
+TEST(Hex8HighLow, PacketSpecialLaneDoesNotChangeOtherColumns) {
+  fesa::hex8_internal::MaterialPacketInput input{};
+  input.weight = {0.125};
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    input.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
+    input.deviatoric_right.values[lane] = {2, std::ldexp(1.0, -55)};
+    input.shear_left.values[lane] = {0.5};
+    input.shear_right.values[lane] = {-0.25};
+    input.tensor_left.values[lane] = {0.75};
+    input.gradient_dot.values[lane] = {0.125};
+  }
+  const auto reference = FrozenMaterialPacket(input);
+  const std::vector<HighLow> special{{}, {-0.0, 0}, {1},
+      {std::numeric_limits<double>::max()},
+      {std::numeric_limits<double>::denorm_min()}, {1, 1},
+      {std::numeric_limits<double>::quiet_NaN(), 0}};
+  for (std::size_t lane = 0; lane < 4; ++lane)
+    for (const auto value : special) {
+      auto changed = input;
+      changed.deviatoric_left.values[lane] = RawPacketPair(value);
+      ExpectFrozenMaterialPacket(changed);
+      fesa::hex8_internal::PairPacket actual;
+      fesa::hex8_internal::SelectMaterialPacketEvaluator()(changed, &actual);
+      for (std::size_t other = 0; other < 4; ++other)
+        if (other != lane) {
+          EXPECT_EQ(PairBits(actual.values[other].high), PairBits(reference.values[other].high));
+          EXPECT_EQ(PairBits(actual.values[other].low), PairBits(reference.values[other].low));
+        }
+    }
+}
+
+TEST(Hex8HighLow, UniformDyadicPacketsExerciseEachLowZeroShape) {
+  std::size_t checks = 0, avx2_checks = 0;
+  for (int shape = 0; shape < 4; ++shape)
+    for (int low_weight = 0; low_weight < 2; ++low_weight)
+      for (int zero_accumulator = 0; zero_accumulator < 2; ++zero_accumulator) {
+        fesa::hex8_internal::MaterialPacketInput input{};
+        input.weight = {0.5, low_weight != 0 ? std::ldexp(1.0, -57) : 0.0};
+        for (std::size_t lane = 0; lane < 4; ++lane) {
+          const bool left_low = (shape & 1) != 0;
+          const bool right_low = (shape & 2) != 0;
+          // Every lane has the same normalized nonzero dyadic inputs. Component
+          // LBE sums stay above -200 and all Add intermediates are bounded by 4.
+          // Thus each Multiply low-shape is uniform and its component products
+          // bypass zero/unsupported-tail fallback. No mixed-zero Add fallback
+          // is possible; the initial zero accumulator is uniform across lanes.
+          input.deviatoric_left.values[lane] = {1, left_low ? std::ldexp(1.0, -54) : 0};
+          input.deviatoric_right.values[lane] = {2, right_low ? std::ldexp(1.0, -54) : 0};
+          input.shear_left.values[lane] = {0.5, left_low ? std::ldexp(1.0, -56) : 0};
+          input.shear_right.values[lane] = {0.25, right_low ? std::ldexp(1.0, -57) : 0};
+          input.tensor_left.values[lane] = {0.75, left_low ? std::ldexp(1.0, -56) : 0};
+          input.gradient_dot.values[lane] = {0.125, right_low ? std::ldexp(1.0, -58) : 0};
+          input.accumulated.values[lane] = zero_accumulator != 0
+              ? fesa::hex8_internal::PacketPair{0, 0}
+              : fesa::hex8_internal::PacketPair{0.0625, 0};
+        }
+        avx2_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
+        ++checks;
+      }
+  RecordProperty("uniform_low_shape_checks", static_cast<int>(checks));
+  RecordProperty("avx2_uniform_low_shape_checks", static_cast<int>(avx2_checks));
+}
+
+TEST(Hex8HighLow, OverflowLaneIsNeutralizedBeforeResidualFma) {
+  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  if (evaluator == fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "AVX2 masked environment is not selected; routing unverified.";
+  std::fenv_t original;
+  ASSERT_EQ(std::fegetenv(&original), 0);
+  for (std::size_t bad = 0; bad < 4; ++bad) {
+    fesa::hex8_internal::MaterialPacketInput input{};
+    input.weight = {0.5};
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      input.deviatoric_left.values[lane] = {
+          lane == bad ? std::numeric_limits<double>::max() : 1};
+      input.deviatoric_right.values[lane] = {2};
+      input.shear_left.values[lane] = {0.5};
+      input.shear_right.values[lane] = {0.25};
+      input.tensor_left.values[lane] = {0.75};
+      input.gradient_dot.values[lane] = {0.125};
+      input.accumulated.values[lane] = {0.0625};
+    }
+    // All operands are finite/nonzero and max*2 has LBE sum 972, so Product
+    // reaches its vector multiply; only the post-product finite mask rejects it.
+    const auto expected = FrozenMaterialPacket(input);
+    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
+    fesa::hex8_internal::PairPacket output;
+    evaluator(input, &output);
+    const bool extra_invalid = (std::fetestexcept(FE_INVALID) & FE_INVALID) != 0;
+    EXPECT_EQ(std::fesetenv(&original), 0);
+    EXPECT_FALSE(extra_invalid);
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      EXPECT_EQ(PairBits(output.values[lane].high), PairBits(expected.values[lane].high));
+      EXPECT_EQ(PairBits(output.values[lane].low), PairBits(expected.values[lane].low));
+    }
+  }
+  RecordProperty("avx2_product_overflow_positions", 4);
+}
+
+#if defined(_M_X64) || defined(__SSE2__)
+// Runs scalar and selected packet from identical FP state. The frozen original
+// supplies raw-bit results; the unchanged baseline supplies FE/MXCSR status.
+void ExpectPacketFlagsAndBits(
+    const fesa::hex8_internal::MaterialPacketInput& input, bool prior_status) {
+  std::fenv_t original;
+  ASSERT_EQ(std::fegetenv(&original), 0);
+  const auto original_csr = _mm_getcsr();
+  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  struct Observation {
+    fesa::hex8_internal::PairPacket output;
+    int flags;
+    unsigned int status;
+  };
+  const auto observe = [&](fesa::hex8_internal::MaterialPacketEvaluator backend) {
+    EXPECT_EQ(std::fesetenv(&original), 0);
+    _mm_setcsr(original_csr & ~0x3fU);
+    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
+    if (prior_status) {
+      EXPECT_EQ(std::feraiseexcept(FE_INEXACT), 0);
+      _mm_setcsr(_mm_getcsr() | 0x02U);  // Preserve prior denormal status too.
+    }
+    Observation value{};
+    backend(input, &value.output);
+    value.flags = std::fetestexcept(FE_ALL_EXCEPT);
+    value.status = _mm_getcsr() & 0x3fU;
+    return value;
+  };
+  const auto scalar = observe(fesa::hex8_internal::EvaluateMaterialPacketScalar);
+  const auto selected = observe(evaluator);
+  const auto expected = FrozenMaterialPacket(input);
+  // Restore both x87 and MXCSR before assertions, including NaN bit assertions.
+  EXPECT_EQ(std::fesetenv(&original), 0);
+  _mm_setcsr(original_csr);
+  EXPECT_EQ(selected.flags, scalar.flags);
+  EXPECT_EQ(selected.status, scalar.status);
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    SCOPED_TRACE(lane);
+    EXPECT_EQ(PairBits(selected.output.values[lane].high),
+              PairBits(expected.values[lane].high));
+    EXPECT_EQ(PairBits(selected.output.values[lane].low),
+              PairBits(expected.values[lane].low));
+    EXPECT_EQ(PairBits(scalar.output.values[lane].high),
+              PairBits(expected.values[lane].high));
+    EXPECT_EQ(PairBits(scalar.output.values[lane].low),
+              PairBits(expected.values[lane].low));
+  }
+}
+
+fesa::hex8_internal::MaterialPacketInput DyadicPacketFixture() {
+  fesa::hex8_internal::MaterialPacketInput input{};
+  input.weight = {0.5, 0};
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    input.deviatoric_left.values[lane] = {1, 0};
+    input.deviatoric_right.values[lane] = {2, 0};
+    input.shear_left.values[lane] = {0.5, 0};
+    input.shear_right.values[lane] = {0.25, 0};
+    input.tensor_left.values[lane] = {0.75, 0};
+    input.gradient_dot.values[lane] = {0.125, 0};
+  }
+  return input;
+}
+
+TEST(Hex8HighLow, PacketVectorGuardsPreserveScalarFlagsAndBits) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "AVX2 masked environment is not selected.";
+  using fesa::hex8_internal::PacketPair;
+  const double maximum = std::numeric_limits<double>::max();
+  const double snan = std::numeric_limits<double>::signaling_NaN();
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const std::vector<PacketPair> invalid_left{
+      {-0.0, 0}, {0, -0.0}, {1, 1}, {snan, 0}, {0, snan}};
+  const std::vector<PacketPair> dangerous_right{
+      {maximum, maximum}, {snan, 0}, {0, snan}};
+  std::size_t checks = 0;
+  for (const auto left : invalid_left)
+    for (const auto right : dangerous_right)
+      for (std::size_t bad = 0; bad < 4; ++bad) {
+        auto input = DyadicPacketFixture();
+        input.deviatoric_left.values[bad] = left;
+        input.deviatoric_right.values[bad] = right;
+        // Later lanes must still contribute the unchanged bridge's RN status.
+        const auto later = (bad + 1) % 4;
+        input.deviatoric_left.values[later] = {1, std::ldexp(1.0, -54)};
+        input.deviatoric_right.values[later] = {tiny, 0};
+        ExpectPacketFlagsAndBits(input, false);
+        ExpectPacketFlagsAndBits(input, true);
+        checks += 2;
+      }
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    auto mixed = DyadicPacketFixture();
+    mixed.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
+    mixed.deviatoric_left.values[(lane + 1) % 4] = {maximum, maximum};
+    ExpectPacketFlagsAndBits(mixed, false);
+    ExpectPacketFlagsAndBits(mixed, true);
+    auto signed_low = DyadicPacketFixture();
+    signed_low.deviatoric_left.values[lane].low = -0.0;
+    signed_low.accumulated.values[lane] = {-0.0, 0};
+    ExpectPacketFlagsAndBits(signed_low, false);
+    checks += 3;
+  }
+  RecordProperty("avx2_guard_flags_and_bit_checks", static_cast<int>(checks));
+}
+
+TEST(Hex8HighLow, PacketTailScreenRetainsExactBoundaryAndFallback) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "AVX2 masked environment is not selected.";
+  const auto odd = [](int exponent) {
+    return std::nextafter(std::ldexp(1.0, exponent),
+                          std::numeric_limits<double>::infinity());
+  };
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const std::array<std::array<double, 2>, 8> operands{{
+      {odd(-500), odd(-470)},      // screen true, exact LBE sum -1074
+      {odd(-500), odd(-471)},      // screen false, exact rejection -1075
+      {std::ldexp(1.0, -500), std::ldexp(1.0, -471)}, // false but supported
+      {tiny, odd(52)}, {tiny, odd(51)}, {3 * tiny, odd(52)},
+      {std::numeric_limits<double>::max(), 2}, {0, odd(52)}}};
+  std::size_t checks = 0;
+  for (const auto& pair : operands)
+    for (std::size_t special = 0; special < 4; ++special) {
+      fesa::hex8_internal::MaterialPacketInput input{};
+      input.weight = {1, 0};
+      for (std::size_t lane = 0; lane < 4; ++lane) {
+        const double sign = lane % 2 == 0 ? 1.0 : -1.0;
+        input.deviatoric_left.values[lane] =
+            {sign * (lane == special ? pair[0] : odd(-500)), 0};
+        input.deviatoric_right.values[lane] =
+            {lane == special ? pair[1] : odd(-470), 0};
+      }
+      ExpectPacketFlagsAndBits(input, false);
+      ExpectPacketFlagsAndBits(input, true);
+      checks += 2;
+    }
+  RecordProperty("avx2_tail_screen_flags_and_bit_checks", static_cast<int>(checks));
+}
+#endif
+
+#if defined(_M_X64) || defined(__SSE2__)
+fesa::hex8_internal::GeometryPacketOutput FrozenGeometryPacket(
+    const fesa::hex8_internal::GeometryPacketInput& input) {
+  const auto pair = [](fesa::hex8_internal::PacketPair value) {
+    return fesa::high_low_original::HighLow{value.high, value.low};
+  };
+  fesa::hex8_internal::GeometryPacketOutput output{};
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    fesa::high_low_original::HighLow dot;
+    for (std::size_t i = 0; i < 3; ++i)
+      dot = fesa::high_low_original::Add(dot,
+          fesa::high_low_original::Multiply(pair(input.gradient_a[i]),
+              pair(input.gradient_b[i].values[lane])));
+    fesa::high_low_original::HighLow geo;
+    for (std::size_t j = 0; j < 3; ++j)
+      for (std::size_t l = 0; l < 3; ++l)
+        geo = fesa::high_low_original::Add(geo,
+            fesa::high_low_original::Multiply(pair(input.gradient_stress[j][l]),
+                pair(input.gradient_b[l].values[lane])));
+    const auto weighted = fesa::high_low_original::Multiply(geo, pair(input.weight));
+    output.gradient_dot.values[lane] = {dot.high, dot.low};
+    output.weighted_geometric.values[lane] = {weighted.high, weighted.low};
+  }
+  return output;
+}
+
+fesa::hex8_internal::GeometryPacketInput GeometryDyadicFixture() {
+  fesa::hex8_internal::GeometryPacketInput input{};
+  input.weight = {0.125, 0};
+  for (std::size_t j = 0; j < 3; ++j) {
+    input.gradient_a[j] = {static_cast<double>(j + 1) * 0.25, 0};
+    for (std::size_t l = 0; l < 3; ++l) {
+      // Deliberately asymmetric first-products and distinct b/component mapping.
+      input.gradient_stress[j][l] =
+          {static_cast<double>(1 + 3 * j + l) * 0.125, 0};
+      for (std::size_t lane = 0; lane < 4; ++lane)
+        input.gradient_b[l].values[lane] =
+            {static_cast<double>(1 + 4 * l + lane) * 0.0625, 0};
+    }
+  }
+  return input;
+}
+
+// Original observation repeats geo*weight for three independent accumulators.
+// Selected observation reuses the weighted pair but retains each checked Add.
+void ExpectGeometryPacketBitsAndFlags(
+    const fesa::hex8_internal::GeometryPacketInput& input, bool prior_status) {
+  SCOPED_TRACE(prior_status);
+  std::fenv_t original;
+  ASSERT_EQ(std::fegetenv(&original), 0);
+  const auto original_csr = _mm_getcsr();
+  struct Observation {
+    fesa::hex8_internal::GeometryPacketOutput output;
+    fesa::hex8_internal::PairPacket accumulators[3];
+    int flags;
+    unsigned int status;
+  };
+  const auto current_pair = [](fesa::hex8_internal::PacketPair value) {
+    return HighLow{value.high, value.low};
+  };
+  const auto observe = [&](bool original_order) {
+    EXPECT_EQ(std::fesetenv(&original), 0);
+    _mm_setcsr(original_csr & ~0x3fU);
+    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
+    if (prior_status) {
+      EXPECT_EQ(std::feraiseexcept(FE_INEXACT), 0);
+      _mm_setcsr(_mm_getcsr() | 0x02U);
+    }
+    Observation value{};
+    if (!original_order)
+      fesa::hex8_internal::EvaluateGeometryPacketAvx2(input, &value.output);
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      HighLow geo;
+      if (original_order) {
+        HighLow dot;
+        for (std::size_t i = 0; i < 3; ++i)
+          dot = fesa::math_internal::Add(dot,
+              fesa::math_internal::Multiply(current_pair(input.gradient_a[i]),
+                  current_pair(input.gradient_b[i].values[lane])));
+        for (std::size_t j = 0; j < 3; ++j)
+          for (std::size_t l = 0; l < 3; ++l)
+            geo = fesa::math_internal::Add(geo,
+                fesa::math_internal::Multiply(
+                    current_pair(input.gradient_stress[j][l]),
+                    current_pair(input.gradient_b[l].values[lane])));
+        value.output.gradient_dot.values[lane] = {dot.high, dot.low};
+      }
+      const HighLow accumulators[3]{{},
+          {static_cast<double>(lane + 1), std::ldexp(1.0, -54)}, {-1, -0.0}};
+      for (std::size_t component = 0; component < 3; ++component) {
+        const auto weighted = original_order
+            ? fesa::math_internal::Multiply(geo, current_pair(input.weight))
+            : current_pair(value.output.weighted_geometric.values[lane]);
+        if (original_order)
+          value.output.weighted_geometric.values[lane] =
+              {weighted.high, weighted.low};
+        const auto accumulated = fesa::math_internal::Add(
+            accumulators[component], weighted);
+        value.accumulators[component].values[lane] =
+            {accumulated.high, accumulated.low};
+      }
+    }
+    value.flags = std::fetestexcept(FE_ALL_EXCEPT);
+    value.status = _mm_getcsr() & 0x3fU;
+    return value;
+  };
+  const auto baseline = observe(true);
+  const auto selected = observe(false);
+  fesa::hex8_internal::GeometryPacketOutput scalar;
+  fesa::hex8_internal::EvaluateGeometryPacketScalar(input, &scalar);
+  const auto expected = FrozenGeometryPacket(input);
+  EXPECT_EQ(std::fesetenv(&original), 0);
+  _mm_setcsr(original_csr);
+  EXPECT_EQ(selected.flags, baseline.flags);
+  EXPECT_EQ(selected.status, baseline.status);
+  const auto same_bits = [](const fesa::hex8_internal::PairPacket& a,
+                            const fesa::hex8_internal::PairPacket& b) {
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      EXPECT_EQ(PairBits(a.values[lane].high), PairBits(b.values[lane].high));
+      EXPECT_EQ(PairBits(a.values[lane].low), PairBits(b.values[lane].low));
+    }
+  };
+  same_bits(selected.output.gradient_dot, expected.gradient_dot);
+  same_bits(selected.output.weighted_geometric, expected.weighted_geometric);
+  same_bits(scalar.gradient_dot, expected.gradient_dot);
+  same_bits(scalar.weighted_geometric, expected.weighted_geometric);
+  for (std::size_t component = 0; component < 3; ++component)
+    same_bits(selected.accumulators[component], baseline.accumulators[component]);
+}
+
+TEST(Hex8HighLow, GeometryPacketPreservesAsymmetricLanesAndUniformShapes) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
+  RecordProperty("private_geometry_packet_backend", "avx2");
+  static_assert(std::is_trivial_v<fesa::hex8_internal::GeometryPacketInput>);
+  static_assert(std::is_trivial_v<fesa::hex8_internal::GeometryPacketOutput>);
+  std::size_t checks = 0;
+  for (int shape = 0; shape < 4; ++shape) {
+    auto input = GeometryDyadicFixture();
+    for (std::size_t j = 0; j < 3; ++j) {
+      if ((shape & 1) != 0) {
+        input.gradient_a[j].low = std::ldexp(input.gradient_a[j].high, -55);
+        for (std::size_t l = 0; l < 3; ++l)
+          input.gradient_stress[j][l].low =
+              std::ldexp(input.gradient_stress[j][l].high, -55);
+      }
+      if ((shape & 2) != 0)
+        for (std::size_t lane = 0; lane < 4; ++lane)
+          input.gradient_b[j].values[lane].low =
+              std::ldexp(input.gradient_b[j].values[lane].high, -55);
+    }
+    ExpectGeometryPacketBitsAndFlags(input, false);
+    ExpectGeometryPacketBitsAndFlags(input, true);
+    checks += 2;
+  }
+  RecordProperty("avx2_geometry_uniform_bit_flag_checks", static_cast<int>(checks));
+}
+
+TEST(Hex8HighLow, GeometryPacketRetainsSpecialLanesAndExactTailBoundary) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
+  using fesa::hex8_internal::PacketPair;
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const double maximum = std::numeric_limits<double>::max();
+  const std::vector<PacketPair> special{{0, 0}, {-0.0, 0}, {1, -0.0}, {1, 1},
+      {tiny, 0}, {maximum, maximum},
+      {std::numeric_limits<double>::signaling_NaN(), 0},
+      {0, std::numeric_limits<double>::signaling_NaN()}};
+  std::size_t checks = 0;
+  for (const auto pair : special)
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      auto input = GeometryDyadicFixture();
+      input.gradient_b[1].values[lane] = pair;
+      ExpectGeometryPacketBitsAndFlags(input, false);
+      ExpectGeometryPacketBitsAndFlags(input, true);
+      checks += 2;
+    }
+  const auto odd = [](int exponent) {
+    return std::nextafter(std::ldexp(1.0, exponent),
+                          std::numeric_limits<double>::infinity());
+  };
+  for (int exponent : {-470, -471})
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      auto input = GeometryDyadicFixture();
+      input.gradient_stress[0][0] = {odd(-500), 0};
+      input.gradient_b[0].values[lane] = {odd(exponent), 0};
+      ExpectGeometryPacketBitsAndFlags(input, false);
+      checks++;
+    }
+  RecordProperty("avx2_geometry_special_tail_checks", static_cast<int>(checks));
+}
+
+TEST(Hex8HighLow, GeometryPacketPreservesOrderedOverflowAndPreparedInvalid) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
+  const double maximum = std::numeric_limits<double>::max();
+  std::size_t checks = 0;
+  for (bool overflow_order : {false, true}) {
+    SCOPED_TRACE(overflow_order);
+    fesa::hex8_internal::GeometryPacketInput input{};
+    input.weight = {0.5, 0};
+    for (std::size_t l = 0; l < 3; ++l)
+      for (std::size_t lane = 0; lane < 4; ++lane)
+        input.gradient_b[l].values[lane] = {1, 0};
+    input.gradient_stress[0][0] = {maximum, 0};
+    input.gradient_stress[0][1] = {overflow_order ? maximum : -maximum, 0};
+    input.gradient_stress[0][2] = {overflow_order ? -maximum : maximum, 0};
+    ExpectGeometryPacketBitsAndFlags(input, false);
+    ExpectGeometryPacketBitsAndFlags(input, true);
+    checks += 2;
+  }
+  auto input = GeometryDyadicFixture();
+  const auto invalid_first = fesa::high_low_original::Multiply(
+      {std::ldexp(1.0, -600), 0}, {std::ldexp(1.0, -500), 0});
+  input.gradient_stress[0][0] = {invalid_first.high, invalid_first.low};
+  for (std::size_t lane = 0; lane < 4; ++lane)
+    input.gradient_b[0].values[lane] = {std::ldexp(1.0, 600), 0};
+  ExpectGeometryPacketBitsAndFlags(input, false);
+  ExpectGeometryPacketBitsAndFlags(input, true);
+  checks += 2;
+  RecordProperty("avx2_geometry_ordered_invalid_checks", static_cast<int>(checks));
+}
+
+TEST(Hex8HighLow, GeometryPacketSumFallbackPreservesNearMaxAndSubnormalStatus) {
+  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
+      fesa::hex8_internal::EvaluateMaterialPacketScalar)
+    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
+  const double maximum = std::numeric_limits<double>::max();
+  const double normal = std::ldexp(1.0, -1022);
+  const double tiny = std::ldexp(1.0, -1074);
+  const double quarter = std::ldexp(1.0, 969);
+  std::size_t checks = 0;
+  for (double tail : {0.0, normal, tiny})
+    for (std::size_t bad = 0; bad < 4; ++bad) {
+      SCOPED_TRACE(tail);
+      SCOPED_TRACE(bad);
+      fesa::hex8_internal::GeometryPacketInput input{};
+      input.weight = {0.5, 0};
+      // Both pairs are normalized. Head sum is finite; middle contributes
+      // 2^970 and the joined sum overflows. Other lanes retain finite values.
+      input.gradient_stress[0][0] = {maximum, quarter};
+      input.gradient_stress[0][1] = {quarter, tail};
+      for (std::size_t lane = 0; lane < 4; ++lane) {
+        const double scale = lane == bad ? 1.0 : 0.5;
+        input.gradient_b[0].values[lane] = {scale, 0};
+        input.gradient_b[1].values[lane] = {scale, 0};
+        input.gradient_b[2].values[lane] = {1, 0};
+      }
+      ExpectGeometryPacketBitsAndFlags(input, false);
+      ExpectGeometryPacketBitsAndFlags(input, true);
+      checks += 2;
+    }
+  fesa::hex8_internal::GeometryPacketInput finite{};
+  finite.weight = {1, 0};
+  finite.gradient_stress[0][0] = {normal, 0};
+  finite.gradient_stress[0][1] = {normal + tiny, 0};
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    finite.gradient_b[0].values[lane] = {1, 0};
+    finite.gradient_b[1].values[lane] = {1, 0};
+  }
+  // TwoSum(N,N+d) has finite normal inputs and a subnormal error d.
+  ExpectGeometryPacketBitsAndFlags(finite, false);
+  ExpectGeometryPacketBitsAndFlags(finite, true);
+  checks += 2;
+  // Mixed prefix: finite subnormal-error lanes run before another lane's head
+  // overflow. The bridge must preserve the whole packet's status union.
+  for (std::size_t bad = 0; bad < 4; ++bad) {
+    SCOPED_TRACE(bad);
+    fesa::hex8_internal::GeometryPacketInput mixed{};
+    mixed.weight = {1, 0};
+    mixed.gradient_stress[0][2] = {normal, 0};
+    mixed.gradient_stress[1][2] = {normal + tiny, 0};
+    mixed.gradient_stress[2][0] = {maximum, 0};
+    mixed.gradient_stress[2][1] = {maximum, 0};
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      mixed.gradient_b[0].values[lane] = {lane == bad ? 1.0 : 0.0, 0};
+      mixed.gradient_b[1].values[lane] = {lane == bad ? 1.0 : 0.0, 0};
+      mixed.gradient_b[2].values[lane] = {lane == bad ? 0.0 : 1.0, 0};
+    }
+    ExpectGeometryPacketBitsAndFlags(mixed, false);
+    ExpectGeometryPacketBitsAndFlags(mixed, true);
+    checks += 2;
+  }
+  RecordProperty("avx2_geometry_sum_fallback_checks", static_cast<int>(checks));
+}
+#endif
+
+TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) {
+  const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  const int rounding = std::fegetround();
+  ASSERT_EQ(std::fesetround(FE_UPWARD), 0);
+  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
+            fesa::hex8_internal::EvaluateMaterialPacketScalar);
+  EXPECT_EQ(std::fesetround(rounding), 0);
+#if defined(_M_X64) || defined(__SSE2__)
+  const auto control = _mm_getcsr();
+  _mm_setcsr(control | 0x8040U);
+  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
+            fesa::hex8_internal::EvaluateMaterialPacketScalar);
+  _mm_setcsr(control);
+  _mm_setcsr(control & ~0x80U);
+  const auto unmasked = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  _mm_setcsr(control);
+  EXPECT_EQ(unmasked, fesa::hex8_internal::EvaluateMaterialPacketScalar);
+#endif
+  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(), initial);
+}
+
+TEST(Hex8HighLow, PacketScalarBridgesPreserveOriginalRawBits) {
+  using fesa::hex8_internal::PacketPair;
+  using fesa::hex8_internal::PairPacket;
+  static_assert(std::is_trivial_v<PacketPair> &&
+                std::is_standard_layout_v<PacketPair>);
+  static_assert(std::is_trivial_v<PairPacket> &&
+                std::is_standard_layout_v<PairPacket>);
+  static_assert(std::is_trivial_v<fesa::hex8_internal::MaterialPacketInput>);
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const std::vector<PacketPair> pairs{{0, 0}, {-0.0, 0}, {0, -0.0},
+      {1, -0.0}, {1, std::ldexp(1.0, -54)}, {tiny, 0}, {1, tiny},
+      {std::numeric_limits<double>::max(), 0}, {1, 1},
+      {std::numeric_limits<double>::infinity(), 0},
+      {std::numeric_limits<double>::signaling_NaN(), 0}};
+  for (const auto left : pairs)
+    for (const auto right : pairs) {
+      PairPacket a{{left, left, left, left}};
+      PairPacket b{{right, right, right, right}};
+      PairPacket result;
+      const fesa::high_low_original::HighLow x{left.high, left.low};
+      const fesa::high_low_original::HighLow y{right.high, right.low};
+      const auto expect = [&](fesa::high_low_original::HighLow expected) {
+        for (std::size_t lane = 0; lane < 4; ++lane) {
+          EXPECT_EQ(PairBits(result.values[lane].high), PairBits(expected.high));
+          EXPECT_EQ(PairBits(result.values[lane].low), PairBits(expected.low));
+        }
+      };
+      fesa::hex8_internal::ScalarPacketAdd(a, b, &result);
+      expect(fesa::high_low_original::Add(x, y));
+      fesa::hex8_internal::ScalarPacketMultiply(a, b, &result);
+      expect(fesa::high_low_original::Multiply(x, y));
+      const double high_left[4]{left.high, left.high, left.high, left.high};
+      const double high_right[4]{right.high, right.high, right.high, right.high};
+      fesa::hex8_internal::ScalarPacketTwoProduct(high_left, high_right, &result);
+      expect(fesa::high_low_original::TwoProduct(left.high, right.high));
+    }
+}
+
+TEST(Hex8HighLow, PacketCapabilityRejectsMissingIntegerIsaAndOsState) {
+  constexpr std::uint32_t leaf1 =
+      (1U << 12) | (1U << 26) | (1U << 27) | (1U << 28);
+  constexpr std::uint32_t leaf7 = (1U << 3) | (1U << 5) | (1U << 8);
+  EXPECT_TRUE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
+      leaf1, 6, leaf7));
+  for (const auto missing : {1U << 3, 1U << 5, 1U << 8})
+    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
+        leaf1, 6, leaf7 & ~missing));
+  for (const auto missing : {1U << 12, 1U << 26, 1U << 27, 1U << 28})
+    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
+        leaf1 & ~missing, 6, leaf7));
+  for (std::uint64_t state : {0ULL, 2ULL, 4ULL})
+    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
+        leaf1, state, leaf7));
+}
+
+TEST(Hex8HighLow, PacketDoesNotClearPriorMaskedStickyFlags) {
+  std::fenv_t original;
+  ASSERT_EQ(std::fegetenv(&original), 0);
+  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
+  ASSERT_EQ(std::feraiseexcept(FE_INEXACT), 0);
+  fesa::hex8_internal::MaterialPacketInput input{};
+  input.weight = {1};
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    input.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
+    input.deviatoric_right.values[lane] = {2, std::ldexp(1.0, -54)};
+  }
+  fesa::hex8_internal::PairPacket output;
+  evaluator(input, &output);
+  const bool preserved = (std::fetestexcept(FE_INEXACT) & FE_INEXACT) != 0;
+  EXPECT_EQ(std::fesetenv(&original), 0);
+  EXPECT_TRUE(preserved);
+}
+
 fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
   fesa::ElementDofLayout layout;
   layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7};
@@ -27,6 +847,116 @@ fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
   return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value();
 }
 
+void ExpectSamePairResponse(const fesa::Hex8HighLowResponse& actual,
+                            const fesa::Hex8HighLowResponse& expected) {
+  EXPECT_EQ(PairBits(actual.high.strain_energy),
+            PairBits(expected.high.strain_energy));
+  EXPECT_EQ(PairBits(actual.strain_energy_low),
+            PairBits(expected.strain_energy_low));
+  ASSERT_EQ(actual.high.tangent.Rows(), 24U);
+  ASSERT_EQ(actual.tangent_low.Rows(), 24U);
+  for (std::size_t i = 0; i < 24; ++i) {
+    EXPECT_EQ(PairBits(actual.high.internal_force[i]),
+              PairBits(expected.high.internal_force[i]));
+    EXPECT_EQ(PairBits(actual.internal_force_low[i]),
+              PairBits(expected.internal_force_low[i]));
+    for (std::size_t j = 0; j < 24; ++j) {
+      EXPECT_EQ(PairBits(actual.high.tangent(i, j)),
+                PairBits(expected.high.tangent(i, j)));
+      EXPECT_EQ(PairBits(actual.tangent_low(i, j)),
+                PairBits(expected.tangent_low(i, j)));
+    }
+  }
+  for (std::size_t point = 0; point < 8; ++point) {
+    EXPECT_EQ(actual.high.points[point].integration_point,
+              expected.high.points[point].integration_point);
+    for (std::size_t i = 0; i < 6; ++i) {
+      EXPECT_EQ(PairBits(actual.high.points[point].cauchy_stress[i]),
+                PairBits(expected.high.points[point].cauchy_stress[i]));
+      EXPECT_EQ(PairBits(actual.high.points[point].strain[i]),
+                PairBits(expected.high.points[point].strain[i]));
+    }
+  }
+}
+
+TEST(Hex8HighLow, ConcurrentFirstLinearTrialsAcrossCopiesAndMoveKeepBits) {
+  const auto element = UnitBrick();
+  auto copied = element;
+  const auto moved = std::move(copied);
+  fesa::Vector high(24), low_a(24), low_b(24);
+  for (std::size_t i = 0; i < 24; ++i) high[i] = 1;
+  for (auto a : {1, 2, 5, 6}) {
+    low_a[3U * a] = std::ldexp(1.0, -60);
+    low_b[3U * a] = -std::ldexp(1.0, -59);
+  }
+  const auto oracle = UnitBrick();
+  const auto expected_a = oracle.EvaluateLinearHighLow(high, low_a);
+  const auto expected_b = oracle.EvaluateLinearHighLow(high, low_b);
+  ASSERT_TRUE(expected_a.HasValue());
+  ASSERT_TRUE(expected_b.HasValue());
+  std::promise<void> start;
+  const auto signal = start.get_future().share();
+  std::atomic<std::size_t> arrived{0U};
+  std::array<std::future<fesa::Result<fesa::Hex8HighLowResponse>>, 4> calls;
+  for (std::size_t i = 0; i < calls.size(); ++i)
+    calls[i] = std::async(std::launch::async, [&, i] {
+      arrived.fetch_add(1U, std::memory_order_release);
+      signal.wait();
+      return (i % 2U == 0U ? element : moved)
+          .EvaluateLinearHighLow(high, i < 2U ? low_a : low_b);
+    });
+  while (arrived.load(std::memory_order_acquire) != calls.size())
+    std::this_thread::yield();
+  start.set_value();
+  for (std::size_t i = 0; i < calls.size(); ++i) {
+    const auto response = calls[i].get();
+    ASSERT_TRUE(response.HasValue());
+    ExpectSamePairResponse(response.Value(),
+                           i < 2U ? expected_a.Value() : expected_b.Value());
+  }
+}
+
+TEST(Hex8HighLow, FirstLinearPreparationRecoversAfterTrialAndEnvironmentFailure) {
+  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
+  const auto element = UnitBrick();
+  fesa::Vector high(24), low(24);
+  EXPECT_FALSE(element.EvaluateLinearHighLow(fesa::Vector(23), low).HasValue());
+  high[0] = std::numeric_limits<double>::infinity();
+  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
+  high[0] = 0;
+  struct RestoreEnvironment {
+    int rounding = std::fegetround();
+#if defined(_M_X64) || defined(__SSE2__)
+    unsigned csr = _mm_getcsr();
+#endif
+    ~RestoreEnvironment() {
+      std::fesetround(rounding);
+#if defined(_M_X64) || defined(__SSE2__)
+      _mm_setcsr(csr);
+#endif
+    }
+  } restore;
+  ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0);
+  const auto bad_environment = element.EvaluateLinearHighLow(high, low);
+  ASSERT_FALSE(bad_environment.HasValue());
+  EXPECT_EQ(bad_environment.GetStatus().Diagnostics()[0].message,
+            "Hex8 high+low requires nearest rounding and gradual underflow.");
+  ASSERT_EQ(std::fesetround(FE_TONEAREST), 0);
+#if defined(_M_X64) || defined(__SSE2__)
+  _mm_setcsr(restore.csr | 0x8040U);
+  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
+  _mm_setcsr(restore.csr);
+#endif
+  const auto recovered = element.EvaluateLinearHighLow(high, low);
+  ASSERT_TRUE(recovered.HasValue());
+  high[3] = -2;
+  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
+  high[3] = 0;
+  const auto after_failure = element.EvaluateLinearHighLow(high, low);
+  ASSERT_TRUE(after_failure.HasValue());
+  ExpectSamePairResponse(after_failure.Value(), recovered.Value());
+}
+
 TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) {
   auto brick = UnitBrick();
   fesa::Vector high(24), low(24);

tests/unit/math/dense_blas_internal_test.cpp

@@ -4,6 +4,7 @@
 #include <gtest/gtest.h>
 
 #include <array>
+#include <cmath>
 #include <cstddef>
 #include <limits>
 #include <stdexcept>
@@ -426,3 +427,266 @@ TEST(JetLinearSolve, FailedChannelSolvePreservesCandidateAndInputs) {
 }
 }  // namespace
 }  // namespace fesa::math_internal
+
+namespace fesa::dense_blas_internal {
+TEST(DenseBlasInternal, DiagonalProductKeepsEntrywiseCancellationAndChecks) {
+  Vector diagonal{3U};
+  diagonal[0U] = 4;
+  diagonal[1U] = 2;
+  diagonal[2U] = -1;
+  Vector input{3U};
+  const double large = std::ldexp(1.0, 56);
+  input[0U] = 1;
+  input[1U] = -large;
+  input[2U] = 2 * large;
+  const auto product = MultiplyDiagonal(diagonal, input);
+  EXPECT_DOUBLE_EQ(product[0U] + (product[1U] - product[2U]), 4);
+  EXPECT_EQ(MultiplyDiagonal(Vector{0U}, Vector{0U}).Size(), 0U);
+  EXPECT_THROW(static_cast<void>(MultiplyDiagonal(diagonal, Vector{2U})),
+               std::invalid_argument);
+  input[0U] = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(MultiplyDiagonal(diagonal, input)),
+               std::invalid_argument);
+  input[0U] = 1;
+  diagonal[0U] = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(static_cast<void>(MultiplyDiagonal(diagonal, input)),
+               std::invalid_argument);
+}
+
+TEST(DenseBlasInternal, LocalDenseScopeRestoresTlsAfterProductsAndUnwind) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int default_threads = mkl_get_max_threads();
+  const std::array<double, 2U> left{0.0, 2.0};
+  const std::array<double, 2U> right{3.0, 4.0};
+  std::array<double, 4U> product{};
+  product.fill(std::numeric_limits<double>::quiet_NaN());
+  {
+    const ScopedLocalDenseThreads scope{};
+    EXPECT_EQ(mkl_get_max_threads(), 1);
+    EXPECT_TRUE(
+        OuterProduct(left.data(), 2U, right.data(), 2U, product.data()));
+    EXPECT_DOUBLE_EQ(product[0U], 0.0);
+    EXPECT_DOUBLE_EQ(product[3U], 8.0);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  try {
+    const ScopedLocalDenseThreads scope{};
+    OuterProduct(nullptr, 2U, right.data(), 2U, product.data());
+  } catch (const std::invalid_argument&) {
+  }
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  static_cast<void>(mkl_set_num_threads_local(3));
+  {
+    const ScopedLocalDenseThreads scope{};
+    EXPECT_EQ(mkl_get_max_threads(), 3);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), 3);
+}
+
+TEST(DenseBlasInternal, OuterProductLeavesFiniteInputOverflowForOwnerGate) {
+  const std::array<double, 1U> left{std::numeric_limits<double>::max()};
+  const std::array<double, 1U> right{2.0};
+  std::array<double, 1U> product{};
+  EXPECT_TRUE(OuterProduct(left.data(), 1U, right.data(), 1U, product.data()));
+  EXPECT_TRUE(std::isinf(product[0U]));
+}
+
+TEST(DenseBlasInternal,
+     OuterProductUsesRectangularStorageAndPreservesRejection) {
+  const std::array<double, 2U> left{2.0, -3.0};
+  const std::array<double, 3U> right{4.0, -5.0, 0.25};
+  std::array<double, 6U> product{};
+  ASSERT_TRUE(OuterProduct(left.data(), left.size(), right.data(), right.size(),
+                           product.data()));
+  const std::array<double, 6U> expected{8.0, -10.0, 0.5, -12.0, 15.0, -0.75};
+  EXPECT_EQ(product, expected);
+  const auto accepted = product;
+  const std::array<double, 2U> invalid{
+      1.0, std::numeric_limits<double>::quiet_NaN()};
+  EXPECT_FALSE(OuterProduct(invalid.data(), invalid.size(), right.data(),
+                            right.size(), product.data()));
+  EXPECT_EQ(product, accepted);
+  EXPECT_TRUE(OuterProduct(nullptr, 0U, nullptr, 3U, nullptr));
+  EXPECT_TRUE(OuterProduct(nullptr, 2U, nullptr, 0U, nullptr));
+  EXPECT_THROW(OuterProduct(nullptr, 2U, right.data(), 3U, product.data()),
+               std::invalid_argument);
+  const auto too_large =
+      static_cast<std::size_t>((std::numeric_limits<MKL_INT>::max)()) + 1U;
+  EXPECT_THROW(OuterProduct(nullptr, too_large, nullptr, 1U, nullptr),
+               std::length_error);
+  const auto too_many_elements =
+      std::numeric_limits<std::size_t>::max() / sizeof(double) + 1U;
+  EXPECT_THROW(OuterProduct(nullptr, too_many_elements, nullptr, 1U, nullptr),
+               std::length_error);
+}
+
+TEST(DenseBlasInternal, PolicyQueryPreservesDefaultPositiveAndNestedTls) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int global = mkl_get_max_threads();
+  const int domain = mkl_domain_get_max_threads(MKL_DOMAIN_BLAS);
+  const int dynamic = mkl_get_dynamic();
+  for (int repeat = 0; repeat < 2; ++repeat) {
+    const auto policy = ReadScopedLocalDensePolicy();
+    EXPECT_EQ(policy.blas_max_threads, 1);
+    EXPECT_EQ(policy.dynamic, dynamic);
+    EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+    EXPECT_EQ(mkl_get_max_threads(), global);
+    EXPECT_EQ(mkl_domain_get_max_threads(MKL_DOMAIN_BLAS), domain);
+    EXPECT_EQ(mkl_get_dynamic(), dynamic);
+  }
+  static_cast<void>(mkl_set_num_threads_local(2));
+  EXPECT_EQ(ReadScopedLocalDensePolicy().blas_max_threads, 2);
+  EXPECT_EQ(mkl_set_num_threads_local(2), 2);
+  static_cast<void>(mkl_set_num_threads_local(0));
+  try {
+    const ScopedLocalDenseThreads outer{};
+    EXPECT_EQ(ReadScopedLocalDensePolicy().blas_max_threads, 1);
+    EXPECT_EQ(mkl_set_num_threads_local(1), 1);
+    throw std::runtime_error{"policy-query-unwind"};
+  } catch (const std::runtime_error&) {
+  }
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  EXPECT_EQ(mkl_get_dynamic(), dynamic);
+}
+
+TEST(DenseBlasInternal, SingleEquationScopeRestoresTlsAndRespectsOverrides) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int default_threads = mkl_get_max_threads();
+  {
+    const ScopedSingleEquationThreads scope{1U};
+    EXPECT_EQ(mkl_get_max_threads(), 1);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  {
+    const ScopedSingleEquationThreads scope{2U};
+    EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  }
+  try {
+    const ScopedSingleEquationThreads scope{1U};
+    throw std::runtime_error{"scope failure"};
+  } catch (const std::runtime_error&) {
+  }
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  static_cast<void>(mkl_set_num_threads_local(3));
+  {
+    const ScopedSingleEquationThreads scope{1U};
+    EXPECT_EQ(mkl_get_max_threads(), 3);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), 3);
+}
+
+TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) {
+  Matrix left{3U, 2U};
+  Matrix right{3U, 1U};
+  for (std::size_t row = 0; row < 3U; ++row) {
+    left(row, 0U) = static_cast<double>(row + 1U);
+    left(row, 1U) = 2.0;
+    right(row, 0U) = static_cast<double>(row);
+  }
+  const auto product = MultiplyTransposeLeft(left, right);
+  ASSERT_EQ(product.Rows(), 2U);
+  ASSERT_EQ(product.Columns(), 1U);
+  EXPECT_DOUBLE_EQ(product(0U, 0U), 8.0);
+  EXPECT_DOUBLE_EQ(product(1U, 0U), 6.0);
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, Matrix{2U, 1U})),
+               std::invalid_argument);
+  const auto empty = MultiplyTransposeLeft(Matrix{0U, 2U}, Matrix{0U, 3U});
+  EXPECT_EQ(empty.Rows(), 2U);
+  EXPECT_EQ(empty.Columns(), 3U);
+  EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0);
+  left(0U, 0U) = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, right)),
+               std::invalid_argument);
+  left(0U, 0U) = 1.0;
+  right(0U, 0U) = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, right)),
+               std::invalid_argument);
+}
+
+TEST(DenseBlasInternal, TransposeRightProductUsesInnerAndRectangularStrides) {
+  const std::array<double, 6U> left{1, 2, 3, 4, 5, 6};
+  const std::array<double, 12U> right{7,  8,  9,  10, 11, 12,
+                                      13, 14, 15, 16, 17, 18};
+  std::array<double, 8U> product;
+  ASSERT_TRUE(MultiplyTransposeRight(left.data(), 2U, right.data(), 4U, 3U,
+                                     product.data()));
+  const std::array<double, 8U> expected{50, 68, 86, 104, 122, 167, 212, 257};
+  EXPECT_EQ(product, expected);
+  EXPECT_TRUE(
+      MultiplyTransposeRight(nullptr, 2U, nullptr, 4U, 0U, product.data()));
+  for (const auto value : product) EXPECT_DOUBLE_EQ(value, 0.0);
+  auto nonfinite = right;
+  nonfinite[8U] = std::numeric_limits<double>::infinity();
+  product.fill(31.0);
+  EXPECT_FALSE(MultiplyTransposeRight(left.data(), 2U, nonfinite.data(), 4U, 3U,
+                                      product.data()));
+  for (const auto value : product) EXPECT_DOUBLE_EQ(value, 31.0);
+}
+
+// The former lower/inverse output parameters are now owned factor results.
+TEST(DenseBlasInternal, CholeskyAndInverseMatchKnownFactorsAndKeepFailedOutput) {
+  Matrix matrix{2U, 2U};
+  matrix(0, 0) = 4;
+  matrix(0, 1) = matrix(1, 0) = 2;
+  matrix(1, 1) = 5;
+  auto factor = DenseFactor::Cholesky(matrix);
+  ASSERT_TRUE(factor.HasValue());
+  const auto lower = factor.Value().Lower();
+  Matrix identity{2U, 2U};
+  identity(0, 0) = identity(1, 1) = 1;
+  const auto inverse = factor.Value().SolveLower(identity);
+  ASSERT_TRUE(lower.HasValue());
+  ASSERT_TRUE(inverse.HasValue());
+  EXPECT_DOUBLE_EQ(lower.Value()(0, 0), 2);
+  EXPECT_DOUBLE_EQ(lower.Value()(1, 0), 1);
+  EXPECT_DOUBLE_EQ(lower.Value()(1, 1), 2);
+  EXPECT_DOUBLE_EQ(lower.Value()(0, 1), 0);
+  EXPECT_DOUBLE_EQ(inverse.Value()(0, 0), .5);
+  EXPECT_DOUBLE_EQ(inverse.Value()(1, 0), -.25);
+  EXPECT_DOUBLE_EQ(inverse.Value()(1, 1), .5);
+  matrix(1, 1) = 0;
+  EXPECT_FALSE(DenseFactor::Cholesky(matrix).HasValue());
+  EXPECT_DOUBLE_EQ(lower.Value()(1, 1), 2);
+  EXPECT_DOUBLE_EQ(inverse.Value()(1, 1), .5);
+  EXPECT_FALSE(DenseFactor::Cholesky(Matrix{2U, 1U}).HasValue());
+  matrix(0, 0) = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_FALSE(DenseFactor::Cholesky(matrix).HasValue());
+  EXPECT_DOUBLE_EQ(factor.Value().SolveLower(identity).Value()(1, 1), .5);
+}
+
+TEST(DenseBlasInternal,
+     CholeskyFailureReasonPreservesCandidatesAndResetsOnReuse) {
+  Matrix matrix{2U, 2U};
+  matrix(0, 0) = 4;
+  matrix(1, 1) = 9;
+  const auto accepted = DenseFactor::Cholesky(matrix);
+  ASSERT_TRUE(accepted.HasValue());
+  Matrix identity{2U, 2U};
+  identity(0, 0) = identity(1, 1) = 1;
+  matrix(1, 1) = -1;
+  const auto failed = DenseFactor::Cholesky(matrix);
+  ASSERT_FALSE(failed.HasValue());
+  ASSERT_FALSE(failed.GetStatus().Diagnostics().empty());
+  EXPECT_EQ(failed.GetStatus().Diagnostics()[0].code, "dense-factor-pivot");
+  EXPECT_DOUBLE_EQ(accepted.Value().Lower().Value()(1, 1), 3);
+  EXPECT_DOUBLE_EQ(accepted.Value().SolveLower(identity).Value()(1, 1), 1.0 / 3.0);
+  matrix(1, 1) = 9;
+  EXPECT_TRUE(DenseFactor::Cholesky(matrix).HasValue());
+  EXPECT_FALSE(DenseFactor::Cholesky(Matrix{2U, 1U}).HasValue());
+  EXPECT_TRUE(DenseFactor::Cholesky(matrix).HasValue());
+}
+}  // namespace fesa::dense_blas_internal

tests/unit/math/high_low_arithmetic_test.cpp

@@ -1,10 +1,16 @@
 #include "../../../src/fesa/math/high_low_arithmetic.h"
-#include "elements/hex8_pair_simd.h"
 
 #include <gtest/gtest.h>
+#include "elements/hex8_pair_simd.h"
 
 #include <cmath>
+#include <cstdint>
+#include <cstring>
 #include <limits>
+#include <random>
+#include <vector>
+
+#include "high_low_arithmetic_original.h"
 
 namespace {
 using namespace fesa::math_internal;
@@ -126,9 +132,111 @@ TEST(HighLowArithmetic, RejectsChangedRoundingAndFlushToZeroEnvironment) {
   EXPECT_FALSE(ArithmeticEnvironmentSupported());
 #endif
 }
-}  // namespace
 
-namespace {
+std::uint64_t Bits(const double value) {
+  std::uint64_t result = 0;
+  std::memcpy(&result, &value, sizeof(result));
+  return result;
+}
+
+double FromBits(const std::uint64_t bits) {
+  double result = 0;
+  std::memcpy(&result, &bits, sizeof(result));
+  return result;
+}
+
+TEST(HighLowArithmetic, FiniteClassificationMatchesBinary64Boundaries) {
+  const auto check = [](const double value) {
+    const bool expected = std::isfinite(value);
+    EXPECT_EQ(IsFinite(HighLow{value, 0.0}), expected);
+    EXPECT_EQ(IsFinite(HighLow{1.0, value}), expected);
+  };
+  for (std::uint64_t exponent = 0; exponent < 2048U; ++exponent) {
+    for (const auto fraction : {0ULL, 1ULL, 1ULL << 51U,
+                                (1ULL << 52U) - 1ULL}) {
+      for (const auto sign : {0ULL, 1ULL << 63U})
+        check(FromBits(sign | (exponent << 52U) | fraction));
+    }
+  }
+  std::mt19937_64 random{0x5800f17U};
+  for (int sample = 0; sample < 10000; ++sample) check(FromBits(random()));
+}
+
+void ExpectSamePair(const HighLow actual,
+                    const fesa::high_low_original::HighLow expected) {
+  EXPECT_EQ(Bits(actual.high), Bits(expected.high));
+  EXPECT_EQ(Bits(actual.low), Bits(expected.low));
+}
+
+TEST(HighLowArithmetic, BitExponentMatchesFrozenOracleAcrossEveryExponent) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  // Every normal exponent, every possible trailing significand-bit count,
+  // both signs; then every subnormal single-bit and mixed-bit boundary.
+  for (std::uint64_t exponent = 1; exponent < 2047; ++exponent) {
+    for (unsigned trailing = 0; trailing <= 52U; ++trailing) {
+      const auto fraction = trailing == 52U ? 0ULL : 1ULL << trailing;
+      for (const auto sign : {0ULL, 1ULL << 63U}) {
+        const double value = FromBits(sign | (exponent << 52U) | fraction);
+        ASSERT_EQ(LowestBitExponent(value),
+                  fesa::high_low_original::LowestBitExponent(value));
+      }
+    }
+  }
+  for (unsigned bit = 0; bit < 52U; ++bit) {
+    for (const auto significand : {1ULL << bit, (1ULL << 52U) - (1ULL << bit)}) {
+      const double value = FromBits(significand);
+      EXPECT_EQ(LowestBitExponent(value),
+                fesa::high_low_original::LowestBitExponent(value));
+    }
+  }
+}
+
+TEST(HighLowArithmetic, OptimizedArithmeticMatchesFrozenOriginalBits) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  std::vector<double> values{0.0, -0.0,
+      std::numeric_limits<double>::denorm_min(),
+      -std::numeric_limits<double>::denorm_min(),
+      std::numeric_limits<double>::min(),
+      std::nextafter(std::numeric_limits<double>::min(), 0.0),
+      std::numeric_limits<double>::max(),
+      -std::numeric_limits<double>::max(), 1.0, -1.0,
+      std::nextafter(1.0, 2.0), std::ldexp(1.0, -53),
+      std::ldexp(1.0, -54), std::numeric_limits<double>::infinity(),
+      -std::numeric_limits<double>::infinity(),
+      std::numeric_limits<double>::quiet_NaN()};
+  std::mt19937_64 random{0x58a0fU};
+  for (int sample = 0; sample < 4000; ++sample) {
+    values.push_back(FromBits(random()));
+  }
+  for (std::size_t index = 0; index < values.size(); ++index) {
+    const double high = values[index];
+    const double low = values[(index * 13U + 7U) % values.size()];
+    const fesa::high_low_original::HighLow original{high, low};
+    EXPECT_EQ(IsNormalized({high, low}),
+              fesa::high_low_original::IsNormalized(original));
+    // Normalize a raw pair with the frozen operation, never the new predicate.
+    const auto a = fesa::high_low_original::TwoSum(high, low);
+    const auto b = fesa::high_low_original::TwoSum(
+        values[(index + 1U) % values.size()],
+        values[(index * 3U + 2U) % values.size()]);
+    ExpectSamePair(Add({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Add(a, b));
+    ExpectSamePair(Multiply({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Multiply(a, b));
+    ExpectSamePair(Divide({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Divide(a, b));
+    ExpectSamePair(TwoProduct(high, low),
+                   fesa::high_low_original::TwoProduct(high, low));
+  }
+  for (double high : {1.0, std::nextafter(1.0, 2.0), -1.0,
+                      std::numeric_limits<double>::max()}) {
+    for (double low : {0.0, -0.0, std::ldexp(1.0, -53),
+                       -std::ldexp(1.0, -53)}) {
+      EXPECT_EQ(IsNormalized({high, low}),
+                fesa::high_low_original::IsNormalized({high, low}));
+    }
+  }
+}
 TEST(HighLowArithmetic, PacketDispatchRetainsComponentsAndInstructionGates) {
   using namespace fesa::hex8_internal;
   const std::uint32_t leaf1 = (1U << 12) | (1U << 26) | (1U << 27) | (1U << 28);

tests/unit/math/matrix_test.cpp

@@ -167,4 +167,21 @@ TEST(DenseMath, MatrixSelfAssignmentPreservesShapeAndStorage) {
 }
 
 }  // namespace
+TEST(DenseMath, MatrixProductsRejectNonfiniteOperandsAndPreserveValidReuse) {
+  Matrix left{2U, 2U, 1.0};
+  Matrix right{2U, 1U, 2.0};
+  Vector rhs{2U, 2.0};
+  rhs[1U] = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(left.Multiply(rhs)), std::invalid_argument);
+  rhs[1U] = 2.0;
+  left(0U, 1U) = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(static_cast<void>(left.Multiply(rhs)), std::invalid_argument);
+  EXPECT_THROW(static_cast<void>(left.Multiply(right)), std::invalid_argument);
+  left(0U, 1U) = 1.0;
+  right(1U, 0U) = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(left.Multiply(right)), std::invalid_argument);
+  right(1U, 0U) = 2.0;
+  EXPECT_DOUBLE_EQ(left.Multiply(rhs)[0U], 4.0);
+  EXPECT_DOUBLE_EQ(left.Multiply(right)(1U, 0U), 4.0);
+}
 }  // namespace fesa

tests/unit/math/sparse_matrix_test.cpp

@@ -3,6 +3,8 @@
 #include <gtest/gtest.h>
 
 #include <algorithm>
+#include <cmath>
+#include <future>
 #include <limits>
 #include <type_traits>
 #include <utility>
@@ -131,4 +133,104 @@ TEST(SparseAssembly, PreservesExpectedStructuralZeros) {
             4);
 }
 
+TEST(SparseAssembly, MultipliesRectangularCopyMoveAndConcurrentConstCalls) {
+  auto assembled = SparseMatrix::FromCoo(
+      3U, 5U, {{0U, 4U, 3.0, 0U, 0U}, {2U, 1U, -2.0, 1U, 0U}},
+      {{0U, 1U, 1U, 2U}, {4U, 1U}});
+  ASSERT_TRUE(assembled.HasValue());
+  const auto original = assembled.Value();
+  auto copied = original;
+  auto moved = std::move(copied);
+  Vector rhs{5U, 2.0};
+  const auto evaluate = [&rhs](const SparseMatrix& matrix) {
+    for (int repeat = 0; repeat < 20; ++repeat) {
+      const auto product = matrix.Multiply(rhs);
+      if (product.Size() != 3U || product[0U] != 6.0 ||
+          product[1U] != 0.0 || product[2U] != -4.0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, [&] { return evaluate(moved); });
+  auto second = std::async(std::launch::async, [&] { return evaluate(moved); });
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+  EXPECT_TRUE(evaluate(original));
+  auto replacement = SparseMatrix::FromCoo(
+      1U, 1U, {{0U, 0U, 7.0, 0U, 0U}}, {{0U, 1U}, {0U}});
+  ASSERT_TRUE(replacement.HasValue());
+  moved = replacement.Value();
+  EXPECT_DOUBLE_EQ(moved.Multiply(Vector{1U, 3.0})[0U], 21.0);
+  EXPECT_TRUE(evaluate(original));
+}
+
+TEST(SparseAssembly, PreservesZeroShapesAndEmptyRows) {
+  auto no_rows = SparseMatrix::FromCoo(0U, 3U, {}, {{0U}, {}});
+  auto no_columns = SparseMatrix::FromCoo(2U, 0U, {}, {{0U, 0U, 0U}, {}});
+  auto no_entries = SparseMatrix::FromCoo(2U, 3U, {}, {{0U, 0U, 0U}, {}});
+  ASSERT_TRUE(no_rows.HasValue());
+  ASSERT_TRUE(no_columns.HasValue());
+  ASSERT_TRUE(no_entries.HasValue());
+  EXPECT_EQ(no_rows.Value().Multiply(Vector{3U}).Size(), 0U);
+  const auto empty_columns = no_columns.Value().Multiply(Vector{0U});
+  const auto empty_entries = no_entries.Value().Multiply(Vector{3U});
+  ASSERT_EQ(empty_columns.Size(), 2U);
+  ASSERT_EQ(empty_entries.Size(), 2U);
+  for (std::size_t row = 0; row < 2U; ++row) {
+    EXPECT_DOUBLE_EQ(empty_columns[row], 0.0);
+    EXPECT_DOUBLE_EQ(empty_entries[row], 0.0);
+  }
+}
+
+TEST(SparseAssembly, RejectsNonfiniteRhsAndRecoversAfterFiniteOverflow) {
+  auto matrix = SparseMatrix::FromCoo(
+      2U, 2U, {{0U, 0U, 2.0, 0U, 0U}}, {{0U, 1U, 2U}, {0U, 1U}});
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{2U, (std::numeric_limits<double>::infinity)()};
+  EXPECT_THROW(static_cast<void>(matrix.Value().Multiply(rhs)),
+               std::invalid_argument);
+  rhs[1U] = 4.0;
+  rhs[0U] = (std::numeric_limits<double>::max)();
+  EXPECT_TRUE(std::isinf(matrix.Value().Multiply(rhs)[0U]));
+  rhs[0U] = (std::numeric_limits<double>::quiet_NaN)();
+  EXPECT_THROW(static_cast<void>(matrix.Value().Multiply(rhs)),
+               std::invalid_argument);
+  EXPECT_THROW(static_cast<void>(matrix.Value().Multiply(Vector{1U})),
+               std::invalid_argument);
+  rhs[0U] = 3.0;
+  rhs[1U] = 4.0;
+  const auto finite = matrix.Value().Multiply(rhs);
+  EXPECT_DOUBLE_EQ(finite[0U], 6.0);
+  EXPECT_DOUBLE_EQ(finite[1U], 0.0);
+}
+
+TEST(SparseAssembly, MatchesIndependentLongDoubleResidualOnUnevenRows) {
+  const std::vector<double> entries{1.25, -0.375, 4.5, -3.0, 0.125, 2.0};
+  const SparsePattern pattern{{0U, 3U, 3U, 5U, 6U}, {0U, 2U, 4U, 1U, 3U, 2U}};
+  std::vector<CooContribution> contributions;
+  for (std::size_t row = 0; row < 4U; ++row) {
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      contributions.push_back(
+          {row, pattern.column_indices[position], entries[position], row,
+           position});
+    }
+  }
+  auto matrix = SparseMatrix::FromCoo(4U, 5U, contributions, pattern);
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{5U};
+  for (std::size_t column = 0; column < 5U; ++column) {
+    rhs[column] = 0.3 + static_cast<double>(column) * 0.7;
+  }
+  const auto product = matrix.Value().Multiply(rhs);
+  for (std::size_t row = 0; row < 4U; ++row) {
+    long double expected = 0;
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      expected += static_cast<long double>(entries[position]) *
+                  rhs[pattern.column_indices[position]];
+    }
+    EXPECT_NEAR(product[row], static_cast<double>(expected), 2.0e-14);
+  }
+}
+
 }  // namespace

tests/unit/results/result_recovery_components_test.cpp

@@ -6,6 +6,7 @@
 #include <filesystem>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -353,6 +354,41 @@ TEST(ResultRecoveryComponents,
   }
 }
 
+TEST(ResultRecoveryComponents, PartitionedScaleKeepsCancelledInternalTerms) {
+  auto definition = MakeShellDefinition();
+  definition.node_sets = {{"Fixed", std::optional<std::string>{"Shell-1"},
+                           {0U}, {definition.source_path, 60U}}};
+  definition.steps[0].boundaries = {
+      {"Fixed", 1, 1, 1.0, {definition.source_path, 61U}}};
+  const auto fixture = MakeShellFixture(std::move(definition));
+  const auto count = fixture.dofs->FullDofCount();
+  fesa::SparsePattern pattern;
+  pattern.row_offsets = {0U, 1U, 3U};
+  pattern.row_offsets.resize(count + 1U, 3U);
+  pattern.column_indices = {0U, 0U, 1U};
+  auto matrix = fesa::SparseMatrix::FromCoo(
+      count, count,
+      {{0U, 0U, 1.0, 0U, 0U}, {1U, 0U, -1.0, 0U, 1U},
+       {1U, 1U, 1.0, 0U, 2U}}, pattern);
+  ASSERT_TRUE(matrix.HasValue());
+  fesa::Vector displacement{count};
+  displacement[0U] = 1.0;
+  displacement[1U] = 1.0;
+  fesa::Vector external{count};
+  external[1U] = 1.0e-12;
+  const auto result = fesa::results_internal::GlobalEquilibriumRecovery(
+      *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external);
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_DOUBLE_EQ(result.Value().internal_force[1U], 0.0);
+  EXPECT_DOUBLE_EQ(result.Value().normalized_free_residual, 1.0e-12);
+  displacement[1U] = std::numeric_limits<double>::quiet_NaN();
+  const auto rejected = fesa::results_internal::GlobalEquilibriumRecovery(
+      *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external);
+  ASSERT_FALSE(rejected.HasValue());
+  ExpectStatusCode(rejected.GetStatus(), "nonfinite-recovery-value");
+  EXPECT_DOUBLE_EQ(result.Value().reaction[1U], -1.0e-12);
+}
+
 // C-MODULE-002
 TEST(ResultRecoveryComponents,
      BeamResultRecoveryKeepsEndpointGaussAndStressIdentities) {

tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp

@@ -3,6 +3,7 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
+#include <future>
 #include <initializer_list>
 #include <limits>
 #include <string>
@@ -11,6 +12,7 @@
 
 #include "fesa/fem/dof_manager.h"
 #include "fesa/math/sparse_matrix.h"
+#include "math/dense_blas_internal.h"
 
 namespace {
 
@@ -123,6 +125,42 @@ TEST(MklPardisoSolver, ReusesOneFactorizationForRepeatedRhs) {
   EXPECT_LE(NormalizedResidual(matrix, second, rhs_second), 1.0e-10);
 }
 
+TEST(MklPardisoSolver, SingleEquationPolicyPreservesLifecycleAndConcurrentTls) {
+  const auto worker = [] {
+    const int previous = mkl_set_num_threads_local(0);
+    struct Restore {
+      int previous;
+      ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+    } restore{previous};
+    const int default_threads = mkl_get_max_threads();
+    fesa::MklPardisoSolver solver;
+    for (int repeat = 0; repeat < 10; ++repeat) {
+      if (!solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk()) return false;
+      auto solution = MakeVector({9.0});
+      if (!solver.Solve(MakeVector({4.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (!solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk()) return false;
+      if (!solver.Solve(MakeVector({4.0}), solution).IsOk() ||
+          solution[0U] != 1.0) return false;
+      auto invalid_rhs = MakeVector({std::numeric_limits<double>::quiet_NaN()});
+      if (solver.Solve(invalid_rhs, solution).IsOk() || solution[0U] != 1.0)
+        return false;
+      if (!solver.Solve(MakeVector({8.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (solver.Factorize(MakeDenseCsr(1U, {-1.0})).IsOk()) return false;
+      if (solver.Solve(MakeVector({8.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (mkl_get_max_threads() != default_threads ||
+          mkl_set_num_threads_local(0) != 0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, worker);
+  auto second = std::async(std::launch::async, worker);
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+}
+
 TEST(MklPardisoSolver, StillRejectsNonsymmetricPublicCsr) {
   const auto nonsymmetric = MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0});
   fesa::MklPardisoSolver solver;
@@ -155,6 +193,32 @@ TEST(MklPardisoSolver, RefactorizesWithoutLeakingState) {
             1.0e-10);
 }
 
+TEST(MklPardisoSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int default_threads = mkl_get_max_threads();
+  fesa::MklPardisoSolver solver;
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk());
+  auto single = MakeVector({0.0});
+  ASSERT_TRUE(solver.Solve(MakeVector({4.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(2U, {2, 0, 0, 3})).IsOk());
+  fesa::Vector pair{2U};
+  ASSERT_TRUE(solver.Solve(MakeVector({4, 6}), pair).IsOk());
+  EXPECT_NEAR(pair[0U], 2, 1e-14);
+  EXPECT_NEAR(pair[1U], 2, 1e-14);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk());
+  ASSERT_TRUE(solver.Solve(MakeVector({8.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+}
+
 TEST(MklPardisoSolver, ClassifiesSingularIndefiniteAndNonfiniteFailures) {
   fesa::MklPardisoSolver solver;
 

tests/unit/solvers/linear/tangent_solver_test.cpp

@@ -3,6 +3,7 @@
 #include <gtest/gtest.h>
 
 #include <cstddef>
+#include <future>
 #include <initializer_list>
 #include <limits>
 #include <type_traits>
@@ -13,6 +14,7 @@
 #include "fesa/math/sparse_matrix.h"
 #include "fesa/math/vector.h"
 #include "fesa/solvers/linear/tangent_solver.h"
+#include "math/dense_blas_internal.h"
 
 namespace {
 
@@ -129,4 +131,104 @@ TEST(TangentSolver, RejectsSolveBeforeSuccessfulFactorization) {
   EXPECT_DOUBLE_EQ(correction[1U], -9.0);
 }
 
+TEST(TangentSolver, ReusesRhsButFreshValuesAndFailureCannotReuseStaleState) {
+  fesa::MklPardisoTangentSolver solver;
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0}))
+                  .IsOk());
+  fesa::Vector correction{2U};
+  for (double scale : {1.0, -2.0, 0.5}) {
+    ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0 * scale, 6.0 * scale}),
+                                      correction).IsOk());
+    EXPECT_NEAR(correction[0U], scale, 1.0e-14);
+    EXPECT_NEAR(correction[1U], 2.0 * scale, 1.0e-14);
+  }
+  const auto bad_rhs = MakeVector({std::numeric_limits<double>::quiet_NaN(), 0});
+  correction = MakeVector({17, -4});
+  EXPECT_FALSE(solver.SolveCorrection(bad_rhs, correction).IsOk());
+  EXPECT_DOUBLE_EQ(correction[0U], 17);
+  EXPECT_DOUBLE_EQ(correction[1U], -4);
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 6}), correction).IsOk());
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {4, 0, 1, 2})).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 5}), correction).IsOk());
+  EXPECT_NEAR(correction[0U], 1, 1.0e-14);
+  EXPECT_NEAR(correction[1U], 2, 1.0e-14);
+  auto diagonal = fesa::SparseMatrix::FromCoo(
+      2U, 2U, {{0U, 0U, 5, 0U, 0U}, {1U, 1U, 7, 1U, 0U}},
+      {{0U, 1U, 2U}, {0U, 1U}});
+  ASSERT_TRUE(diagonal.HasValue());
+  ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk());
+  EXPECT_NEAR(correction[0U], 2, 1.0e-14);
+  EXPECT_NEAR(correction[1U], -1, 1.0e-14);
+  EXPECT_FALSE(solver.FactorizeFresh(MakeDenseCsr(2U, {0, 0, 0, 0})).IsOk());
+  correction = MakeVector({19, -8});
+  const auto stale = solver.SolveCorrection(MakeVector({10, -7}), correction);
+  ASSERT_FALSE(stale.IsOk());
+  EXPECT_EQ(stale.Diagnostics().front().code, "solver-not-factorized");
+  EXPECT_DOUBLE_EQ(correction[0U], 19);
+  EXPECT_DOUBLE_EQ(correction[1U], -8);
+  ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk());
+  EXPECT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk());
+}
+
+TEST(TangentSolver, SingleEquationScopeRestoresTlsAcrossFreshFailureAndReuse) {
+  const auto worker = [] {
+    const int previous = mkl_set_num_threads_local(0);
+    struct Restore {
+      int previous;
+      ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+    } restore{previous};
+    const int default_threads = mkl_get_max_threads();
+    fesa::MklPardisoTangentSolver solver;
+    for (int repeat = 0; repeat < 10; ++repeat) {
+      auto correction = MakeVector({17.0});
+      if (!solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk()) return false;
+      if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 2.0) return false;
+      if (!solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk()) return false;
+      if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (solver.SolveCorrection(MakeVector({
+              std::numeric_limits<double>::infinity()}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (solver.FactorizeFresh(MakeDenseCsr(1U, {0.0})).IsOk()) return false;
+      if (solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (mkl_get_max_threads() != default_threads ||
+          mkl_set_num_threads_local(0) != 0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, worker);
+  auto second = std::async(std::launch::async, worker);
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+}
+
+TEST(TangentSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int default_threads = mkl_get_max_threads();
+  fesa::MklPardisoTangentSolver solver;
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk());
+  auto single = MakeVector({0.0});
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2, 1, 0, 3})).IsOk());
+  fesa::Vector pair{2U};
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({6, 6}), pair).IsOk());
+  EXPECT_NEAR(pair[0U], 2, 1e-14);
+  EXPECT_NEAR(pair[1U], 2, 1e-14);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({8.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+}
+
 }  // namespace