Merge branch 'main' into codex/verifier-integrity-538-500

Author: morluto

.github/scripts/manage-test-timings

@@ -300,10 +300,16 @@ def merge(
         payload = load_json_bytes(path.read_bytes(), str(path))
         durations = validate_durations(payload, str(path), suite)
         duplicates = merged.keys() & durations.keys()
-        if duplicates:
-            duplicate = min(duplicates)
-            raise ValueError(f"duplicate timing entry across shards: {duplicate}")
-        merged.update(durations)
+        for nodeid in sorted(duplicates):
+            previous = merged[nodeid]
+            incoming = durations[nodeid]
+            if incoming > previous:
+                merged[nodeid] = incoming
+            warning(
+                f"duplicate timing entry across shards for {nodeid}; "
+                f"keeping max({previous}, {incoming})"
+            )
+        merged.update({k: v for k, v in durations.items() if k not in duplicates})
 
     envelope = {
         "version": 1,

docs/reference/capabilities/polynomial/gaussian-polynomial-moments.md

@@ -30,7 +30,7 @@ The request bounds are:
 | Nonzero sparse terms | 1–16 |
 | Total degree of each input term | at most 8 |
 | Fixed moment order \(m\) | 0–16 |
-| Raw ordered expansion paths | `term_count ** m <= 4096` |
+| Raw ordered expansion paths | `term_count ** m <= 65536` |
 | Input rational numerator/denominator | at most 128 decimal digits |
 
 Terms use exponent vectors of length `variable_count`, in strictly increasing
@@ -39,8 +39,10 @@ lexicographic order. A coefficient has separate canonical rational `real` and
 exponents, inconsistent dimensions, and requests beyond the complete-expansion
 bound fail validation before computation or artifact writes.
 
-The producer exactly expands \(P^m\), merges equal exponent vectors, and removes
-zero coefficients. For every remaining exponent vector
+The producer exactly expands \(P^m\) via pinned Python-FLINT `fmpq_mpoly`
+binary exponentiation over a typed real/imaginary coefficient pair,
+merges equal exponent vectors, and removes zero coefficients. For every
+remaining exponent vector
 \(\alpha=(\alpha_1,\ldots,\alpha_n)\), its ledger records:
 
 - the exact expanded coefficient;

docs/reference/scenarios/certified-smith-integral-homology.md

@@ -3,9 +3,11 @@
 [Documentation home](../../index.md)
 
 - Status: Current implementation reference; contracts are experimental
-- Matrix producer backend: bounded standard-library elementary operations
+- Matrix producer backend: SymPy `smith_normal_decomp` over `ZZ`; canonical
+  diagonal and invariant factors, noncanonical unimodular transformations
 - Topology producer backend: exact integer chain arithmetic over canonical
-  simplex bases
+  simplex bases; integral homology derives cycle bases and generator
+  coordinates from the certified Smith transformations
 - Checker backend: isolated standard-library replay with no producer or public
   contract imports
 
@@ -65,13 +67,25 @@ zero dimensions for reuse in chain-complex witnesses. Every output integer is
 bounded to 32,768 digits. Complete request validation occurs before computation
 or artifact writes.
 
-The implementation uses exact elementary row and column operations and has no
-optional-provider availability gate. Python-FLINT 0.9.0's Python API exposes a
-diagonal-only `snf()` operation even though current FLINT C documentation also
-describes `fmpz_mat_snf_transform`; SymPy's public normal-form API likewise
-returns the Smith form without full transformations. Those systems are
-independent test oracles for invariant factors, not hidden runtime
-dependencies of this producer.
+The producer delegates the Smith decomposition to SymPy's
+`smith_normal_decomp` over `ZZ` and has no optional-provider availability gate.
+The canonical diagonal \(D\) and the invariant factors \(d_1,\ldots,d_r\) are
+mathematical invariants: they are determined by the determinantal divisors of
+\(A\) and do not depend on the backend. The unimodular transformations
+\(U\) and \(V\), and every representative derived from them (kernel bases,
+cycle coordinates, generator coordinates, bounding chains), are deterministic
+for the pinned SymPy version but are **not** canonical: a different Smith
+backend or SymPy release may produce different \(U,V\) that satisfy the same
+relation \(D=UAV\). Compatibility is therefore semantic—the relation, both
+unimodular determinants, and the positive divisibility diagonal—rather than
+byte-identical transformations. The producer fail-closed checks verify all
+three before returning, and the independent checker replays them without
+calling the producer.
+
+Python-FLINT 0.9.0's Python API exposes a diagonal-only `snf()` operation
+even though current FLINT C documentation also describes
+`fmpz_mat_snf_transform`. FLINT and SymPy are independent test oracles for
+invariant factors, not hidden runtime dependencies of the checker.
 
 ## Integral simplicial homology
 
@@ -189,8 +203,10 @@ and reconsider the bounds only with new artifact-size and runtime evidence.
 This family does not provide persistent homology, a homology ring, cup
 products, canonical generators independent of the declared simplex
 orientation, manifold recognition, a preferred proof strategy, or conclusions
-beyond the supplied bounded complex. It does not make optional FLINT or SymPy
-providers authoritative checkers.
+beyond the supplied bounded complex. SymPy is the producer backend for the
+Smith decomposition; it is not an authoritative checker. The independent
+checker uses only standard-library integer arithmetic and never imports the
+producer or SymPy.
 
 ## Primary references
 

src/jacobian/contracts/certified_snf.py

@@ -143,8 +143,8 @@ class CertifiedSmithNormalFormResult(ContractModel):
     certificate: SmithNormalFormCertificate
     exactness: Literal["EXACT_INTEGER"] = "EXACT_INTEGER"
     determinism: Literal["DETERMINISTIC"] = "DETERMINISTIC"
-    backend: Literal["jacobian-stdlib-elementary-operations"] = (
-        "jacobian-stdlib-elementary-operations"
+    backend: Literal["jacobian-sympy-smith-normal-decomposition"] = (
+        "jacobian-sympy-smith-normal-decomposition"
     )
     backend_version: Literal["1"] = "1"
     completeness: Literal["FULL_MATRIX_TRANSFORMATIONS"] = "FULL_MATRIX_TRANSFORMATIONS"

src/jacobian/contracts/probability.py

@@ -21,7 +21,7 @@
 MAX_GAUSSIAN_POLYNOMIAL_TERMS = 16
 MAX_GAUSSIAN_TERM_DEGREE = 8
 MAX_GAUSSIAN_MOMENT_ORDER = 16
-MAX_GAUSSIAN_EXPANSION_PATHS = 4096
+MAX_GAUSSIAN_EXPANSION_PATHS = 65536
 MAX_GAUSSIAN_RESULT_RATIONAL_DIGITS = 4096
 MAX_GRAPH_RELIABILITY_VERTICES = 16
 MAX_GRAPH_RELIABILITY_EDGES = 12

src/jacobian/domains/_certified_snf.py

@@ -1,4 +1,12 @@
-"""Bounded elementary-operation Smith reduction used by trusted producers."""
+"""Smith normal decomposition used by trusted producers.
+
+The canonical Smith diagonal and divisibility chain are mathematical
+invariants.  The unimodular transformations ``U`` and ``V`` and every
+representative derived from them are deterministic for the pinned SymPy
+version but are **not** byte-identical across backend versions: compatibility
+is semantic (``D = U A V``, unimodularity, canonical diagonal) rather than
+representational.
+"""
 
 from __future__ import annotations
 
@@ -173,42 +181,24 @@ def inverse_unimodular(matrix: Matrix) -> Matrix:
     return [row[size:] for row in augmented]
 
 
-def _swap_rows(matrix: Matrix, left: int, right: int) -> None:
-    matrix[left], matrix[right] = matrix[right], matrix[left]
-
-
-def _swap_columns(matrix: Matrix, left: int, right: int) -> None:
-    for row in matrix:
-        row[left], row[right] = row[right], row[left]
-
-
-def _add_row_multiple(
-    matrix: Matrix, target: int, source: int, multiplier: int
-) -> None:
-    matrix[target] = [
-        value + multiplier * other
-        for value, other in zip(matrix[target], matrix[source], strict=True)
-    ]
-
-
-def _add_column_multiple(
-    matrix: Matrix, target: int, source: int, multiplier: int
-) -> None:
-    for row in matrix:
-        row[target] += multiplier * row[source]
-
-
-def _negate_row(matrix: Matrix, row: int) -> None:
-    matrix[row] = [-value for value in matrix[row]]
-
-
 def smith_reduce(
     source: Matrix,
     *,
     row_count: int | None = None,
     column_count: int | None = None,
 ) -> SmithReduction:
-    """Return a canonical Smith diagonal and explicit elementary transformations."""
+    """Return a canonical Smith diagonal and explicit unimodular transformations.
+
+    Delegates the decomposition to SymPy's ``smith_normal_decomp`` over ``ZZ``
+    and converts the result to native integer matrices.  The canonical diagonal
+    and invariant factors are mathematical invariants; ``U`` and ``V`` are
+    deterministic for the pinned SymPy version but may differ from other
+    backends.  Fail-closed checks verify ``D = U A V``, unimodularity, and the
+    positive divisibility chain before returning.
+    """
+
+    import sympy
+    from sympy.matrices.normalforms import smith_normal_decomp
 
     rows = len(source) if row_count is None else row_count
     columns = (
@@ -217,99 +207,51 @@ def smith_reduce(
     if len(source) != rows or any(len(row) != columns for row in source):
         raise ValueError("source entries do not match the declared matrix shape")
     original = [row[:] for row in source]
-    matrix = [row[:] for row in source]
-    left = identity_matrix(rows)
-    right = identity_matrix(columns)
-    diagonal_count = min(rows, columns)
 
-    for pivot in range(diagonal_count):
-        selected = min(
-            (
-                (abs(matrix[row][column]), row, column)
-                for row in range(pivot, rows)
-                for column in range(pivot, columns)
-                if matrix[row][column] != 0
-            ),
-            default=None,
-        )
-        if selected is None:
-            break
-        _, selected_row, selected_column = selected
-        if selected_row != pivot:
-            _swap_rows(matrix, pivot, selected_row)
-            _swap_rows(left, pivot, selected_row)
-        if selected_column != pivot:
-            _swap_columns(matrix, pivot, selected_column)
-            _swap_columns(right, pivot, selected_column)
-
-        while True:
-            changed = False
-            for row in range(pivot + 1, rows):
-                while matrix[row][pivot] != 0:
-                    quotient = matrix[row][pivot] // matrix[pivot][pivot]
-                    _add_row_multiple(matrix, row, pivot, -quotient)
-                    _add_row_multiple(left, row, pivot, -quotient)
-                    if matrix[row][pivot] != 0 and abs(matrix[row][pivot]) < abs(
-                        matrix[pivot][pivot]
-                    ):
-                        _swap_rows(matrix, row, pivot)
-                        _swap_rows(left, row, pivot)
-                    changed = True
-            for column in range(pivot + 1, columns):
-                while matrix[pivot][column] != 0:
-                    quotient = matrix[pivot][column] // matrix[pivot][pivot]
-                    _add_column_multiple(matrix, column, pivot, -quotient)
-                    _add_column_multiple(right, column, pivot, -quotient)
-                    if matrix[pivot][column] != 0 and abs(matrix[pivot][column]) < abs(
-                        matrix[pivot][pivot]
-                    ):
-                        _swap_columns(matrix, column, pivot)
-                        _swap_columns(right, column, pivot)
-                    changed = True
-            offender = next(
-                (
-                    (row, column)
-                    for row in range(pivot + 1, rows)
-                    for column in range(pivot + 1, columns)
-                    if matrix[row][column] % matrix[pivot][pivot]
-                ),
-                None,
-            )
-            if offender is not None:
-                _add_row_multiple(matrix, pivot, offender[0], 1)
-                _add_row_multiple(left, pivot, offender[0], 1)
-                changed = True
-                continue
-            if not changed:
-                break
-            if all(matrix[row][pivot] == 0 for row in range(pivot + 1, rows)) and all(
-                matrix[pivot][column] == 0 for column in range(pivot + 1, columns)
-            ):
-                break
-        if matrix[pivot][pivot] < 0:
-            _negate_row(matrix, pivot)
-            _negate_row(left, pivot)
+    # smith_normal_decomp accepts a plain SymPy Matrix and handles every shape,
+    # including 0xm and nx0 matrices, returning identity transformations for
+    # the empty side.
+    if rows and columns:
+        sympy_source = sympy.Matrix([[int(value) for value in row] for row in original])
+    else:
+        sympy_source = sympy.Matrix(rows, columns, [])
+
+    diagonal, left, right = smith_normal_decomp(sympy_source, domain=sympy.ZZ)
 
+    diagonal_matrix = [
+        [int(diagonal[row, column]) for column in range(columns)] for row in range(rows)
+    ]
+    left_matrix = [
+        [int(left[row, column]) for column in range(rows)] for row in range(rows)
+    ]
+    right_matrix = [
+        [int(right[row, column]) for column in range(columns)] for row in range(columns)
+    ]
+
+    diagonal_count = min(rows, columns)
     factors = tuple(
-        matrix[index][index]
+        diagonal_matrix[index][index]
         for index in range(diagonal_count)
-        if matrix[index][index] != 0
+        if diagonal_matrix[index][index] != 0
     )
     if any(value <= 0 for value in factors) or any(
         right_factor % left_factor for left_factor, right_factor in pairwise(factors)
     ):
         raise ArithmeticError("Smith reduction did not produce a canonical diagonal")
-    if matrix_multiply(matrix_multiply(left, original), right) != matrix:
+    if (
+        matrix_multiply(matrix_multiply(left_matrix, original), right_matrix)
+        != diagonal_matrix
+    ):
         raise ArithmeticError("Smith transformations do not bind the source")
-    left_determinant = determinant(left)
-    right_determinant = determinant(right)
+    left_determinant = determinant(left_matrix)
+    right_determinant = determinant(right_matrix)
     if abs(left_determinant) != 1 or abs(right_determinant) != 1:
         raise ArithmeticError("Smith transformations are not unimodular")
     return SmithReduction(
         source=original,
-        diagonal=matrix,
-        left=left,
-        right=right,
+        diagonal=diagonal_matrix,
+        left=left_matrix,
+        right=right_matrix,
         rank=len(factors),
         invariant_factors=factors,
         left_determinant=left_determinant,

src/jacobian/domains/certified_snf/bundle.py

@@ -4,7 +4,7 @@
 from jacobian.domains.certified_snf.checkers import CERTIFIED_SNF_EXACT_REPLAY_CHECKERS
 from jacobian.domains.certified_snf.operations import CERTIFIED_SNF_CAPABILITIES
 from jacobian.operations import DomainBundle, DomainDiagnostics, DomainSemantics
-from jacobian.provider_runtime import known_provider_runtime
+from jacobian.provider_runtime import SYMPY_VERSION, known_provider_runtime
 
 
 def build_certified_snf_bundle() -> DomainBundle:
@@ -31,17 +31,16 @@ def build_certified_snf_bundle() -> DomainBundle:
             },
         ),
         provider_runtime=known_provider_runtime(
-            "jacobian.certified-snf",
+            "jacobian.sympy",
             features=(
                 "exact-integer",
-                "elementary-row-operations",
-                "elementary-column-operations",
+                "sympy-smith-normal-decomposition",
                 "left-unimodular-transformation",
                 "right-unimodular-transformation",
                 "smith-divisibility-chain",
             ),
         ),
-        backend_version="jacobian.certified-snf/1",
+        backend_version=SYMPY_VERSION,
         capabilities=CERTIFIED_SNF_CAPABILITIES,
         diagnostics=DomainDiagnostics(
             invalid_request=CapabilityDiagnostic(
@@ -56,12 +55,15 @@ def build_certified_snf_bundle() -> DomainBundle:
         ),
         scope_description="the complete supplied bounded integer matrix",
         completeness_basis=(
-            "elementary row and column operations produced both full basis changes "
-            "and the complete canonical Smith diagonal"
+            "SymPy's smith_normal_decomp over ZZ produced both full unimodular "
+            "basis changes and the complete canonical Smith diagonal; the producer "
+            "fail-closed checks verified D = U A V, both unimodular determinants, "
+            "and the positive divisibility chain"
         ),
         assurance_basis=(
-            "exact integer computation capped at COMPUTED; independent certificate "
-            "replay is available through matrix.normal_form.smith.certified.verify"
+            "exact SymPy smith_normal_decomp capped at COMPUTED; independent "
+            "certificate replay is available through "
+            "matrix.normal_form.smith.certified.verify"
         ),
         checker_declarations=CERTIFIED_SNF_EXACT_REPLAY_CHECKERS,
     )

src/jacobian/domains/certified_snf/checkers.py

@@ -12,7 +12,7 @@
         entrypoint_module="jacobian_checkers.certified_snf",
         replay_method="independent Smith transformation-certificate replay",
         reason=(
-            "operator-authorized standard-library checker validates D=UAV, both "
+            "operator-authorized independent checker validates D=UAV, both "
             "unimodular determinants, and the complete canonical divisibility chain"
         ),
         verification_capability_id="matrix.normal_form.smith.certified.verify",

src/jacobian/domains/certified_snf/operations.py

@@ -63,7 +63,7 @@ def _certified_smith(
                 },
             ),
         ),
-        version="3",
+        version="4",
     ),
 )