perf(row-reduce): optimize prepared helpers and complete Phase 4 evidence

.github/workflows/ci.yml

@@ -72,6 +72,7 @@ jobs:
           python3 -m unittest scripts/bench/test_primality_negative_sweep.py
           python3 -m unittest scripts/bench/test_primality_mathlib_proof_sweep.py
           python3 -m unittest scripts/oracle/test_flint_mpoly_bench.py
+          python3 -m unittest scripts/oracle/test_flint_fmpq_bench.py
           python3 -m unittest scripts/oracle/test_singular_mpoly_bench.py
           python3 -m unittest scripts/oracle/test_rcf_flint_bench.py
           python3 scripts/ci/check_benches_mathlib_free.py
@@ -172,7 +173,7 @@ jobs:
       - name: Define the hex-dev cache + build target sets
         run: |
           echo "HEX_LIB_TARGETS=HexBasic HexTruncatedSeries HexTruncatedSeriesMathlib HexArith HexPoly HexPolyFast HexMvPoly HexModArith HexGF2 HexPolyZ HexRoots HexResultant HexInterval HexIntervalExperiment HexIntervalMathlib HexIntervalMathlibExperiment HexIntervalReplayProbe HexIntervalMathlibReplayProbe HexRealRootsMathlibReplayProbe HexRCFProofProbe HexPolyFp HexGFqRing HexGFqField HexBerlekamp HexHensel HexConway HexGFq HexPrimality HexIntFactor HexPrimalityKernelProbe HexPrimalityElabProbe HexPrimalityMathlibProofProbe HexIntFactorKernelProbe HexMvGcdKernelProbe HexBerlekampZassenhaus HexRealRoots HexMatrix HexRowReduce HexDeterminant HexBareiss HexHermite HexSmith HexCharPoly HexMinPoly HexPolySmith HexGramSchmidt HexLLL HexMatrixMathlib HexHermiteMathlib HexSmithMathlib HexSmithTests HexCharPolyMathlib HexMinPolyMathlib HexPolySmithMathlib HexGramSchmidtMathlib HexLLLMathlib HexBerlekampZassenhausMathlib HexBerlekampZassenhausMathlibProofProbe HexBerlekampMathlibProofProbe HexPrimalityMathlib HexIntFactorMathlib HexPolyFpMathlib HexFactorizationModules HexRealRootsMathlib HexGF2Mathlib HexGFqMathlib HexMvPolyMathlib HexMvPolyMathlibProofProbe HexRCF HexRCFTests HexRootsMathlib HexResultantMathlib HexNumberField HexNumberFieldMathlib HexNumberFieldTower HexNumberFieldTowerMathlib HexReleaseTests HexReleaseExamples HexAggregateCheck" >> "$GITHUB_ENV"
-          echo "HEX_EXE_TARGETS=hextruncatedseries_bench hexarith_bench hexpoly_bench hexpolyfast_bench hexpolysmith_bench hexmvpoly_bench hexmvgcd_bench hexsparsepoly_bench hexpolyz_bench hexpolyfp_bench hexmodarith_bench hexmodular_bench hexgf2_bench hexgfqring_bench hexgfqfield_bench hexgfq_bench hexhensel_bench hexprimality_bench hexprimality_policy_probe hexprimality_fuel_probe hexintfactor_bench hexberlekamp_bench hexbz_bench hexconway_bench hexmatrix_bench hexstrassen_compare hexdeterminant_bench hexbareiss_bench hexcharpoly_bench hexminpoly_bench hexgramschmidt_bench hexhermite_bench hexsmith_bench hexrealroots_bench hexrcf_bench hexroots_bench hexresultant_bench hexnumberfield_bench hexnumberfieldtower_bench hexinterval_decision_bench hexroots_demo hex_interval_representation_spike hex_interval_center_spike hex_interval_scale_spike hex_interval_scheduler_spike hex_interval_policy_frontier_spike hex_arith_floor hexlll_bench hexlll_gram_bench hexlll_external_reduction hexbz_factor_service" >> "$GITHUB_ENV"
+          echo "HEX_EXE_TARGETS=hextruncatedseries_bench hexarith_bench hexpoly_bench hexpolyfast_bench hexpolysmith_bench hexmvpoly_bench hexmvgcd_bench hexsparsepoly_bench hexpolyz_bench hexpolyfp_bench hexmodarith_bench hexmodular_bench hexgf2_bench hexgfqring_bench hexgfqfield_bench hexgfq_bench hexhensel_bench hexprimality_bench hexprimality_policy_probe hexprimality_fuel_probe hexintfactor_bench hexberlekamp_bench hexbz_bench hexconway_bench hexmatrix_bench hexstrassen_compare hexrowreduce_bench hexdeterminant_bench hexbareiss_bench hexcharpoly_bench hexminpoly_bench hexgramschmidt_bench hexhermite_bench hexsmith_bench hexrealroots_bench hexrcf_bench hexroots_bench hexresultant_bench hexnumberfield_bench hexnumberfieldtower_bench hexinterval_decision_bench hexroots_demo hex_interval_representation_spike hex_interval_center_spike hex_interval_scale_spike hex_interval_scheduler_spike hex_interval_policy_frontier_spike hex_arith_floor hexlll_bench hexlll_gram_bench hexlll_external_reduction hexbz_factor_service" >> "$GITHUB_ENV"
       # Shared build. The libraries, bench exes, conformance #guard drivers, and
       # emit-fixture exes are all elaborated here so the two verification tails
       # below only *run* things, never rebuild them -- which is what lets the
@@ -242,7 +243,7 @@ jobs:
             hexgf2_bench hexgfqring_bench hexgfqfield_bench hexgfq_bench \
             hexhensel_bench hexprimality_bench hexintfactor_bench \
             hexberlekamp_bench hexbz_bench \
-            hexconway_bench hexdeterminant_bench hexmatrix_bench \
+            hexconway_bench hexdeterminant_bench hexmatrix_bench hexrowreduce_bench \
             hexhermite_bench hexsmith_bench \
             hexcharpoly_bench hexminpoly_bench \
             hexgramschmidt_bench hexlll_gram_bench \

HexBerlekamp/PackedKernel.lean

@@ -875,6 +875,7 @@ theorem nullspaceArray_eq {q : UInt32} (hq : q.toNat = p)
     by_cases hjf : (⟨j, hj⟩ : Fin m) = hE.toIsEchelonForm.freeCols.get ⟨k, hklt⟩
     · rw [ite_eq_left hjf, dite_eq_left hjf]
     · rw [ite_eq_right hjf, dite_eq_right hjf, pivotRowOf_eq (Matrix.rowReduce E) ⟨j, hj⟩]
+      rw [Matrix.IsRowReduced.pivotRows_get]
       cases hpi : Matrix.IsRowReduced.pivotIndex? (Matrix.rowReduce E) (⟨j, hj⟩ : Fin m) with
       | none => rfl
       | some i =>

HexRowReduce/Nullspace.lean

@@ -116,16 +116,28 @@ private theorem pivotIndex?_free_none (E : IsEchelonForm M D) (k : Fin (m - D.ra
   intro i
   exact E.pivotCols_disjoint_freeCols i k
 
+omit [Mul R] [Add R] [OfNat R 0] [OfNat R 1] in
+/-- Reading the materialized pivot-row table agrees with `pivotIndex?`. -/
+@[simp] theorem pivotRows_get (j : Fin m) :
+    (Vector.ofFn fun j => pivotIndex? D j).get j = pivotIndex? D j := by
+  change (Vector.ofFn fun j => pivotIndex? D j)[j] = pivotIndex? D j
+  simp
+
 /-- Nullspace basis vectors assembled as columns indexed by the free variables. -/
 @[expose]
 def nullspaceMatrix [Lean.Grind.Ring R] (E : IsRowReduced M D) :
     Matrix R m (m - D.rank) :=
   let freeCols := E.toIsEchelonForm.freeCols
+  -- Materialize the column-to-pivot-row lookup once.  `Matrix.ofFn` then
+  -- performs constant-time vector reads instead of scanning `pivotCols` for
+  -- every output entry.
+  let pivotRows : Vector (Option (Fin D.rank)) m :=
+    Vector.ofFn fun j => pivotIndex? D j
   Matrix.ofFn fun j k =>
     if hFree : j = freeCols.get k then
       1
     else
-      match pivotIndex? D j with
+      match pivotRows.get j with
       | some i =>
           -D.echelon[(IsEchelonForm.pivotRow E.toIsEchelonForm i, freeCols.get k)]
       | none => 0
@@ -145,7 +157,7 @@ def nullspaceMatrix [Lean.Grind.Ring R] (E : IsRowReduced M D) :
   have hne : E.toIsEchelonForm.freeCols.get l ≠ E.toIsEchelonForm.freeCols.get k := by
     intro h
     exact hkl ((E.toIsEchelonForm.freeCols_injective h).symm)
-  simp [hne, pivotIndex?_free_none E.toIsEchelonForm l]
+  simp [hne, pivotRows_get, pivotIndex?_free_none E.toIsEchelonForm l]
 
 /-- In a pivot-column row, a nullspace-matrix entry is the negative RREF entry in
 the matching pivot row and free column. -/
@@ -154,7 +166,7 @@ the matching pivot row and free column. -/
     E.nullspaceMatrix[D.pivotCols.get i][k] =
       -(D.echelon[(IsEchelonForm.pivotRow E.toIsEchelonForm i)][E.toIsEchelonForm.freeCols.get k]) := by
   simp only [nullspaceMatrix, getElem_ofFn]
-  simp [E.toIsEchelonForm.pivotCols_disjoint_freeCols i k,
+  simp [pivotRows_get, E.toIsEchelonForm.pivotCols_disjoint_freeCols i k,
     pivotIndex?_pivot E.toIsEchelonForm i]
 
 /-- The individual nullspace basis vectors.

HexRowReduce/RowEchelon/Contracts.lean

@@ -128,6 +128,127 @@ Mathlib-free layer). -/
 def freeColsList (_E : IsEchelonForm M D) : List (Fin m) :=
   (List.finRange m).filter fun j => j ∉ D.pivotCols.toList
 
+/-- Linear merge used by the compiled implementation of `freeColsList`.
+
+The second list is a sorted sublist of the first.  Matching heads are
+discarded; otherwise the first head belongs to the complement.  Both recursive
+calls consume a column, and a matching call also consumes a pivot, so the loop
+is linear in the two input lengths. -/
+@[expose]
+def sortedComplement [DecidableEq α] : List α → List α → List α
+  | xs, [] => xs
+  | [], _ => []
+  | x :: xs, y :: ys =>
+      if x = y then
+        sortedComplement xs ys
+      else
+        x :: sortedComplement xs (y :: ys)
+
+private theorem sortedComplement_eq_filter [DecidableEq α]
+    {pivots cols : List α} (hsub : pivots.Sublist cols) (hnodup : cols.Nodup) :
+    sortedComplement cols pivots = cols.filter fun x => x ∉ pivots := by
+  induction hsub with
+  | slnil => simp [sortedComplement]
+  | @cons pivots cols x hsub ih =>
+      have hnot : x ∉ pivots := by
+        intro hx
+        exact (List.nodup_cons.mp hnodup).1 (hsub.subset hx)
+      cases pivots with
+      | nil =>
+          rw [sortedComplement]
+          symm
+          apply List.filter_eq_self.mpr
+          intro z hz
+          simp
+      | cons y ys =>
+          have hxy : x ≠ y := by
+            intro h
+            subst y
+            exact hnot (by simp)
+          rw [sortedComplement, if_neg hxy]
+          rw [ih (List.nodup_cons.mp hnodup).2]
+          rw [List.filter_cons]
+          have hp : (decide (x ∉ y :: ys) : Bool) = true := decide_eq_true hnot
+          rw [hp]
+          simp
+  | @cons_cons pivots cols x hsub ih =>
+      rw [sortedComplement, if_pos rfl, List.filter_cons_of_neg (by simp)]
+      rw [ih (List.nodup_cons.mp hnodup).2]
+      apply List.filter_congr
+      intro z hz
+      have hzx : z ≠ x := by
+        intro h
+        subst z
+        exact (List.nodup_cons.mp hnodup).1 hz
+      simp [hzx]
+
+private theorem pairwise_perm_eq {xs ys : List (Fin m)}
+    (hxs : List.Pairwise (fun a b : Fin m => a < b) xs)
+    (hys : List.Pairwise (fun a b : Fin m => a < b) ys)
+    (hperm : xs.Perm ys) : xs = ys := by
+  induction xs generalizing ys with
+  | nil => exact hperm.symm.eq_nil.symm
+  | cons x xs ih =>
+      cases ys with
+      | nil =>
+          have hlen := hperm.length_eq
+          simp at hlen
+      | cons y ys =>
+          have hxmem : x ∈ y :: ys := hperm.subset (by simp)
+          have hymem : y ∈ x :: xs := hperm.symm.subset (by simp)
+          have hxy : x = y := by
+            rcases List.mem_cons.mp hxmem with h | hx
+            · exact h
+            · rcases List.mem_cons.mp hymem with h | hy
+              · exact h.symm
+              · have hyx : y < x := (List.pairwise_cons.mp hys).1 x hx
+                have hxy : x < y := (List.pairwise_cons.mp hxs).1 y hy
+                omega
+          subst y
+          simp only [List.cons.injEq, true_and]
+          exact ih (List.pairwise_cons.mp hxs).2 (List.pairwise_cons.mp hys).2
+            hperm.cons_inv
+
+/-- Runtime implementation of `freeColsList`, merging the sorted pivot columns
+with the complete sorted column range. -/
+@[expose]
+def freeColsListImpl (_E : IsEchelonForm M D) : List (Fin m) :=
+  sortedComplement (List.finRange m) D.pivotCols.toList
+
+private theorem pivotCols_sublist_finRange (E : IsEchelonForm M D) :
+    D.pivotCols.toList.Sublist (List.finRange m) := by
+  let p : Fin m → Bool := fun j => decide (j ∈ D.pivotCols.toList)
+  have hfilterPairs : List.Pairwise (fun a b : Fin m => a < b)
+      ((List.finRange m).filter p) := by
+    exact List.Pairwise.filter p (List.pairwise_lt_finRange m)
+  have hfilterNodup : ((List.finRange m).filter p).Nodup := by
+    rw [List.nodup_iff_pairwise_ne]
+    exact hfilterPairs.imp (fun hlt heq => by subst heq; omega)
+  have hperm : D.pivotCols.toList.Perm ((List.finRange m).filter p) := by
+    rw [List.perm_ext_iff_of_nodup E.pivotCols_nodup hfilterNodup]
+    intro a
+    constructor
+    · intro ha
+      rw [List.mem_filter]
+      exact ⟨List.mem_finRange a, show p a = true from by exact decide_eq_true ha⟩
+    · intro ha
+      rw [List.mem_filter] at ha
+      exact of_decide_eq_true ha.2
+  have heq : D.pivotCols.toList = (List.finRange m).filter p :=
+    pairwise_perm_eq E.pivotCols_pairwise hfilterPairs hperm
+  rw [heq]
+  exact List.filter_sublist
+
+/-- Register the linear merge as the compiled implementation of
+`freeColsList`.  The proof keeps the specification and all theorem-facing
+unfolding unchanged. -/
+@[csimp]
+theorem freeColsList_eq_impl : @freeColsList = @freeColsListImpl := by
+  funext R n m _ _ _ _ M D E
+  unfold freeColsList freeColsListImpl
+  exact (sortedComplement_eq_filter (E.pivotCols_sublist_finRange)
+    (List.nodup_finRange m)).symm
+
 /-- The number of free columns is the ambient column count minus the rank. -/
 theorem freeColsList_length (E : IsEchelonForm M D) :
     E.freeColsList.length = m - D.rank := by

HexRowReduce/SPEC/hex-row-reduce.md

@@ -106,11 +106,46 @@ with `colPartition` (free columns telescope to `v[freeCols[l]]`; pivot columns
 follow from `pivot_one` / `above_pivot_zero` / `below_pivot_zero` / `zero_row`);
 package into `E.nullspaceMatrix * c = v`.
 
+## Complexity and benchmark contract
+
+For an `n × m` matrix of rank `r`, Gauss--Jordan reduction performs at most
+`r` pivot stages and visits `O(n(n + m))` entries per stage across the echelon
+and transform matrices.  The arithmetic-operation bound is therefore
+`O(rn(n + m))`, cubic on the fixed-aspect square benchmark families.  Exact
+rational bit cost additionally depends on numerator and denominator growth.
+
+The compiled implementation uses a proved linear sorted merge for the
+free-column complement.  Prepared nullspace construction materializes the
+column-to-pivot-row lookup once and then writes `m(m - r)` output entries with
+constant-time lookup per entry; its fixed-aspect bound is quadratic.
+
+`bench/HexRowReduce/Bench.lean` gives all 12 advertised executable operations
+direct mode-1 coverage:
+
+| Operations | Prepared state | Model |
+| --- | --- | --- |
+| `Matrix.rowReduce`, `Matrix.rowReduce_rank` | dense input | `n³` |
+| `Matrix.spanCoeffs`, `Matrix.spanContains` | dense input | `n³` |
+| `IsEchelonForm.spanCoeffs`, `IsEchelonForm.spanContains` | proved RREF | `n²` |
+| `IsEchelonForm.echelonCoeffs`, `IsEchelonForm.freeCols` | proved RREF | `n` |
+| `Matrix.nullspaceBasisMatrix`, `Matrix.nullspace` | deficient input | `n³` |
+| `IsRowReduced.nullspaceMatrix`, `IsRowReduced.nullspace` | proved RREF | `n²` |
+
+The dense family is `I + J`, so every pivot fires while coefficient heights
+remain controlled.  The deficient family has rank and nullity `n / 2`.
+Preparation is outside the timed region and result forcing is inside it.
+
 ## External comparators
 
-The `rank`, `rowReduce`, and `nullspace` operations are cross-checked for correctness
-against python-flint's `fmpz_mat` / `fmpq_mat` through the conformance oracle
-(`scripts/oracle/matrix_flint.py`, driven by `hexrowreduce_emit_fixtures`).
-There is no Phase-4 performance comparator: row reduction is an exact rational
-computation validated for correctness, not timed against an external tool. See
-`reports/hex-row-reduce-performance.md`.
+The identical constant-size rank result is compared informationally with
+python-flint's `fmpq_mat.rref()` through the shared persistent driver.  Both
+arms use the same dense `I + J` family; construction is cached during warmup,
+and each timed request returns only the integer rank.
+
+The remaining operations declare
+`no-comparable-surface-in-named-comparator`: Hex `rowReduce` returns the row
+transform that `fmpq_mat.rref()` omits; python-flint 0.9.0's `fmpq_mat` has no
+native nullspace callable; and span coefficients are transform-dependent
+witnesses.  A comparator-specific derived algorithm would not be the same
+callable surface.  Full RREF and nullspace correctness remain cross-checked by
+`scripts/oracle/matrix_flint.py`, driven by `hexrowreduce_emit_fixtures`.