feat(Schonhage): execute native queried thin matrix products

Complexity/Algorithms/MatrixMultiplication/Schonhage/BandBounds.lean

@@ -0,0 +1,102 @@
+/-
+Copyright (c) 2026 Erik Demaine. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Erik Demaine
+-/
+
+module
+-- upstream: candidate
+
+public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Layout
+import Mathlib.Data.Nat.Choose.Bounds
+import Complexity.Algorithms.MatrixMultiplication.Schonhage.Tiling
+
+/-!
+# Bounds for guarded band addressing
+
+## Main results
+
+* `bandSize_le_seven_pow`: a complete band fits within the tensor input dimension.
+* `band_position_bounds`: all padded row/column arithmetic stays below `N+7^L` before
+  the native guard. Row-major matrix addresses are computed only after that guard.
+-/
+
+@[expose] public section
+
+namespace Cslib.Algorithms.MatrixMultiplication.Schonhage
+
+private theorem mul_add_lt_mul {a b c d : ℕ} (ha : a < b) (hc : c < d) :
+    a * d + c < b * d := by
+  have h := Nat.mul_le_mul_right d (Nat.succ_le_of_lt ha)
+  rw [Nat.succ_mul] at h
+  omega
+
+/-- The entire band dimension is bounded by the flat tensor input dimension. -/
+theorem bandSize_le_seven_pow (L m : ℕ) : bandSize L m ≤ 7 ^ L := by
+  have hK : K0 L m ≤ 2 ^ L := (Nat.sqrt_le_self _).trans (Nat.choose_le_two_pow _ _)
+  have hN : N0 L m ≤ 3 ^ L := Nat.pow_le_pow_right (by decide) (Nat.sub_le _ _)
+  calc bandSize L m ≤ 2 ^ L * 3 ^ L := Nat.mul_le_mul hK hN
+    _ = 6 ^ L := by rw [← Nat.mul_pow]
+    _ ≤ 7 ^ L := Nat.pow_le_pow_left (by decide) _
+
+/-- Even the last padded band has small block arithmetic. This bound precedes the `r<N`
+guard; it does not assert that a padded row/column is a real matrix position. -/
+theorem band_position_bounds {L m N β q o : ℕ} (hβ : β < numBands L m N)
+    (hq : q < K0 L m) (ho : o < N0 L m) :
+    β * K0 L m < N + 7 ^ L ∧
+      β * K0 L m + q < N + 7 ^ L ∧
+      (β * K0 L m + q) * N0 L m < N + 7 ^ L ∧
+      (β * K0 L m + q) * N0 L m + o < N + 7 ^ L := by
+  have hb : β * K0 L m + q < numBands L m N * K0 L m := mul_add_lt_mul hβ hq
+  have hr : (β * K0 L m + q) * N0 L m + o < numBands L m N * K0 L m * N0 L m :=
+    mul_add_lt_mul hb ho
+  have hpad : numBands L m N * bandSize L m ≤ N + bandSize L m - 1 :=
+    Nat.div_mul_le_self _ _
+  have hsize := bandSize_le_seven_pow L m
+  have hblock : β * K0 L m + q ≤ (β * K0 L m + q) * N0 L m :=
+    Nat.le_mul_of_pos_right _ (N0_pos L m)
+  have heq : numBands L m N * K0 L m * N0 L m = numBands L m N * bandSize L m := by
+    rw [bandSize, Nat.mul_assoc]
+  rw [heq] at hr
+  omega
+
+/-- Nonempty matrices need at least one shared band. -/
+theorem numBands_pos {L m N : ℕ} (hm : m ≤ L) (hN : 0 < N) : 0 < numBands L m N := by
+  have hb : 0 < bandSize L m := Nat.mul_pos (K0_pos hm) (N0_pos L m)
+  change 0 < (N + bandSize L m - 1) / bandSize L m
+  exact Nat.div_pos (by omega) hb
+
+/-- A full band's original matrix entries fit within one encoding-sized block. This bound
+is valid for every `m≤L`, without the thin-regime specialization. -/
+theorem bandSize_mul_D_le_ten_pow {L m : ℕ} (hm : m ≤ L) :
+    bandSize L m * D m ≤ 10 ^ L := by
+  have hK : K0 L m ≤ 2 ^ L := (Nat.sqrt_le_self _).trans (Nat.choose_le_two_pow _ _)
+  have hN : N0 L m ≤ 4 ^ (L - m) := Nat.pow_le_pow_left (by decide) _
+  calc bandSize L m * D m ≤ (2 ^ L * 4 ^ (L - m)) * 4 ^ m :=
+      Nat.mul_le_mul_right _ (Nat.mul_le_mul hK hN)
+    _ = 8 ^ L := by
+      rw [Nat.mul_assoc, ← Nat.pow_add, Nat.sub_add_cancel hm, ← Nat.mul_pow]
+    _ ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _
+
+/-- All original matrix entries are dominated by the shared band catalogue size; this
+charges input copying/restoration without a separate thin-regime estimate. -/
+theorem matrix_entries_le_band_encodings {L m : ℕ} (hm : m ≤ L) (N : ℕ) :
+    N * D m ≤ numBands L m N * 10 ^ L := by
+  calc N * D m ≤ padN L m N * D m := Nat.mul_le_mul_right _ (le_padN hm N)
+    _ = numBands L m N * (bandSize L m * D m) := by rw [padN, Nat.mul_assoc]
+    _ ≤ numBands L m N * 10 ^ L := Nat.mul_le_mul_left _ (bandSize_mul_D_le_ten_pow hm)
+
+/-- For nonempty matrices, table initialization and the original matrix input are all
+dominated by a constant times the shared encoding entry count. -/
+theorem preparation_shape_le {L m N : ℕ} (hm : m ≤ L) (hN : 0 < N) :
+    1 + N * D m + 2 ^ L + 2 * numBands L m N * 10 ^ L ≤
+      5 * (numBands L m N * 10 ^ L) := by
+  have hinput := matrix_entries_le_band_encodings hm N
+  have hmask : 2 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _
+  have hpos : 1 ≤ 10 ^ L := Nat.one_le_pow L 10 (by decide)
+  have hbands : 10 ^ L ≤ numBands L m N * 10 ^ L :=
+    Nat.le_mul_of_pos_left _ (numBands_pos hm hN)
+  rw [Nat.mul_assoc]
+  omega
+
+end Cslib.Algorithms.MatrixMultiplication.Schonhage

Complexity/Algorithms/MatrixMultiplication/Schonhage/Batch.lean

@@ -105,7 +105,8 @@ private theorem exists_of_summand_ne_zero {L m : ℕ} (Q : Finset (Fin L)) (hQ :
     (r c : OuterStr L m) {u : LeftStr L} {v : RightStr L}
     (hγ : ∀ ℓ, coefficient (u ℓ) (v ℓ) (outStrOf Q hQ r c ℓ) ≠ 0) (hu : (innerSetL u).card = m)
     (hv : (innerSetR v).card = m) : ∃ π, u = leftStrOf Q hQ r π ∧ v = rightStrOf Q hQ π c := by
-  simp only [outStrOf, forall_level_iff Q hQ, glue_innerLevel, glue_outerLevel] at hγ
+  simp only [outStrOf, forall_iff_orderEmbOfFin_and_compl Q hQ, glue_orderEmbOfFin,
+    glue_orderEmbOfFin_compl] at hγ
   obtain ⟨hγin, hγout⟩ := hγ
   -- at a level ℓ ∈ Q we have w_ℓ = z₀ […] So (u_ℓ, v_ℓ) = (p_ij, q_ij) for some i, j.
   have hin : ∀ k, ∃ ij : Fin 2 × Fin 2,
@@ -117,11 +118,11 @@ private theorem exists_of_summand_ne_zero {L m : ℕ} (Q : Finset (Fin L)) (hQ :
   -- Thus u and v have an inner variable at every level of Q, and their inner sets are therefore
   -- exactly Q.  So u_ℓ and v_ℓ are outer at the levels outside Q.
   have huQ : innerSetL u = Q :=
-    filter_eq_of_forall_innerLevel Q hQ (fun k => by rw [(hπ k).1]; trivial) hu
+    filter_eq_of_forall_orderEmbOfFin Q hQ (fun k => by rw [(hπ k).1]; trivial) hu
   have hvQ : innerSetR v = Q :=
-    filter_eq_of_forall_innerLevel Q hQ (fun k => by rw [(hπ k).2]; trivial) hv
-  have huout := ((filter_eq_iff_forall_level Q hQ).mp huQ).2
-  have hvout := ((filter_eq_iff_forall_level Q hQ).mp hvQ).2
+    filter_eq_of_forall_orderEmbOfFin Q hQ (fun k => by rw [(hπ k).2]; trivial) hv
+  have huout := ((filter_eq_iff_orderEmbOfFin_and_compl Q hQ).mp huQ).2
+  have hvout := ((filter_eq_iff_orderEmbOfFin_and_compl Q hQ).mp hvQ).2
   -- at a level ℓ ∉ Q we have w_ℓ = z_ij for some i, j […] so (u_ℓ, v_ℓ) = (x_i, y_j).
   have hout : ∀ k, u ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k) = .x (r k) ∧ v
     ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k) = .y (c k) :=
@@ -135,8 +136,8 @@ private theorem prod_gamma_eq_one {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card
     ∏ ℓ, coefficient (leftStrOf Q hQ r π ℓ) (rightStrOf Q hQ π c ℓ) (outStrOf Q hQ r c ℓ) = 1 := by
   have hone : ∀ ℓ,
       coefficient (leftStrOf Q hQ r π ℓ) (rightStrOf Q hQ π c ℓ) (outStrOf Q hQ r c ℓ) = 1 := by
-    simp only [forall_level_iff Q hQ, leftStrOf, rightStrOf, outStrOf, glue_innerLevel,
-      glue_outerLevel]
+    simp only [forall_iff_orderEmbOfFin_and_compl Q hQ, leftStrOf, rightStrOf, outStrOf,
+      glue_orderEmbOfFin, glue_orderEmbOfFin_compl]
     exact ⟨fun _ => (coefficient_z0 _ _).trans (ite_eq_left ⟨_, _, rfl, rfl⟩),
       fun _ => (coefficient_x_y_z _ _ _ _).trans (ite_eq_left ⟨rfl, rfl⟩)⟩
   exact prod_eq_one fun ℓ _ => hone ℓ

Complexity/Algorithms/MatrixMultiplication/Schonhage/Levels.lean

@@ -187,35 +187,36 @@ theorem innerSetO_outStrOf (r c : OuterStr L m) : innerSetO (outStrOf Q hQ r c)
 theorem outerPartL_leftStrOf (r : OuterStr L m) (π : InnerStr m)
     (h : (innerSetL (leftStrOf Q hQ r π)).card = m) : outerPartL (leftStrOf Q hQ r π) h = r :=
   funext fun k => congrArg LeftVar.outerIndex
-    (glue_outerLevel_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k)
+    (glue_orderEmbOfFin_compl_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k)
 
 /-- The inner part of `leftStrOf Q hQ r π` is `π`. -/
 @[simp]
 theorem innerPartL_leftStrOf (r : OuterStr L m) (π : InnerStr m)
     (h : (innerSetL (leftStrOf Q hQ r π)).card = m) : innerPartL (leftStrOf Q hQ r π) h = π :=
   funext fun k => congrArg LeftVar.innerIndex
-    (glue_innerLevel_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k)
+    (glue_orderEmbOfFin_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k)
 
 /-- The inner part of `rightStrOf Q hQ π c` is `π`. -/
 @[simp]
 theorem innerPartR_rightStrOf (π : InnerStr m) (c : OuterStr L m)
     (h : (innerSetR (rightStrOf Q hQ π c)).card = m) : innerPartR (rightStrOf Q hQ π c) h = π :=
   funext fun k => congrArg RightVar.innerIndex
-    (glue_innerLevel_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k)
+    (glue_orderEmbOfFin_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k)
 
 /-- The outer part of `rightStrOf Q hQ π c` is `c`. -/
 @[simp]
 theorem outerPartR_rightStrOf (π : InnerStr m) (c : OuterStr L m)
     (h : (innerSetR (rightStrOf Q hQ π c)).card = m) : outerPartR (rightStrOf Q hQ π c) h = c :=
   funext fun k => congrArg RightVar.outerIndex
-    (glue_outerLevel_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k)
+    (glue_orderEmbOfFin_compl_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k)
 
 /-- The variable of `outStrOf Q hQ r c` at the `k`-th lowest level outside its inner set. -/
 private theorem outStrOf_outerLevel (r c : OuterStr L m)
     (h : (innerSetO (outStrOf Q hQ r c)).card = m) (k : Fin (L - m)) :
     outStrOf Q hQ r c (((innerSetO (outStrOf Q hQ r c))ᶜ.orderEmbOfFin (card_compl_of_card_eq
       (innerSetO (outStrOf Q hQ r c)) h)) k) = .z (r k) (c k) := by
-  rw [outerLevel_congr Q hQ (innerSetO_outStrOf Q hQ r c) h, outStrOf, glue_outerLevel]
+  simp only [innerSetO_outStrOf] at h ⊢
+  exact glue_orderEmbOfFin_compl Q hQ _ _ k
 
 /-- The row of `outStrOf Q hQ r c` is `r`. -/
 @[simp]
@@ -236,7 +237,8 @@ part. -/
 theorem leftStr_eq_leftStrOf {u : LeftStr L} (hu : innerSetL u = Q) (h : (innerSetL u).card = m) :
     u = leftStrOf Q hQ (outerPartL u h) (innerPartL u h) := by
   subst hu
-  refine glue_index _ h LeftVar.IsInner (inner := fun π : Fin 2 × Fin 2 => .p π.1 π.2) ?_ ?_ rfl
+  refine eq_glue_of_reconstruct _ h LeftVar.IsInner
+    (fromSubset := fun π : Fin 2 × Fin 2 => .p π.1 π.2) ?_ ?_ rfl
   · rintro (i | ⟨i, j⟩) hs
     exacts [hs.elim, rfl]
   · rintro (i | ⟨i, j⟩) hs
@@ -247,7 +249,8 @@ part. -/
 theorem rightStr_eq_rightStrOf {v : RightStr L} (hv : innerSetR v = Q)
     (h : (innerSetR v).card = m) : v = rightStrOf Q hQ (innerPartR v h) (outerPartR v h) := by
   subst hv
-  refine glue_index _ h RightVar.IsInner (inner := fun π : Fin 2 × Fin 2 => .q π.1 π.2) ?_ ?_ rfl
+  refine eq_glue_of_reconstruct _ h RightVar.IsInner
+    (fromSubset := fun π : Fin 2 × Fin 2 => .q π.1 π.2) ?_ ?_ rfl
   · rintro (j | ⟨i, j⟩) ht
     exacts [ht.elim, rfl]
   · rintro (j | ⟨i, j⟩) ht
@@ -257,9 +260,9 @@ theorem rightStr_eq_rightStrOf {v : RightStr L} (hv : innerSetR v = Q)
 theorem outStr_eq_outStrOf {w : OutStr L} (hw : innerSetO w = Q) (h : (innerSetO w).card = m) :
     w = outStrOf Q hQ (rowO w h) (colO w h) := by
   subst hw
-  refine glue_index _ h OutVar.IsInner (inner := fun _ : Unit => .z0)
-    (outer := fun rc : Fin 3 × Fin 3 => .z rc.1 rc.2) (innerIndex := fun _ => ())
-    (outerIndex := fun z => (z.rowIndex, z.colIndex)) ?_ ?_ rfl
+  refine eq_glue_of_reconstruct _ h OutVar.IsInner (fromSubset := fun _ : Unit => .z0)
+    (fromCompl := fun rc : Fin 3 × Fin 3 => .z rc.1 rc.2) (toSubset := fun _ => ())
+    (toCompl := fun z => (z.rowIndex, z.colIndex)) ?_ ?_ rfl
   · rintro (⟨i, j⟩ | _) hz
     exacts [hz.elim, rfl]
   · rintro (⟨i, j⟩ | _) hz