diff --git a/Complexity.lean b/Complexity.lean index ef65ccf388..6206aeb426 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -31,6 +31,7 @@ public import Complexity.Algorithms.Lean.VanEmdeBoas.Insert public import Complexity.Algorithms.Lean.VanEmdeBoas.Layout public import Complexity.Algorithms.Lean.VanEmdeBoas.Successor public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Alphabets +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.BandBounds public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Batch public import Complexity.Algorithms.MatrixMultiplication.Schonhage.CoordinateTransform public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Full @@ -38,10 +39,12 @@ public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Identity public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Layout public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Levels public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Magnitude +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.MaskTable public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Orders public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Parameters public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Pruned public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Sparse +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.SplitDigits public import Complexity.Algorithms.MatrixMultiplication.Schonhage.ThinProduct public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Tiling public import Complexity.Algorithms.Reachability.BreadthFirst @@ -112,6 +115,28 @@ public import Complexity.Algorithms.WordRAM.Karger.FastCut public import Complexity.Algorithms.WordRAM.Karger.Graph public import Complexity.Algorithms.WordRAM.KargerMinCut public import Complexity.Algorithms.WordRAM.KargerStein +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandAll +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCatalogue +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCollect +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCopy +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCost +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCount +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandEncode +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandExec +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandFlat +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandHeader +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandInput +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandMeaning +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandMemory +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandParameters +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPopulation +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPrepare +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandProgram +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandRead +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandRoot +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandSplit +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandTables +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandThin public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.Cell public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.Controller public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.Copy @@ -1759,6 +1784,7 @@ public import Complexity.Computability.WordRAM.Multiword.Workspace public import Complexity.Computability.WordRAM.Native public import Complexity.Computability.WordRAM.Notation public import Complexity.Computability.WordRAM.PackedCount +public import Complexity.Computability.WordRAM.Power public import Complexity.Computability.WordRAM.RAM public import Complexity.Computability.WordRAM.RAMInput public import Complexity.Computability.WordRAM.Randomized.Defs @@ -2221,6 +2247,7 @@ public import Complexity.Foundations.Data.Finset.Expect public import Complexity.Foundations.Data.Finset.Finpartition public import Complexity.Foundations.Data.Finset.Markov public import Complexity.Foundations.Data.Finset.OrderGlue +public import Complexity.Foundations.Data.Finset.OrderedFilter public import Complexity.Foundations.Data.Finset.Quantile public import Complexity.Foundations.Data.Finset.Quotient public import Complexity.Foundations.Data.Finset.Sampling @@ -2243,6 +2270,7 @@ public import Complexity.Foundations.Data.List.Lex public import Complexity.Foundations.Data.List.Lookup public import Complexity.Foundations.Data.List.OfFn public import Complexity.Foundations.Data.List.Range +public import Complexity.Foundations.Data.List.RowMajor public import Complexity.Foundations.Data.List.SelectDistinct public import Complexity.Foundations.Data.Matrix.Binary public import Complexity.Foundations.Data.Matrix.IntegerBits @@ -2278,6 +2306,7 @@ public import Complexity.Foundations.Data.Nat.Div public import Complexity.Foundations.Data.Nat.Factorial public import Complexity.Foundations.Data.Nat.IntermediateValue public import Complexity.Foundations.Data.Nat.Log +public import Complexity.Foundations.Data.Nat.OfDigits public import Complexity.Foundations.Data.Nat.PolynomialBound.Closure public import Complexity.Foundations.Data.Nat.PolynomialBound.Exponential public import Complexity.Foundations.Data.Nat.PolynomialBound.Log diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/BandBounds.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/BandBounds.lean new file mode 100644 index 0000000000..4fc53bbe64 --- /dev/null +++ b/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 finFunctionFinEquiv.symm a i = (1 : Fin 2) + +/-- The integer mask of a finite level set, with level zero least significant. -/ +def maskCode {L : ℕ} (Q : Finset (Fin L)) : Fin (2 ^ L) := + finFunctionFinEquiv fun i => if i ∈ Q then (1 : Fin 2) else 0 + +/-- Decoding the explicitly numbered mask recovers the original set. -/ +theorem maskSet_maskCode {L : ℕ} (Q : Finset (Fin L)) : maskSet (maskCode Q) = Q := by + ext i + simp [maskSet, maskCode] + +/-- Encoding a decoded bounded mask recovers the original integer. -/ +theorem maskCode_maskSet {L : ℕ} (a : Fin (2 ^ L)) : maskCode (maskSet a) = a := by + apply finFunctionFinEquiv.symm.injective + funext i + simp only [maskCode, Equiv.symm_apply_apply, maskSet, mem_filter, mem_univ, true_and] + have h := (finFunctionFinEquiv.symm a i).isLt + split <;> apply Fin.ext <;> simp_all only [Fin.ext_iff, Fin.val_one, Fin.val_zero] + omega + +/-- Binary-mask numbering is an explicit bijection, not a chosen enumeration. -/ +def maskEquiv (L : ℕ) : Fin (2 ^ L) ≃ Finset (Fin L) where + toFun := maskSet + invFun := maskCode + left_inv := maskCode_maskSet + right_inv := maskSet_maskCode + +/-- Number of one digits among exactly the first `L` binary digits. -/ +def population (L a : ℕ) : ℕ := ∑ i : Fin L, a / 2 ^ i.val % 2 + +/-- A zero-depth scan has population zero. -/ +theorem population_zero (a : ℕ) : population 0 a = 0 := by simp [population] + +/-- One native remainder/division step removes the lowest digit. -/ +theorem population_succ (L a : ℕ) : + population (L + 1) a = a % 2 + population L (a / 2) := by + simp only [population, Fin.sum_univ_succ, Fin.val_zero, pow_zero, Nat.div_one, + Fin.val_succ, pow_succ', ← Nat.div_div_eq_div_mul] + +/-- There are at most `L` one digits, even for an unbounded natural argument. -/ +theorem population_le (L a : ℕ) : population L a ≤ L := by + induction L generalizing a with + | zero => rw [population_zero] + | succ L ih => + rw [population_succ] + have hm := Nat.mod_lt a (by decide : 0 < 2) + have ht := ih (a / 2) + omega + +/-- The loop's digit sum is the cardinality of the represented level set. -/ +theorem population_eq_card {L : ℕ} (a : Fin (2 ^ L)) : + population L a.val = (maskSet a).card := by + rw [maskSet, card_eq_sum_ones, sum_filter] + apply sum_congr rfl + intro i _ + have h : (finFunctionFinEquiv.symm a i).val = a.val / 2 ^ i.val % 2 := rfl + have hb := (finFunctionFinEquiv.symm a i).isLt + split <;> simp_all only [Fin.ext_iff, Fin.val_one] + omega + +/-- All valid masks, in the ordinary linear order on bounded naturals. -/ +def validMasks (L m : ℕ) : Finset (Fin (2 ^ L)) := + univ.filter fun a => population L a.val = m + +/-- Every size-`m` subset occurs once in the bounded mask scan. -/ +theorem card_validMasks (L m : ℕ) : (validMasks L m).card = K L m := by + calc + (validMasks L m).card = ((univ : Finset (Fin L)).powersetCard m).card := by + apply card_bij (fun a _ => maskSet a) + · intro a ha + rw [mem_powersetCard_univ] + rw [← population_eq_card] + exact (mem_filter.mp ha).2 + · intro a _ b _ h + exact (maskEquiv L).injective h + · intro Q hQ + refine ⟨maskCode Q, ?_, maskSet_maskCode Q⟩ + simp only [validMasks, mem_filter, mem_univ, true_and, population_eq_card, + maskSet_maskCode] + exact mem_powersetCard_univ.mp hQ + _ = K L m := by rw [card_powersetCard, card_univ, Fintype.card_fin]; rfl + +/-- Number of valid masks strictly below the scan counter. -/ +def prefixCount (L m t : ℕ) : ℕ := ((range t).filter fun a => population L a = m).card + +/-- The scan begins without any accepted mask. -/ +theorem prefixCount_zero (L m : ℕ) : prefixCount L m 0 = 0 := by simp [prefixCount] + +/-- A scan accepts at most one mask in each iteration. -/ +theorem prefixCount_le (L m t : ℕ) : prefixCount L m t ≤ t := + (card_filter_le _ _).trans (by rw [card_range]) + +/-- The actual scalar counter update for one mask. -/ +theorem prefixCount_succ (L m t : ℕ) : + prefixCount L m (t + 1) = prefixCount L m t + if population L t = m then 1 else 0 := by + by_cases h : population L t = m <;> simp [prefixCount, range_add_one, filter_insert, h] + +/-- Extending a scan cannot decrease the accepted-mask counter. -/ +theorem prefixCount_mono (L m : ℕ) {a b : ℕ} (h : a ≤ b) : + prefixCount L m a ≤ prefixCount L m b := + card_le_card (filter_subset_filter _ (range_mono h)) + +/-- Earlier valid masks have strictly smaller accepted-mask ranks. -/ +theorem prefixCount_lt {L m a b : ℕ} (ha : population L a = m) (hab : a < b) : + prefixCount L m a < prefixCount L m b := by + have h := prefixCount_mono L m (show a + 1 ≤ b by omega) + rw [prefixCount_succ] at h + simp only [ha, ↓reduceIte] at h + omega + +/-- The finite complete scan has exactly the binomial count. -/ +theorem prefixCount_full (L m : ℕ) : prefixCount L m (2 ^ L) = K L m := by + rw [← card_validMasks] + apply card_bij (fun a ha => (⟨a, mem_range.mp (mem_filter.mp ha).1⟩ : Fin (2 ^ L))) + · intro a ha + exact mem_filter.mpr ⟨mem_univ _, (mem_filter.mp ha).2⟩ + · intro a _ b _ h + exact congrArg Fin.val h + · intro a ha + exact ⟨a.val, mem_filter.mpr ⟨mem_range.mpr a.isLt, (mem_filter.mp ha).2⟩, rfl⟩ + +private theorem card_below_orderEmb {α : Type*} [LinearOrder α] (s : Finset α) + {n : ℕ} (h : s.card = n) (i : Fin n) : + (s.filter fun a => a < s.orderEmbOfFin h i).card = i.val := by + rw [← card_range i.val] + symm + apply card_bij (fun j hj => s.orderEmbOfFin h + ⟨j, (mem_range.mp hj).trans i.isLt⟩) + · intro j hj + exact mem_filter.mpr ⟨s.orderEmbOfFin_mem h _, + (s.orderEmbOfFin h).strictMono (mem_range.mp hj)⟩ + · intro j hj k hk hjk + exact congrArg Fin.val ((s.orderEmbOfFin h).injective hjk) + · intro a ha + let j := (s.orderIsoOfFin h).symm ⟨a, (mem_filter.mp ha).1⟩ + have hj : s.orderEmbOfFin h j = a := + congrArg Subtype.val ((s.orderIsoOfFin h).apply_symm_apply ⟨a, (mem_filter.mp ha).1⟩) + have hji : j < i := (s.orderEmbOfFin h).lt_iff_lt.mp (hj ▸ (mem_filter.mp ha).2) + exact ⟨j.val, mem_range.mpr hji, hj⟩ + +/-- The scalar scan count equals the rank among the bounded valid masks. -/ +theorem prefixCount_eq_bounded (L m t : ℕ) (ht : t ≤ 2 ^ L) : + prefixCount L m t = ((validMasks L m).filter fun a => a.val < t).card := by + apply card_bij (fun a ha => + (⟨a, (mem_range.mp (mem_filter.mp ha).1).trans_le ht⟩ : Fin (2 ^ L))) + · intro a ha + exact mem_filter.mpr ⟨mem_filter.mpr ⟨mem_univ _, (mem_filter.mp ha).2⟩, + mem_range.mp (mem_filter.mp ha).1⟩ + · intro a _ b _ h + exact congrArg Fin.val h + · intro a ha + exact ⟨a.val, mem_filter.mpr ⟨mem_range.mpr (mem_filter.mp ha).2, + (mem_filter.mp (mem_filter.mp ha).1).2⟩, rfl⟩ + +/-- Row-major block-pair rank in the full list of size-`m` masks. -/ +def pairRank {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : Fin (K L m) := + Fin.castLE (by simpa only [K0] using Nat.sqrt_le (K L m)) (finProdFinEquiv gh) + +/-- The selected mask at a block pair is the corresponding mask in increasing order. -/ +def selectedMask {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : Fin (2 ^ L) := + (validMasks L m).orderEmbOfFin (card_validMasks L m) (pairRank gh) + +/-- The computed prefix counter at a selected mask is its row-major block-pair rank. -/ +theorem prefixCount_selectedMask {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : + prefixCount L m (selectedMask gh).val = (pairRank gh).val := by + rw [prefixCount_eq_bounded _ _ _ (selectedMask gh).isLt.le] + exact card_below_orderEmb (validMasks L m) (card_validMasks L m) (pairRank gh) + +/-- The dense inverse-table word: zero for an unselected mask, otherwise rank plus one. -/ +def lookup (L m a : ℕ) : ℕ := + if a < 2 ^ L ∧ population L a = m ∧ prefixCount L m a < K0 L m * K0 L m + then prefixCount L m a + 1 else 0 + +/-- Exactly the selected valid masks have nonzero inverse-table words. -/ +theorem lookup_ne_zero_iff (L m a : ℕ) : lookup L m a ≠ 0 ↔ + a < 2 ^ L ∧ population L a = m ∧ prefixCount L m a < K0 L m * K0 L m := by + simp [lookup] + +/-- Every inverse-table word is bounded by its original mask plus one. -/ +theorem lookup_le (L m a : ℕ) : lookup L m a ≤ a + 1 := by + unfold lookup + split + · exact Nat.add_le_add_right (prefixCount_le L m a) 1 + · exact Nat.zero_le _ + +/-- The inverse-table word at every selected block pair has the promised rank. -/ +theorem lookup_selectedMask {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : + lookup L m (selectedMask gh).val = (pairRank gh).val + 1 := by + have hm := (validMasks L m).orderEmbOfFin_mem (card_validMasks L m) (pairRank gh) + have hp := (mem_filter.mp hm).2 + change population L (selectedMask gh).val = m at hp + have hr : (pairRank gh).val < K0 L m * K0 L m := (finProdFinEquiv gh).isLt + simp only [lookup, (selectedMask gh).isLt, hp, prefixCount_selectedMask, + hr, and_self, ↓reduceIte] + +/-- A nonzero inverse-table word belongs to one of the selected block pairs. -/ +theorem exists_selectedMask_of_lookup {L m a : ℕ} (h : lookup L m a ≠ 0) : + ∃ gh : Fin (K0 L m) × Fin (K0 L m), (selectedMask gh).val = a := by + unfold lookup at h + split_ifs at h with hc + · obtain ⟨ha, hp, hr⟩ := hc + let aa : Fin (2 ^ L) := ⟨a, ha⟩ + have hmem : aa ∈ validMasks L m := mem_filter.mpr ⟨mem_univ _, hp⟩ + let j := (validMasks L m).orderIsoOfFin (card_validMasks L m) |>.symm ⟨aa, hmem⟩ + have hj : (validMasks L m).orderEmbOfFin (card_validMasks L m) j = aa := + congrArg Subtype.val ((validMasks L m).orderIsoOfFin (card_validMasks L m) + |>.apply_symm_apply ⟨aa, hmem⟩) + have hjrank : j.val = prefixCount L m a := by + rw [prefixCount_eq_bounded L m a ha.le] + have hbelow := card_below_orderEmb (validMasks L m) (card_validMasks L m) j + rw [hj] at hbelow + exact hbelow.symm + let gh := (@finProdFinEquiv (K0 L m) (K0 L m)).symm ⟨prefixCount L m a, hr⟩ + have hprval : (pairRank gh).val = prefixCount L m a := + congrArg Fin.val ((@finProdFinEquiv (K0 L m) (K0 L m)).apply_symm_apply + ⟨prefixCount L m a, hr⟩) + have hpr : pairRank gh = j := Fin.ext (hprval.trans hjrank.symm) + refine ⟨gh, ?_⟩ + change ((validMasks L m).orderEmbOfFin (card_validMasks L m) (pairRank gh)).val = a + rw [hpr, hj] + · exact (h rfl).elim + +/-- Inverse-table division and remainder recover exactly the two block indices. -/ +theorem lookup_blockPair {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : + (lookup L m (selectedMask gh).val - 1) / K0 L m = gh.1.val ∧ + (lookup L m (selectedMask gh).val - 1) % K0 L m = gh.2.val := by + rw [lookup_selectedMask, Nat.add_sub_cancel] + have hpos : 0 < K0 L m := by have := gh.2.isLt; omega + change (gh.2.val + K0 L m * gh.1.val) / K0 L m = gh.1.val ∧ + (gh.2.val + K0 L m * gh.1.val) % K0 L m = gh.2.val + rw [Nat.add_mul_div_left _ _ hpos, Nat.div_eq_of_lt gh.2.isLt, Nat.zero_add, + Nat.add_mul_mod_self_left, Nat.mod_eq_of_lt gh.2.isLt] + exact ⟨rfl, rfl⟩ + +/-- Both decoded block indices are in range whenever the inverse-table word is nonzero. -/ +theorem lookup_block_lt {L m a : ℕ} (h : lookup L m a ≠ 0) (rightSide : Bool) : + (if rightSide then (lookup L m a - 1) % K0 L m + else (lookup L m a - 1) / K0 L m) < K0 L m := by + obtain ⟨gh, hs⟩ := exists_selectedMask_of_lookup h + have hp := lookup_blockPair gh + rw [hs] at hp + cases rightSide + · simpa only [Bool.false_eq_true, ↓reduceIte, hp.1] using gh.1.isLt + · simpa only [↓reduceIte, hp.2] using gh.2.isLt + +/-- Selected masks always have exactly the required cardinality. -/ +theorem card_selectedMask {L m : ℕ} (gh : Fin (K0 L m) × Fin (K0 L m)) : + (maskSet (selectedMask gh)).card = m := by + have h := (validMasks L m).orderEmbOfFin_mem (card_validMasks L m) (pairRank gh) + exact (population_eq_card _).symm.trans ((mem_filter.mp h).2) + +/-- Distinct block pairs receive distinct level sets. -/ +theorem selectedMask_injective (L m : ℕ) : + Function.Injective (fun gh : Fin (K0 L m) × Fin (K0 L m) => maskSet (selectedMask gh)) := by + intro gh gh' h + have hmask := (maskEquiv L).injective h + have hrank := (validMasks L m).orderEmbOfFin (card_validMasks L m) |>.injective hmask + apply finProdFinEquiv.injective + apply Fin.ext + exact congrArg (fun i : Fin (K L m) => i.val) hrank + +/-- The thin-dimension pair digit uses row-major order, as do the encoder's inner variables. -/ +def innerDigits : Fin 4 ≃ Fin 2 × Fin 2 := (@finProdFinEquiv 2 2).symm + +/-- Explicit base-three outer strings and base-four strings of index pairs, with the first +`K0²` binary masks in increasing order as the shared block table. -/ +def layout {L m : ℕ} (hmL : m ≤ L) : Layout L m where + hmL := hmL + rowIdx := finFunctionFinEquiv.symm + colIdx := finFunctionFinEquiv.symm + innerIdx := finFunctionFinEquiv.symm.trans (Equiv.piCongrRight fun _ => innerDigits) + table gh := maskSet (selectedMask gh) + table_card := card_selectedMask + table_injective := selectedMask_injective L m + +end Cslib.Algorithms.MatrixMultiplication.Schonhage.MaskTable diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/SplitDigits.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/SplitDigits.lean new file mode 100644 index 0000000000..a610801d1a --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/SplitDigits.lean @@ -0,0 +1,103 @@ +/- +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.MaskTable +public import Complexity.Foundations.Data.Nat.OfDigits +public import Complexity.Foundations.Data.Finset.OrderedFilter + +/-! +# Splitting flat tensor digits into band indices + +Outer digits `0,1,2` contribute to a base-three offset. Inner digits `3,4,5,6` +contribute to a base-four offset and set the corresponding binary mask bit. All positions +are least significant first, in the encoder's existing numbering. + +## Main definitions + +* `BandDigits.parts`: the three indices computed by the native digit loop. +* `BandDigits.sevenDigits`: exactly `L` digits, retaining leading zeros. + +## Main results + +* `BandDigits.parts_bounds`: component bounds derived from the original depth. +* `BandDigits.parts_lists`: the compressed digit-list interpretation of all three indices. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage.BandDigits + +/-- Binary level mask and the two compressed offsets. -/ +structure Parts where + /-- Inner-level mask. -/ + mask : ℕ + /-- Base-three outer offset. -/ + outer : ℕ + /-- Base-four inner offset. -/ + inner : ℕ + deriving DecidableEq, Repr + +/-- Remove one base-seven digit, recursively interpreting the remaining higher levels. -/ +def parts : ℕ → ℕ → Parts + | 0, _ => ⟨0, 0, 0⟩ + | L + 1, a => + let tail := parts L (a / 7) + if a % 7 < 3 then ⟨2 * tail.mask, a % 7 + 3 * tail.outer, tail.inner⟩ + else ⟨1 + 2 * tail.mask, tail.outer, (a % 7 - 3) + 4 * tail.inner⟩ + +/-- Component fits depend only on the input depth, even for an unbounded argument. -/ +theorem parts_bounds (L a : ℕ) : + (parts L a).mask < 2 ^ L ∧ (parts L a).outer < 3 ^ L ∧ (parts L a).inner < 4 ^ L := by + induction L generalizing a with + | zero => simp [parts] + | succ L ih => + obtain ⟨hm, ho, hi⟩ := ih (a / 7) + have hd : a % 7 < 7 := Nat.mod_lt _ (by decide) + rw [parts] + split <;> dsimp only + · simp only [pow_succ] + constructor + · omega + constructor <;> omega + · simp only [pow_succ] + constructor + · omega + constructor <;> omega + +/-- A fixed number of digits, so that zero digits at the end are retained. -/ +def sevenDigits : ℕ → ℕ → List ℕ + | 0, _ => [] + | L + 1, a => a % 7 :: sevenDigits L (a / 7) + +/-- These are exactly the digits in Mathlib's finite-function numbering. -/ +theorem sevenDigits_eq_ofFn (L a : ℕ) : + sevenDigits L a = List.ofFn (fun i : Fin L => a / 7 ^ i.val % 7) := by + induction L generalizing a with + | zero => simp [sevenDigits] + | succ L ih => + rw [sevenDigits, ih, List.ofFn_succ] + simp only [Fin.val_zero, pow_zero, Nat.div_one, Fin.val_succ, pow_succ', + ← Nat.div_div_eq_div_mul] + +/-- The native indices are precisely the binary mask and the two compressed digit lists. -/ +theorem parts_lists (L a : ℕ) : + (parts L a).mask = Nat.ofDigits 2 ((sevenDigits L a).map fun d => if d < 3 then 0 else 1) ∧ + (parts L a).outer = Nat.ofDigits 3 ((sevenDigits L a).filter fun d => d < 3) ∧ + (parts L a).inner = Nat.ofDigits 4 + (((sevenDigits L a).filter fun d => ¬d < 3).map fun d => d - 3) := by + induction L generalizing a with + | zero => simp [parts, sevenDigits] + | succ L ih => + obtain ⟨hm, ho, hi⟩ := ih (a / 7) + by_cases hd : a % 7 < 3 + · simp [parts, sevenDigits, hd, Nat.ofDigits, hm, ho, hi] + · have hg : 3 ≤ a % 7 := by omega + simp [parts, sevenDigits, hd, hg, Nat.ofDigits, hm, ho, hi] + +end Cslib.Algorithms.MatrixMultiplication.Schonhage.BandDigits diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandAll.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandAll.lean new file mode 100644 index 0000000000..1cea6cc669 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandAll.lean @@ -0,0 +1,171 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandEncode + +/-! +# Every shared band, computed once + +## Main results + +* `BandPreparation.exec_encodeBands`: the complete native loop prepares, encodes and collects + all bands of one matrix side. Catalogue cells retain the existing flat leaf numbering. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +-- Static preservation checks reduce the finite persistent-register set. +set_option maxRecDepth 4096 + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Catalogue blocks are ordered by band index, without duplicating them per tile. -/ +theorem catalogue_address_lt {L m N β γ p : ℕ} (side : Bool) (hβγ : β < γ) + (hp : p < 10 ^ L) : + catalogueAddress L m N side β p < catalogueAddress L m N side γ 0 := by + have ht : 0 < 10 ^ L := Nat.pow_pos (by decide) + unfold catalogueAddress + nlinarith only [hβγ, hp, ht] + +/-- The full one-side catalogue interval ends at the block immediately after its last band. -/ +theorem catalogue_block_inside {L m N β : ℕ} (side : Bool) (hβ : β < numBands L m N) : + catalogueAddress L m N side 0 0 ≤ catalogueAddress L m N side β 0 ∧ + catalogueAddress L m N side β 0 + 10 ^ L ≤ + catalogueAddress L m N side (numBands L m N) 0 := by + have ht : 0 < 10 ^ L := Nat.pow_pos (by decide) + unfold catalogueAddress + constructor <;> nlinarith only [hβ, ht, Nat.zero_le β] + +/-- All ordinary band counters fit from the original explicit allocation budget. -/ +theorem bands_fits {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) : + Fits w (numBands L m N) := by + apply hfit.of_le + unfold wordBudget + have ht : 0 < 10 ^ L := Nat.pow_pos (by decide) + nlinarith only [ht, Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (numBands L m N)] + +/-- The actual native loop calls the complete encoder once for every shared band, charges +all loop control and counter increments, and preserves any outside caller memory. -/ +theorem exec_encodeBands {w L m N : ℕ} (hw : 0 < w) (side : Bool) (x : Data 96) + (hx : BandParameters L m N x) + (htable : ∀ t < 2 ^ L, x.mem (inverseBase L m N + t) = MaskTable.lookup L m t) + (matrix : ℕ → ℕ → ℤ) + (hmat : ∀ r < N, ∀ k < D m, + x.mem (matrixAddress L m N side r k) = encode (matrix r k)) + (A : ℤ) (hA : 0 ≤ A) (ha : ∀ r < N, ∀ k < D m, |matrix r k| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + ∃ y, ExecWithin (w + 2) (encodeBands side) x y + (numBands L m N * (bandTime L + 3) + 3) ∧ + BandParameters L m N y ∧ + (∀ β < numBands L m N, ∀ p < 10 ^ L, + y.mem (catalogueAddress L m N side β p) = + encode (catalogueValue L m N side matrix β p)) ∧ + (∀ b, ¬(1 ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + ¬(catalogueAddress L m N side 0 0 ≤ b ∧ + b < catalogueAddress L m N side (numBands L m N) 0) → + y.mem b = x.mem b) := by + let start := x.set band 0 + have hs : ExecWithin (B := .standard) (w + 2) (ram[ band ← #0 ]) x start 1 := + ExecWithin.const (fits_zero _) + have hp : BandParameters L m N start := parameters_reframe hx (by + intro i hi + have hne : i ≠ band := by + have hp : ∀ i ∈ parameterRegisters, i ≠ band := by decide + exact hp i hi + simp only [start, Data.set_regs, Function.update_of_ne hne]) + let I (y : Data 96) : Prop := BandParameters L m N y ∧ + y.regs band ≤ numBands L m N ∧ + (∀ β < y.regs band, ∀ p < 10 ^ L, + y.mem (catalogueAddress L m N side β p) = + encode (catalogueValue L m N side matrix β p)) ∧ + (∀ b, ¬(1 ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + ¬(catalogueAddress L m N side 0 0 ≤ b ∧ + b < catalogueAddress L m N side (numBands L m N) 0) → + y.mem b = x.mem b) + have hstart : I start := by + refine ⟨hp, ?_, ?_, ?_⟩ + · simp [start] + · simp [start] + · intro b _ _; rfl + have hbody (y : Data 96) (hy : I y) (hc : (Cond.lt band bands).Holds y) : + ∃ z, ExecWithin (w + 2) (bandBody side) y z (bandTime L + 1) ∧ I z ∧ + numBands L m N - z.regs band < numBands L m N - y.regs band := by + let β := y.regs band + have hβ : β < numBands L m N := by simpa only [Cond.holds_lt, hy.1.bands] using hc + have hlow := workspace_le_snapshot L m N + have hcat := catalogue_le_address L m N 0 side + have htable' (t : ℕ) (ht : t < 2 ^ L) : + y.mem (inverseBase L m N + t) = MaskTable.lookup L m t := by + rw [hy.2.2.2 _ (by unfold inverseBase; omega) (by + unfold catalogueBase tableBase at hcat + omega)] + exact htable t ht + have hmat' (r : ℕ) (hr : r < N) (k : ℕ) (hk : k < D m) : + y.mem (matrixAddress L m N side r k) = encode (matrix r k) := by + have hbound := matrixAddress_lt_inverse (L := L) side hr hk + rw [hy.2.2.2 _ (by unfold matrixAddress; omega) (by + unfold catalogueBase tableBase at hcat + have := Nat.zero_le (2 ^ L) + omega)] + exact hmat r hr k hk + obtain ⟨t, he, htp, htb, hvalues, hframe⟩ := exec_band hw side β y hy.1 hβ rfl + htable' matrix hmat' A hA ha hmag hfit + let z := t.set band (β + 1) + have hz : ExecWithin (B := .standard) (w + 2) (ram[ band ← band + one ]) t z 1 := by + ram_eval [z, htb, htp.fixed.one, (bands_fits hfit).of_le (Nat.succ_le_of_lt hβ), + band, one, Fin.mk.injEq] + have hzb : z.regs band = β + 1 := by simp [z] + refine ⟨z, he.seq hz, ?_, ?_⟩ + · refine ⟨parameters_reframe htp ?_, ?_, ?_, ?_⟩ + · intro i hi + have hp : ∀ i ∈ parameterRegisters, i ≠ band := by decide + simp only [z, Data.set_regs, Function.update_of_ne (hp i hi)] + · rw [hzb]; omega + · intro γ hγ p hp + rw [hzb] at hγ + change t.mem _ = _ + by_cases hγβ : γ = β + · subst γ; exact hvalues p hp + · have hγlt : γ < β := by omega + have hprior := catalogue_address_lt (m := m) (N := N) side hγlt hp + rw [hframe _ ?_ (by omega)] + · exact hy.2.2.1 γ hγlt p hp + · have haddr := catalogue_le_address L m N γ side + unfold catalogueAddress at haddr + have hsep : Encoding.scratchBase L + 10 ^ L ≤ catalogueBase L m N := by + unfold catalogueBase tableBase inverseBase + have := Nat.zero_le (2 ^ L) + omega + unfold catalogueAddress + omega + · intro b hb hc + change t.mem b = x.mem b + have hinside := catalogue_block_inside side hβ + rw [hframe b hb (by omega)] + exact hy.2.2.2 b hb hc + · rw [hzb]; change numBands L m N - (β + 1) < numBands L m N - β; omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => numBands L m N - y.regs band) hbody hstart + have hb : y.regs band = numBands L m N := by + simp only [Cond.holds_lt, hy.1.bands] at hc + have := hy.2.1 + omega + refine ⟨y, ?_, hy.1, ?_, hy.2.2.2⟩ + · apply (hs.seq he).mono + have hz : start.regs band = 0 := by simp [start] + simp only [hz, Nat.sub_zero, Nat.add_assoc, Nat.reduceAdd] + omega + · simpa only [hb] using hy.2.2.1 + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCatalogue.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCatalogue.lean new file mode 100644 index 0000000000..163d03c14a --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCatalogue.lean @@ -0,0 +1,70 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandFlat + +/-! +# Shared catalogue inputs for an existing mathematical tile + +## Main definitions + +* `BandPreparation.tileInputBases`: the two shared band pointers of a tile in `tiles`. + +## Main results + +* `BandPreparation.tile_input_values`: both pointed arrays contain the corresponding + unchanged public band encodings, in `Encoding.leafIndex` order. +* `BandPreparation.table_word_layout`: the computed forward table encodes the actual + deterministic `Layout` table, rather than caller-supplied indexing data. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- A tile `(rowBand,columnBand)` reads these two shared `10^L` arrays. The left pointer +depends only on its row band, and the right pointer only on its column band. -/ +def tileInputBases (L m N : ℕ) (T : ℕ × ℕ) : ℕ × ℕ := + (catalogueAddress L m N false T.1 0, catalogueAddress L m N true T.2 0) + +/-- The catalogue results returned by `runs` directly supply both inputs of every original +mathematical tile. Repeated row/column bands refer to the same stored arrays. -/ +theorem tile_input_values {L m N : ℕ} (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (y : Data 96) + (hleft : ∀ β < numBands L m N, ∀ τ, + y.mem (catalogueAddress L m N false β (Encoding.leafIndex τ).val) = + encode (encodingL (bandArrayL (MaskTable.layout hm) X β) τ)) + (hright : ∀ β < numBands L m N, ∀ τ, + y.mem (catalogueAddress L m N true β (Encoding.leafIndex τ).val) = + encode (encodingR (bandArrayR (MaskTable.layout hm) Y β) τ)) + (T : ℕ × ℕ) (hT : T ∈ tiles L m N) (τ : Leaf L) : + y.mem ((tileInputBases L m N T).1 + (Encoding.leafIndex τ).val) = + encode (encodingL (bandArrayL (MaskTable.layout hm) X T.1) τ) ∧ + y.mem ((tileInputBases L m N T).2 + (Encoding.leafIndex τ).val) = + encode (encodingR (bandArrayR (MaskTable.layout hm) Y T.2) τ) := by + have ht : T.1 < numBands L m N ∧ T.2 < numBands L m N := by + simpa only [tiles, Finset.mem_product, Finset.mem_range] using hT + simpa only [tileInputBases, catalogueAddress, Nat.add_zero] using + And.intro (hleft T.1 ht.1 τ) (hright T.2 ht.2 τ) + +/-- The retained native forward words encode exactly the computed mathematical layout. -/ +theorem table_word_layout {L m N : ℕ} (hm : m ≤ L) (y : Data 96) + (hforward : ∀ gh : Fin (K0 L m) × Fin (K0 L m), + y.mem (tableBase L m N + (MaskTable.pairRank gh).val) = + (MaskTable.selectedMask gh).val) + (gh : Fin (K0 L m) × Fin (K0 L m)) : + y.mem (tableBase L m N + (MaskTable.pairRank gh).val) = + (MaskTable.maskCode ((MaskTable.layout hm).table gh)).val := by + rw [hforward] + exact congrArg Fin.val (MaskTable.maskCode_maskSet (MaskTable.selectedMask gh)).symm + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCollect.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCollect.lean new file mode 100644 index 0000000000..ba061da3e0 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCollect.lean @@ -0,0 +1,119 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandCopy +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandParameters + +/-! +# Charged copies into the shared band catalogue + +## Main results + +* `BandPreparation.catalogue_block_le`: each band has one disjoint catalogue block. +* `BandPreparation.exec_collectBand`: the real native encoder-output copy, including + every address calculation and the unchanged copy kernel's instruction count. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +-- Static register-frame proofs reduce the fixed finite code and register bank. +set_option maxRecDepth 4096 + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Every complete band block lies in the catalogue allocation. -/ +theorem catalogue_block_le {L m N β : ℕ} (side : Bool) (hβ : β < numBands L m N) : + catalogueAddress L m N side β 0 + 10 ^ L ≤ memoryEnd L m N := by + have ht : 0 < 10 ^ L := Nat.pow_pos (by decide) + unfold catalogueAddress memoryEnd + cases side <;> simp only [Bool.false_eq_true, ↓reduceIte, Nat.zero_add, Nat.add_zero] + all_goals nlinarith only [hβ, ht] + +/-- No band output is placed in the protected snapshot, tables or native workspace. -/ +theorem catalogue_le_address (L m N β : ℕ) (side : Bool) : + catalogueBase L m N ≤ catalogueAddress L m N side β 0 := by + unfold catalogueAddress + omega + +/-- Address arithmetic and the full copy are charged. The copy's source and destination are +disjoint by the concrete allocation; arbitrary words outside its destination are retained. -/ +theorem exec_collectBand {w L m N : ℕ} (side : Bool) (β : ℕ) (x : Data 96) + (hx : BandParameters L m N x) (hβ : β < numBands L m N) (hb : x.regs band = β) + (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w (collectBand side) x y (7 * 10 ^ L + 10) ∧ + y.mem = ArrayCopy.copyPrefix x.mem (Encoding.outputBase L) + (catalogueAddress L m N side β 0) (10 ^ L) ∧ + BandParameters L m N y ∧ y.regs band = β := by + let b := (if side then numBands L m N else 0) + β + let dest := catalogueAddress L m N side β 0 + have hdest : Fits w (dest + 10 ^ L) := (memory_fits hfit).of_le + (catalogue_block_le side hβ) + have hd0 : Fits w dest := hdest.of_le (Nat.le_add_right _ _) + have hp : Fits w (b * 10 ^ L) := hd0.of_le (by + dsimp only [dest, catalogueAddress, b] + omega) + have hbfit : Fits w b := hp.of_le (Nat.le_mul_of_pos_right _ + (Nat.pow_pos (by decide))) + have hs : Fits w (7 ^ L + 2) := (encoder_fits hfit).of_le (by + have h7 := Nat.pow_le_pow_left (by decide : 7 ≤ 10) L + unfold Encoding.addressBound + nlinarith only [h7, Nat.zero_le L, Nat.zero_le (10 ^ L)]) + have hsrc : Fits w (Encoding.outputBase L + 10 ^ L) := + ((memory_fits hfit).of_le (snapshot_le_memoryEnd L m N)).of_le (by + have h := workspace_le_snapshot L m N + unfold Encoding.scratchBase at h + have := Nat.zero_le (10 ^ L) + omega) + let start := ((((((x.set copySrc (7 ^ L + 2)).set temporary b).set temporary + (b * 10 ^ L)).set temporary dest).set copyDst dest).set copyLen (10 ^ L)) + dsimp only [dest, catalogueAddress] at hd0 + have hhead : ExecWithin w (collectHead side) x start 6 := by + cases side <;> + simp only [b, Bool.false_eq_true, ↓reduceIte, Nat.zero_add, Nat.add_zero] at hbfit hp hd0 <;> + ram_eval [collectHead, start, dest, b, catalogueAddress, hx.cells, hx.leaves, + hx.catalogue, hx.bands, hx.fixed.two, hb, hs, hd0, hp, hbfit, + copySrc, copyDst, copyLen, temporary, cells, two, leaves, catalogue, bands, + band, Fin.mk.injEq, Nat.zero_add, Nat.add_zero] + obtain ⟨y, he, hmem, _⟩ := exec_copyKernel start (small_fits hfit (by decide)) + (by simpa [start, copySrc, copyLen, temporary, Encoding.outputBase, Nat.add_comm] using hsrc) + (by simpa [start, copyDst, copyLen, temporary] using hdest) + (Or.inl (by + have hlow := workspace_le_snapshot L m N + have hcat : snapshotBase L m N ≤ catalogueBase L m N := by + unfold catalogueBase tableBase inverseBase + have := Nat.zero_le (2 ^ L) + omega + have hsep : Encoding.outputBase L + 10 ^ L ≤ dest := by + have hn := Nat.zero_le (10 ^ L) + unfold Encoding.scratchBase at hlow + exact (by omega : Encoding.outputBase L + 10 ^ L ≤ snapshotBase L m N).trans + (hcat.trans (catalogue_le_address L m N β side)) + simpa [start, copySrc, copyDst, copyLen, temporary, Encoding.outputBase, Nat.add_comm] + using hsep)) + have hwhole : ExecWithin w (collectBand side) x y (7 * 10 ^ L + 10) := by + apply (hhead.seq he).mono + have hn : start.regs copyLen = 10 ^ L := by simp [start] + rw [hn] + omega + refine ⟨y, hwhole, ?_, parameters_reframe hx ?_, ?_⟩ + · simpa [start, copySrc, copyDst, copyLen, temporary, Encoding.outputBase, Nat.add_comm] + using hmem + · intro i hi + apply hwhole.regs_eq + have hp : ∀ side i, i ∈ parameterRegisters → (collectBand side).Preserves i := by + intro side + cases side <;> simp +decide only [] + exact hp side i hi + · exact (hwhole.regs_eq (by + cases side <;> simp +decide only [])).trans hb + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCopy.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCopy.lean new file mode 100644 index 0000000000..151981f38d --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCopy.lean @@ -0,0 +1,134 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandMemory +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandHeader + +/-! +# Native copies used by the band-preparation caller + +## Main results + +* `BandPreparation.exec_copyKernel`: reuse the actual array-copy program in its fixed low + register bank, preserving every persistent register and permitting both address orders. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM + +/-- The seven low copy registers are embedded injectively. -/ +theorem copyRegisters_injective : Function.Injective copyRegisters := by + intro i j h + exact Fin.ext (congrArg (fun i : Fin 96 => i.val) h) + +/-- The unchanged encoder occupies precisely its original low register bank. -/ +theorem encoderRegisters_injective : Function.Injective encoderRegisters := by + intro i j h + exact Fin.ext (congrArg (fun i : Fin 96 => i.val) h) + +/-- Copying uses only the original seven-register routine and its actual native cost. +The persistent high register bank is untouched; arbitrary outside memory is preserved +by the explicit `copyPrefix` result. -/ +theorem exec_copyKernel {w : ℕ} (x : Data 96) (h1 : Fits w 1) + (hs : Fits w (x.regs copySrc + x.regs copyLen)) + (hd : Fits w (x.regs copyDst + x.regs copyLen)) + (hsep : x.regs copySrc + x.regs copyLen ≤ x.regs copyDst ∨ + x.regs copyDst + x.regs copyLen ≤ x.regs copySrc ∨ x.regs copySrc = x.regs copyDst) : + ∃ y, ExecWithin w (ArrayCopy.program.rename copyRegisters) x y (7 * x.regs copyLen + 4) ∧ + y.mem = ArrayCopy.copyPrefix x.mem (x.regs copySrc) (x.regs copyDst) (x.regs copyLen) ∧ + (∀ i, 7 ≤ i.val → y.regs i = x.regs i) := by + obtain ⟨z, hz, hm, _, _, _⟩ := ArrayCopy.exec_of_separated_or_eq (x.comap copyRegisters) + h1 hs hd hsep + have he := hz.rename copyRegisters_injective + refine ⟨z.extend copyRegisters x, he, hm, ?_⟩ + intro i hi + apply he.regs_eq + have hf : ∀ i : Fin 96, 7 ≤ i.val → (ArrayCopy.program.rename copyRegisters).Preserves i := by + simp only [ArrayCopy.program, ArrayCopy.body, Code.rename, Code.Preserves, + Op.rename, Op.bin, Op.dest, copyRegisters] + decide + exact hf i hi + +private theorem lowSetup_regs (x : Data 96) (s d n : ℕ) (i : Fin 96) (hi : 7 ≤ i.val) : + (((x.set copySrc s).set copyDst d).set copyLen n).regs i = x.regs i := by + have hs : i ≠ copySrc := by intro h; subst i; norm_num [copySrc] at hi + have hd : i ≠ copyDst := by intro h; subst i; norm_num [copyDst] at hi + have hn : i ≠ copyLen := by intro h; subst i; norm_num [copyLen] at hi + simp only [Data.set_regs, Function.update_of_ne hs, Function.update_of_ne hd, + Function.update_of_ne hn] + +/-- Protect the entire original initialized prefix, including header zero, with the actual +native copy. Source words and all memory outside the snapshot block remain untouched. -/ +theorem exec_saveInput {w L m N : ℕ} (x : Data 96) (hx : HeaderPost L m N x) + (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w saveInput x y (7 * inputEnd m N + 7) ∧ + y.mem = ArrayCopy.copyPrefix x.mem 0 (snapshotBase L m N) (inputEnd m N) ∧ + (∀ i, 7 ≤ i.val → y.regs i = x.regs i) := by + let start := ((x.set copySrc 0).set copyDst (snapshotBase L m N)).set copyLen (inputEnd m N) + have hb : Fits w (inputEnd m N) := (memory_fits hfit).of_le (input_le_memoryEnd L m N) + have hi : Fits w (inverseBase L m N) := + (memory_fits hfit).of_le (inverse_le_memoryEnd L m N) + obtain ⟨y, hy, hm, hf⟩ := exec_copyKernel start (small_fits hfit (by decide)) + (by simpa [start, copySrc, copyDst, copyLen] using hb) + (by simpa [start, copySrc, copyDst, copyLen, inverseBase] using hi) + (Or.inl (by simpa [start, copySrc, copyDst, copyLen] using input_le_snapshot L m N)) + have h0 : ExecWithin (B := .standard) w (.op (.const copySrc 0)) + x (x.set copySrc 0) 1 := ExecWithin.const (fits_zero w) + have h1 : ExecWithin (B := .standard) w (.op (.copy copyDst snapshot)) + (x.set copySrc 0) ((x.set copySrc 0).set copyDst (snapshotBase L m N)) 1 := by + ram_eval [hx.snapshot, copySrc, copyDst, snapshot, Fin.mk.injEq] + have h2 : ExecWithin (B := .standard) w (.op (.copy copyLen prefixSize)) + ((x.set copySrc 0).set copyDst (snapshotBase L m N)) start 1 := by + ram_eval [start, hx.prefixSize, copySrc, copyDst, copyLen, prefixSize, Fin.mk.injEq] + refine ⟨y, ?_, ?_, ?_⟩ + · apply (h0.seq (h1.seq (h2.seq hy))).mono + have hn : start.regs copyLen = inputEnd m N := by simp [start] + rw [hn] + omega + · simpa [start, copySrc, copyDst, copyLen] using hm + · intro i hi + exact (hf i hi).trans (lowSetup_regs x _ _ _ i hi) + +/-- Restore the complete original input using the reverse disjoint-copy case. This copies +the actual protected words rather than reconstructing or assuming the input is unchanged. -/ +theorem exec_restoreInput {w L m N : ℕ} (x : Data 96) (hx : HeaderPost L m N x) + (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w restoreInput x y (7 * inputEnd m N + 7) ∧ + y.mem = ArrayCopy.copyPrefix x.mem (snapshotBase L m N) 0 (inputEnd m N) ∧ + (∀ i, 7 ≤ i.val → y.regs i = x.regs i) := by + let start := ((x.set copySrc (snapshotBase L m N)).set copyDst 0).set copyLen (inputEnd m N) + have hb : Fits w (inputEnd m N) := (memory_fits hfit).of_le (input_le_memoryEnd L m N) + have hi : Fits w (inverseBase L m N) := + (memory_fits hfit).of_le (inverse_le_memoryEnd L m N) + obtain ⟨y, hy, hm, hf⟩ := exec_copyKernel start (small_fits hfit (by decide)) + (by simpa [start, copySrc, copyDst, copyLen, inverseBase] using hi) + (by simpa [start, copySrc, copyDst, copyLen] using hb) + (Or.inr (Or.inl (by + simpa [start, copySrc, copyDst, copyLen] using input_le_snapshot L m N))) + have h0 : ExecWithin (B := .standard) w (.op (.copy copySrc snapshot)) + x (x.set copySrc (snapshotBase L m N)) 1 := by ram_eval [hx.snapshot] + have h1 : ExecWithin (B := .standard) w (.op (.const copyDst 0)) + (x.set copySrc (snapshotBase L m N)) + ((x.set copySrc (snapshotBase L m N)).set copyDst 0) 1 := ExecWithin.const (fits_zero w) + have h2 : ExecWithin (B := .standard) w (.op (.copy copyLen prefixSize)) + ((x.set copySrc (snapshotBase L m N)).set copyDst 0) start 1 := by + ram_eval [start, hx.prefixSize, copySrc, copyDst, copyLen, prefixSize, Fin.mk.injEq] + refine ⟨y, ?_, ?_, ?_⟩ + · apply (h0.seq (h1.seq (h2.seq hy))).mono + have hn : start.regs copyLen = inputEnd m N := by simp [start] + rw [hn] + omega + · simpa [start, copySrc, copyDst, copyLen] using hm + · intro i hi + exact (hf i hi).trans (lowSetup_regs x _ _ _ i hi) + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCost.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCost.lean new file mode 100644 index 0000000000..f8c74f3fa7 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCost.lean @@ -0,0 +1,70 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.BandBounds + +/-! +# Numerical cost bounds for shared band preparation + +## Main definitions + +* `BandPreparation.timeBound`: a uniform budget for all setup, tables, input copying, + band preparation, shared encoder calls and catalogue copies. + +## Main results + +* `BandPreparation.thin_timeBound`: its explicit finite thin-regime bound. This numerical + lemma is separate from the theorem connecting the budget to the actual native execution. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Generic numerical instruction budget; a final halt is charged separately. -/ +def timeBound (L m N : ℕ) : ℕ := + 4096 * (L + 1) * (1 + N * D m + 2 ^ L + 2 * numBands L m N * 10 ^ L) + +/-- At `L=19m`, all original input and table overhead is absorbed by the saving for shared +encodings. The displayed constants and polynomial exponent are fixed and positive. -/ +theorem thin_timeBound (m N : ℕ) (hN : D m ^ 18 ≤ N) : + (timeBound (19 * m) m N : ℝ) ≤ + 40960 * ((19 * m + 1 : ℕ) : ℝ) * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2 := by + have hNpos : 0 < N := (Nat.pow_pos (show 0 < D m from Nat.pow_pos (by decide))).trans_le hN + have hs := preparation_shape_le (show m ≤ 19 * m by omega) hNpos + have hsR : ((1 + N * D m + 2 ^ (19 * m) + + 2 * numBands (19 * m) m N * 10 ^ (19 * m) : ℕ) : ℝ) ≤ + 5 * ((numBands (19 * m) m N * 10 ^ (19 * m) : ℕ) : ℝ) := by exact_mod_cast hs + have he := encodings_total_le m N hN + have hshape : ((1 + N * D m + 2 ^ (19 * m) + + 2 * numBands (19 * m) m N * 10 ^ (19 * m) : ℕ) : ℝ) ≤ + 10 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2 := by + push_cast at hsR he ⊢ + nlinarith only [hsR, he] + calc (timeBound (19 * m) m N : ℝ) + = 4096 * ((19 * m + 1 : ℕ) : ℝ) * ((1 + N * D m + 2 ^ (19 * m) + + 2 * numBands (19 * m) m N * 10 ^ (19 * m) : ℕ) : ℝ) := by + simp only [timeBound, Nat.cast_mul, Nat.cast_ofNat] + _ ≤ 4096 * ((19 * m + 1 : ℕ) : ℝ) * + (10 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2) := + mul_le_mul_of_nonneg_left hshape (by positivity) + _ = _ := by ring + +/-- Equivalent division form of the thin-regime bound requested by the matrix algorithm. -/ +theorem thin_timeBound_div (m N : ℕ) (hN : D m ^ 18 ≤ N) : + (timeBound (19 * m) m N : ℝ) ≤ + 40960 * ((19 * m + 1 : ℕ) : ℝ) * (N : ℝ) ^ 2 / (D m : ℝ) ^ (1 / 18 : ℝ) := by + have h := thin_timeBound m N hN + rw [Real.rpow_neg (by positivity)] at h + simpa only [div_eq_mul_inv, mul_assoc, mul_left_comm, mul_comm] using h + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCount.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCount.lean new file mode 100644 index 0000000000..6fb297b6e9 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandCount.lean @@ -0,0 +1,112 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandPopulation + +/-! +# Native counting of all valid subset masks + +## Main results + +* `BandPreparation.exec_countMasks`: the complete binary-mask scan computes the binomial + coefficient, with its inner population loops and all outer controls charged. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open MaskTable + +/-- Scanning all masks computes `K L m`; neither that count nor a valid table is supplied +as an input. Every mask-counter and accepted-counter fit follows from the scan range. -/ +theorem exec_countMasks {w : ℕ} (L m : ℕ) (x : Data 96) + (hL : x.regs levels = L) (hm : x.regs inner = m) (hM : x.regs masks = 2 ^ L) + (h1 : x.regs one = 1) (h2 : x.regs two = 2) (hfit : Fits w (2 ^ L + L + 4)) : + ∃ y, ExecWithin w countMasks x y ((6 * L + 11) * 2 ^ L + 4) ∧ + y.regs count = K L m ∧ y.regs mask = 2 ^ L ∧ y.mem = x.mem ∧ + (∀ i, i ≠ mask → i ≠ count → i ≠ residue → i ≠ population → i ≠ counter → + i ≠ digit → y.regs i = x.regs i) := by + have h0 : Fits w 0 := hfit.of_le (by omega) + let start := (x.set mask 0).set count 0 + have hs : ExecWithin (B := .standard) w (ram[ mask ← #0; count ← #0 ]) x start 2 := by + ram_eval [start, h0, mask, count, Fin.mk.injEq] + let I (y : Data 96) : Prop := + y.mem = x.mem ∧ y.regs levels = L ∧ y.regs inner = m ∧ y.regs masks = 2 ^ L ∧ + y.regs one = 1 ∧ y.regs two = 2 ∧ y.regs mask ≤ 2 ^ L ∧ + y.regs count = prefixCount L m (y.regs mask) + have hi : I start := by + simp [I, start, hL, hm, hM, h1, h2, levels, inner, masks, one, two, mask, count, + prefixCount_zero] + have hb (y : Data 96) (hy : I y) (hc : (Cond.lt mask masks).Holds y) : + ∃ z, ExecWithin w countBody y z (6 * L + 9) ∧ I z ∧ + 2 ^ L - z.regs mask < 2 ^ L - y.regs mask := by + obtain ⟨hmem, hL, hm, hM, h1, h2, ht, hcount⟩ := hy + have htlt : y.regs mask < 2 ^ L := by simpa only [Cond.holds_lt, hM] using hc + obtain ⟨u, hu, hpop, humem, huf⟩ := exec_population L (y.regs mask) y hL rfl h1 h2 + (hfit.of_le (by omega)) + have hul : u.regs levels = L := + (huf levels (by decide) (by decide) (by decide) (by decide)).trans hL + have hui : u.regs inner = m := + (huf inner (by decide) (by decide) (by decide) (by decide)).trans hm + have huM : u.regs masks = 2 ^ L := + (huf masks (by decide) (by decide) (by decide) (by decide)).trans hM + have hu1 : u.regs one = 1 := + (huf one (by decide) (by decide) (by decide) (by decide)).trans h1 + have hu2 : u.regs two = 2 := + (huf two (by decide) (by decide) (by decide) (by decide)).trans h2 + have hut : u.regs mask = y.regs mask := + huf mask (by decide) (by decide) (by decide) (by decide) + have huc : u.regs count = prefixCount L m (y.regs mask) := + (huf count (by decide) (by decide) (by decide) (by decide)).trans hcount + have hfnext : Fits w (y.regs mask + 1) := hfit.of_le (by omega) + have hfc : Fits w (prefixCount L m (y.regs mask) + 1) := + hfit.of_le (by have := prefixCount_le L m (y.regs mask); omega) + let z := (u.set count (prefixCount L m (y.regs mask + 1))).set mask (y.regs mask + 1) + have hz : ExecWithin w (ram[ + if population = inner then [ count ← count + one ] else [ count ← count ]; + mask ← mask + one ]) u z 4 := by + by_cases h : MaskTable.population L (y.regs mask) = m <;> + ram_eval [z, prefixCount_succ, h, hpop, hui, hu1, hut, huc, hfnext, hfc, + population, inner, count, one, mask, Fin.mk.injEq] + refine ⟨z, hu.seq hz, ?_, ?_⟩ + · simpa +decide [I, z, levels, inner, masks, one, two, mask, count, + Function.update_apply] using + (show u.mem = x.mem ∧ u.regs levels = L ∧ u.regs inner = m ∧ + u.regs masks = 2 ^ L ∧ u.regs one = 1 ∧ u.regs two = 2 ∧ + y.regs mask + 1 ≤ 2 ^ L from + ⟨humem.trans hmem, hul, hui, huM, hu1, hu2, by omega⟩) + · have hzt : z.regs mask = y.regs mask + 1 := by simp [z] + rw [hzt] + omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => 2 ^ L - y.regs mask) hb hi + have ht : y.regs mask = 2 ^ L := by + simp only [Cond.holds_lt, hy.2.2.2.1] at hc + have := hy.2.2.2.2.2.2.1 + omega + have hwhole : ExecWithin w countMasks x y ((6 * L + 11) * 2 ^ L + 4) := by + apply (hs.seq he).mono + have hs0 : start.regs mask = 0 := by simp [start, mask, count] + rw [hs0, Nat.sub_zero] + nlinarith + refine ⟨y, hwhole, ?_, ht, hy.1, ?_⟩ + · have h := hy.2.2.2.2.2.2.2 + simpa only [ht, prefixCount_full] using h + · intro i hm hc hr hp hj hd + apply hwhole.regs_eq + have hf : ∀ i : Fin 96, i ≠ mask → i ≠ count → i ≠ residue → i ≠ population → + i ≠ counter → i ≠ digit → countMasks.Preserves i := by + simp only [countMasks, countBody, populationCode, Code.Preserves, Op.dest, Op.bin] + decide + exact hf i hm hc hr hp hj hd + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandEncode.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandEncode.lean new file mode 100644 index 0000000000..a3831dde56 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandEncode.lean @@ -0,0 +1,137 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandPrepare +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCollect + +/-! +# Actual native encoder calls on the prepared bands + +## Main definitions + +* `BandPreparation.catalogueValue`: flat shared-transform values, using the existing + coordinate transform and the computed, padded band array. + +## Main results + +* `BandPreparation.exec_band`: one complete native preparation, unchanged encoder call and + catalogue copy, with charged instructions and an explicit memory frame. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +-- The persistent-register checks reduce the fixed 96-register parameter bank. +set_option maxRecDepth 4096 + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The shared catalogue uses exactly the existing flat coordinate-transform values. +`Encoding.value_final_left/right` and `bandValue_left/right` identify these with the public +`encodingL/R` of the original mathematical `bandArrayL/R`. -/ +def catalogueValue (L m N : ℕ) (side : Bool) (matrix : ℕ → ℕ → ℤ) (β p : ℕ) : ℤ := + CoordinateTransform.value (CoordinateTransform.table (Encoding.coefficients side)) + L 0 (bandValue L m N β side matrix) p 0 + +/-- The guarded band only reads real bounded matrix entries; all other entries are zero. -/ +theorem abs_bandValue_le {L m N β : ℕ} (side : Bool) (matrix : ℕ → ℕ → ℤ) + (A : ℤ) (hA : 0 ≤ A) (hmat : ∀ r < N, ∀ k < D m, |matrix r k| ≤ A) + (a : Fin (7 ^ L)) : |bandValue L m N β side matrix a.val| ≤ A := by + cases side <;> unfold bandValue <;> simp only [Bool.false_eq_true, ↓reduceIte] + all_goals + split_ifs with hq hr + · simpa only [abs_zero] using hA + · exact hmat _ hr _ (selected_parts_bounds a hq).2 + · simpa only [abs_zero] using hA + +/-- A real call of the unchanged encoder in its low 32-register bank. Parameters in the +other registers, arbitrary input headers and arbitrary outside memory are preserved. -/ +theorem exec_encoderKernel {w L m N : ℕ} (hw : 0 < w) (side : Bool) + (arr : ℕ → ℤ) (x : Data 96) (hheader : x.mem 1 = L) + (hinput : ∀ i < 7 ^ L, x.mem (2 + i) = encode (arr i)) + (A : ℤ) (hA : 0 ≤ A) (ha : ∀ i < 7 ^ L, |arr i| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + ∃ y, ExecWithin (w + 2) ((Encoding.program side).rename encoderRegisters) x y + (Encoding.timeBound L) ∧ + (∀ p < 10 ^ L, y.mem (Encoding.outputBase L + p) = encode + (CoordinateTransform.value (CoordinateTransform.table (Encoding.coefficients side)) + L 0 arr p 0)) ∧ + (∀ b, ¬(Encoding.outputBase L ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + y.mem b = x.mem b) ∧ (∀ i, 32 ≤ i.val → y.regs i = x.regs i) := by + obtain ⟨z, hz, hv, hf, _, _⟩ := Encoding.exec hw side arr (x.comap encoderRegisters) + hheader hinput A hA ha hmag (encoder_fits hfit) + refine ⟨z.extend encoderRegisters x, hz.rename encoderRegisters_injective, hv, hf, ?_⟩ + intro i hi + change Function.extend encoderRegisters z.regs x.regs i = x.regs i + apply Function.extend_apply' + rintro ⟨j, hj⟩ + have hv := congrArg (fun k : Fin 96 => k.val) hj + change j.val = i.val at hv + omega + +/-- Native cost of preparing, encoding and collecting one shared band. -/ +def bandTime (L : ℕ) : ℕ := + (13 * L + 39) * 7 ^ L + 5 + Encoding.timeBound L + (7 * 10 ^ L + 10) + +/-- Every band invokes the actual encoder code and copies its complete output. Only the low +workspace and this single shared catalogue block can change. Original matrix/table data and +arbitrary caller query regions outside those blocks remain untouched. -/ +theorem exec_band {w L m N : ℕ} (hw : 0 < w) (side : Bool) (β : ℕ) (x : Data 96) + (hx : BandParameters L m N x) (hβ : β < numBands L m N) (hb : x.regs band = β) + (htable : ∀ t < 2 ^ L, x.mem (inverseBase L m N + t) = MaskTable.lookup L m t) + (matrix : ℕ → ℕ → ℤ) + (hmat : ∀ r < N, ∀ k < D m, + x.mem (matrixAddress L m N side r k) = encode (matrix r k)) + (A : ℤ) (hA : 0 ≤ A) (ha : ∀ r < N, ∀ k < D m, |matrix r k| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + ∃ y, ExecWithin (w + 2) + (.seq (prepareBand side) (.seq ((Encoding.program side).rename encoderRegisters) + (collectBand side))) x y (bandTime L) ∧ + BandParameters L m N y ∧ y.regs band = β ∧ + (∀ p < 10 ^ L, y.mem (catalogueAddress L m N side β p) = + encode (catalogueValue L m N side matrix β p)) ∧ + (∀ b, ¬(1 ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + ¬(catalogueAddress L m N side β 0 ≤ b ∧ + b < catalogueAddress L m N side β 0 + 10 ^ L) → y.mem b = x.mem b) := by + obtain ⟨p, hep, hpp, hpb, hheader, hinput, hpf⟩ := + exec_prepareBand side β x hx hβ hb htable matrix hmat hfit + obtain ⟨e, hee, hev, hef, her⟩ := exec_encoderKernel hw side + (bandValue L m N β side matrix) p hheader hinput A hA + (fun i hi => abs_bandValue_le side matrix A hA ha ⟨i, hi⟩) hmag hfit + have hep' : BandParameters L m N e := parameters_reframe hpp (by + intro i hi + apply her + have hp : ∀ i ∈ parameterRegisters, 32 ≤ i.val := by decide + exact hp i hi) + have heb : e.regs band = β := (her band (by decide)).trans hpb + obtain ⟨y, hec, hmem, hyp, hyb⟩ := exec_collectBand side β e hep' hβ heb hfit + refine ⟨y, ?_, hyp, hyb, ?_, ?_⟩ + · simpa only [bandTime, Nat.add_assoc] using hep.seq (hee.seq hec) + · intro p hp + have heq : catalogueAddress L m N side β p = + catalogueAddress L m N side β 0 + p := by simp [catalogueAddress] + rw [heq, hmem, ArrayCopy.copyPrefix_read hp] + exact hev p hp + · intro b hb hc + rw [hmem, ArrayCopy.copyPrefix_outside hc] + calc + e.mem b = p.mem b := hef b (by + unfold Encoding.outputBase + have := Nat.zero_le (7 ^ L) + omega) + _ = x.mem b := hpf b (by + unfold Encoding.scratchBase Encoding.outputBase at hb + have := Nat.zero_le (10 ^ L) + omega) + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandExec.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandExec.lean new file mode 100644 index 0000000000..c299842d75 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandExec.lean @@ -0,0 +1,245 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandAll +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandTables +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCost + +/-! +# Complete native shared-band execution + +## Main results + +* `BandPreparation.exec`: compute both tables and every shared encoding, from real scalar + headers and matrix words, restoring the full original input prefix and every outside cell. + All setup, copies, loops and encoder calls contribute to the instruction bound. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +-- Static frame checks reduce the fixed parameter bank and finite native table code. +set_option maxRecDepth 4096 + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The two tables and the protected input all precede the output catalogue. -/ +theorem snapshot_le_catalogue (L m N : ℕ) : snapshotBase L m N ≤ catalogueBase L m N := by + unfold catalogueBase tableBase inverseBase + have := Nat.zero_le (2 ^ L) + omega + +/-- The one-side catalogue end is within the complete allocation. -/ +theorem catalogue_side_end_le (L m N : ℕ) (side : Bool) : + catalogueAddress L m N side (numBands L m N) 0 ≤ memoryEnd L m N := by + unfold catalogueAddress memoryEnd + cases side <;> simp only [Bool.false_eq_true, ↓reduceIte, Nat.zero_add, Nat.add_zero] + all_goals nlinarith only [Nat.zero_le (numBands L m N * 10 ^ L)] + +private theorem band_frame_protected {L m N : ℕ} {side : Bool} {μ ν : ℕ → ℕ} + (h : ∀ b, ¬(1 ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + ¬(catalogueAddress L m N side 0 0 ≤ b ∧ + b < catalogueAddress L m N side (numBands L m N) 0) → ν b = μ b) + (a : ℕ) (ha : snapshotBase L m N ≤ a ∧ a < catalogueBase L m N) : ν a = μ a := by + apply h + · have := workspace_le_snapshot L m N + omega + · have := catalogue_le_address L m N 0 side + omega + +private theorem band_frame_above {L m N : ℕ} {side : Bool} {μ ν : ℕ → ℕ} + (h : ∀ b, ¬(1 ≤ b ∧ b < Encoding.scratchBase L + 10 ^ L) → + ¬(catalogueAddress L m N side 0 0 ≤ b ∧ + b < catalogueAddress L m N side (numBands L m N) 0) → ν b = μ b) + (a : ℕ) (ha : memoryEnd L m N ≤ a) : ν a = μ a := by + apply h + · have := workspace_le_snapshot L m N + have := (snapshot_le_catalogue L m N).trans (catalogue_le_memoryEnd L m N) + omega + · have := catalogue_side_end_le L m N side + omega + +/-- A coarse uniform bound for one actual shared-band call and its outer loop control. -/ +theorem bandTime_le (L : ℕ) : bandTime L + 3 ≤ 512 * (L + 1) * 10 ^ L := by + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + have ht : 1 ≤ 10 ^ L := Nat.one_le_pow L 10 (by decide) + have hs := Nat.mul_le_mul_left (13 * L + 39) h7 + unfold bandTime Encoding.timeBound + nlinarith only [hs, ht, Nat.zero_le L, Nat.zero_le (L * 10 ^ L)] + +/-- The sum of the previously proved native component charges fits the generic bound. -/ +theorem execution_time_le (L m N : ℕ) : + (16 * L + 49) + ((7 * inputEnd m N + 7) + + (((6 * L + 16) * 2 ^ L + 14) + (((6 * L + 17) * 2 ^ L + 4) + + ((numBands L m N * (bandTime L + 3) + 3) + + ((numBands L m N * (bandTime L + 3) + 3) + (7 * inputEnd m N + 7)))))) ≤ + timeBound L m N := by + have hb := Nat.mul_le_mul_left (numBands L m N) (bandTime_le L) + have h2 : 1 ≤ 2 ^ L := Nat.one_le_pow L 2 (by decide) + unfold timeBound inputEnd + nlinarith only [hb, h2, Nat.zero_le L, Nat.zero_le (N * D m), + Nat.zero_le (L * (N * D m)), Nat.zero_le (L * 2 ^ L), + Nat.zero_le (numBands L m N * 10 ^ L), + Nat.zero_le (L * (numBands L m N * 10 ^ L))] + +/-- The full fixed native program. The caller supplies only ordinary `L,m,N` headers and +signed matrix words; both tables and every padded/encoded band are computed. Header zero +and arbitrary outside memory need not be zero or describe a standalone input list. -/ +theorem exec {w L m N : ℕ} (hw : 0 < w) (hm : m ≤ L) (x : Data 96) + (hL : x.mem 1 = L) (hm' : x.mem 2 = m) (hN : x.mem 3 = N) + (matrix : Bool → ℕ → ℕ → ℤ) + (hmat : ∀ side r, r < N → ∀ k, k < D m → + x.mem (matrixOffset m N side r k + 4) = encode (matrix side r k)) + (A : ℤ) (hA : 0 ≤ A) (ha : ∀ side r, r < N → ∀ k, k < D m → |matrix side r k| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + ∃ y, ExecWithin (w + 2) program x y (timeBound L m N) ∧ + BandParameters L m N y ∧ + (∀ side β, β < numBands L m N → ∀ p, p < 10 ^ L → + y.mem (catalogueAddress L m N side β p) = + encode (catalogueValue L m N side (matrix side) β p)) ∧ + (∀ a < 2 ^ L, y.mem (inverseBase L m N + a) = MaskTable.lookup L m a) ∧ + (∀ gh : Fin (K0 L m) × Fin (K0 L m), + y.mem (tableBase L m N + (MaskTable.pairRank gh).val) = + (MaskTable.selectedMask gh).val) ∧ + (∀ b, ¬(inputEnd m N ≤ b ∧ b < memoryEnd L m N) → y.mem b = x.mem b) := by + obtain ⟨h, heh, hhp, hhmem⟩ := exec_header hm x hL hm' hN hfit + obtain ⟨s, hes, hsmem, hsr⟩ := exec_saveInput h hhp hfit + have hsp : HeaderPost L m N s := headerPost_reframe hhp (by + intro i hi + have hp : ∀ i ∈ headerRegisters, 7 ≤ i.val := by decide + exact hsr i (hp i hi)) + obtain ⟨c, hec, hcp, hcmem⟩ := exec_countCode s hsp hfit + have htableParams : TableParameters L m (inverseBase L m N) (tableBase L m N) c := + ⟨hcp.fixed.levels, hcp.fixed.inner, hcp.masks, hcp.inverse, hcp.table, hcp.square, + hcp.fixed.zero, hcp.fixed.one, hcp.fixed.two⟩ + obtain ⟨t, het, _, hlookup, hforward, htframe⟩ := exec_tables L m + (inverseBase L m N) (tableBase L m N) c htableParams (Nat.le_refl _) + (tables_fit hfit) + have htp : BandParameters L m N t := parameters_reframe hcp (by + intro i hi + apply het.regs_eq + have hp : ∀ i ∈ parameterRegisters, tablesCode.Preserves i := by + simp +decide only [] + exact hp i hi) + have hinv : inverseBase L m N ≤ catalogueBase L m N := by + unfold catalogueBase tableBase + have := Nat.zero_le (2 ^ L) + omega + have hsnapshot (a : ℕ) (ha : a < inputEnd m N) : + t.mem (snapshotBase L m N + a) = x.mem a := by + rw [htframe _ (by unfold inverseBase; omega) (by unfold tableBase inverseBase; omega), + hcmem, hsmem, ArrayCopy.copyPrefix_read ha, hhmem, Nat.zero_add] + have hmatrix (side : Bool) (r : ℕ) (hr : r < N) (k : ℕ) (hk : k < D m) : + t.mem (matrixAddress L m N side r k) = encode (matrix side r k) := by + have hb := matrixAddress_lt_inverse (L := L) side hr hk + unfold matrixAddress inverseBase at hb + have ho : matrixOffset m N side r k + 4 < inputEnd m N := by omega + rw [matrixAddress, hsnapshot _ ho] + exact hmat side r hr k hk + obtain ⟨l, hel, hlp, hlvalues, hlframe⟩ := exec_encodeBands hw false t htp + hlookup (matrix false) (hmatrix false) A hA (ha false) hmag hfit + have hlprotected := band_frame_protected hlframe + have hinverse (a : ℕ) (ha : a < 2 ^ L) : + snapshotBase L m N ≤ inverseBase L m N + a ∧ + inverseBase L m N + a < catalogueBase L m N := by + have hs : snapshotBase L m N ≤ inverseBase L m N := Nat.le_add_right _ _ + have ht : inverseBase L m N + a < tableBase L m N := by + change inverseBase L m N + a < inverseBase L m N + 2 ^ L + omega + exact ⟨hs.trans (Nat.le_add_right _ _), ht.trans_le (Nat.le_add_right _ _)⟩ + have hlookup' (a : ℕ) (ha : a < 2 ^ L) : + l.mem (inverseBase L m N + a) = MaskTable.lookup L m a := by + rw [hlprotected _ (hinverse a ha)] + exact hlookup a ha + have hmatrix' (r : ℕ) (hr : r < N) (k : ℕ) (hk : k < D m) : + l.mem (matrixAddress L m N true r k) = encode (matrix true r k) := by + have hb := matrixAddress_lt_inverse (L := L) true hr hk + rw [hlprotected (matrixAddress L m N true r k) ⟨Nat.le_add_right _ _, hb.trans_le hinv⟩] + exact hmatrix true r hr k hk + obtain ⟨r, her, hrp, hrvalues, hrframe⟩ := exec_encodeBands hw true l hlp + hlookup' (matrix true) hmatrix' A hA (ha true) hmag hfit + have hrprotected := band_frame_protected hrframe + obtain ⟨y, hey, hymem, hyr⟩ := exec_restoreInput r hrp.toHeaderPost hfit + have hyp : BandParameters L m N y := parameters_reframe hrp (by + intro i hi + have hp : ∀ i ∈ parameterRegisters, 7 ≤ i.val := by decide + exact hyr i (hp i hi)) + have hrestore (a : ℕ) (ha : inputEnd m N ≤ a) : y.mem a = r.mem a := by + rw [hymem, ArrayCopy.copyPrefix_outside (by omega)] + have hprotected (a : ℕ) (ha : snapshotBase L m N ≤ a ∧ a < catalogueBase L m N) : + r.mem a = t.mem a := (hrprotected a ha).trans (hlprotected a ha) + have hcatlow : Encoding.scratchBase L + 10 ^ L ≤ catalogueBase L m N := + (workspace_le_snapshot L m N).trans (snapshot_le_catalogue L m N) + have hinputcat : inputEnd m N ≤ catalogueBase L m N := + (input_le_snapshot L m N).trans (snapshot_le_catalogue L m N) + refine ⟨y, ?_, hyp, ?_, ?_, ?_, ?_⟩ + · apply (heh.seq (hes.seq (hec.seq (het.seq (hel.seq (her.seq hey)))))).mono + exact execution_time_le L m N + · intro side β hβ p hp + have hlow : inputEnd m N ≤ catalogueAddress L m N side β p := by + have := catalogue_le_address L m N β side + unfold catalogueAddress at * + omega + rw [hrestore _ hlow] + cases side + · rw [hrframe _ ?_ ?_] + · exact hlvalues β hβ p hp + · unfold catalogueAddress + simp only [Bool.false_eq_true, ↓reduceIte] + omega + · have hlt : catalogueAddress L m N false β p < + catalogueAddress L m N true 0 0 := by + unfold catalogueAddress + simp only [Bool.false_eq_true, ↓reduceIte, Nat.zero_add, Nat.add_zero] + have ht : 0 < 10 ^ L := Nat.pow_pos (by decide) + nlinarith only [hβ, hp, ht] + omega + · exact hrvalues β hβ p hp + · intro a ha + have hb := hinverse a ha + rw [hrestore _ ((input_le_snapshot L m N).trans hb.1), hprotected _ hb] + exact hlookup a ha + · intro gh + have hK : K L m ≤ 2 ^ L := Nat.choose_le_two_pow _ _ + have hi := (MaskTable.pairRank gh).isLt + have hb : snapshotBase L m N ≤ tableBase L m N + (MaskTable.pairRank gh).val ∧ + tableBase L m N + (MaskTable.pairRank gh).val < catalogueBase L m N := by + unfold catalogueBase tableBase inverseBase + omega + rw [hrestore _ ((input_le_snapshot L m N).trans hb.1), hprotected _ hb] + exact hforward gh + · intro b hb + by_cases hbin : b < inputEnd m N + · have hread : y.mem b = r.mem (snapshotBase L m N + b) := by + rw [hymem] + simpa only [Nat.zero_add] using (ArrayCopy.copyPrefix_read + (mem := r.mem) (src := snapshotBase L m N) (dst := 0) hbin) + rw [hread, hprotected _ (by + unfold catalogueBase tableBase inverseBase + have := Nat.zero_le (2 ^ L) + omega)] + exact hsnapshot b hbin + · have hbig : memoryEnd L m N ≤ b := by omega + rw [hrestore b (by omega), band_frame_above hrframe b hbig, + band_frame_above hlframe b hbig, htframe b ?_ ?_, hcmem, hsmem, + ArrayCopy.copyPrefix_outside ?_, hhmem] + · have hend := inverse_le_memoryEnd L m N + change snapshotBase L m N + inputEnd m N ≤ memoryEnd L m N at hend + omega + · have hend := table_le_memoryEnd L m N + change inverseBase L m N + 2 ^ L ≤ memoryEnd L m N at hend + omega + · have hend := catalogue_le_memoryEnd L m N + change tableBase L m N + 2 ^ L ≤ memoryEnd L m N at hend + omega + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandFlat.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandFlat.lean new file mode 100644 index 0000000000..1b356a8f0e --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandFlat.lean @@ -0,0 +1,118 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandInput + +/-! +# Flat native preparation of every public band encoding + +## Main results + +* `BandPreparation.runs`: actual flat `State.init` run, all public left/right encodings for + the computed deterministic layout, retained computed tables, input/outside frame and charged halt. +* `BandPreparation.program_length_le`: one fixed finite program, independent of every parameter. + +This is a bounded preparation component under the ordinary native 3SUM objective, not a +whole 3SUM algorithm or a pruned-query evaluation theorem. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The left catalogue retains the original mathematical band and leaf interfaces. -/ +theorem catalogueValue_left {L m N : ℕ} (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (matrix : ℕ → ℕ → ℤ) + (hX : ∀ i k, matrix i.val k.val = X i k) (β : ℕ) (τ : Leaf L) : + catalogueValue L m N false matrix β (Encoding.leafIndex τ).val = + encodingL (bandArrayL (MaskTable.layout hm) X β) τ := by + apply Encoding.value_final_left + intro u + simpa only [Equiv.symm_apply_apply] using + (bandValue_left hm X matrix hX β (finFunctionFinEquiv u)) + +/-- The right catalogue uses the same original flat leaf order. -/ +theorem catalogueValue_right {L m N : ℕ} (hm : m ≤ L) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (matrix : ℕ → ℕ → ℤ) + (hY : ∀ k i, matrix i.val k.val = Y k i) (β : ℕ) (τ : Leaf L) : + catalogueValue L m N true matrix β (Encoding.leafIndex τ).val = + encodingR (bandArrayR (MaskTable.layout hm) Y β) τ := by + apply Encoding.value_final_right + intro u + simpa only [Equiv.symm_apply_apply] using + (bandValue_right hm Y matrix (fun i k => hY k i) β (finFunctionFinEquiv u)) + +set_option maxRecDepth 4096 in +-- Static size arithmetic traverses the new finite blocks, reusing both encoder size bounds. +/-- One finite code-size bound for the 96-register program and its final halt. -/ +theorem program_length_le : (program.program catalogue).length ≤ 16384 := by + have hl := Encoding.program_length_le false + have hr := Encoding.program_length_le true + simp only [Code.program, List.length_append, Code.length_assemble, + List.length_cons, List.length_nil] at hl hr ⊢ + simp only [program, header, headerLoads, headerPowers, headerLayout, powerCode, + saveInput, countCode, countMasks, countBody, populationCode, rootCode, rootBody, + bandCount, tablesCode, tableBody, tableWrite, encodeBands, bandBody, prepareBand, + prepareHead, prepareBody, splitCode, splitHead, splitBody, readCell, collectBand, + collectHead, restoreInput, Code.size, Code.size_rename, Power.program, Power.body, + ArrayCopy.program, ArrayCopy.body, Bool.false_eq_true, ↓reduceIte, Nat.reduceAdd] + omega + +set_option maxRecDepth 4096 in +-- This finite structural check covers the same fixed instruction list, not an execution. +/-- The uniform program does not inspect its target word width. -/ +theorem not_readsWidth : ¬program.ReadsWidth := by decide + +/-- Every row and column band is encoded exactly once into its shared flat block. The +actual input contains only the three scalar parameters and signed matrices. `w` is the signed +source magnitude width and `w+2` is the target width used by the existing native kernels. +The prefix run exposes the actual final memory; its next instruction is the charged halt. +No unrelated caller memory or header is required to be zero in the general `exec` theorem. -/ +theorem runs {w L m N : ℕ} (hw : 0 < w) (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (A : ℤ) (hA : 0 ≤ A) (hX : ∀ i k, |X i k| ≤ A) (hY : ∀ k j, |Y k j| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + InputFits (w + 2) (input L m N X Y) ∧ + ∃ y n, n ≤ timeBound L m N ∧ + Runs (w + 2) (program.program catalogue) (State.init (input L m N X Y)) + ⟨program.size, y⟩ n ∧ + Executes (w + 2) (program.program catalogue) (State.init (input L m N X Y)) + (catalogueBase L m N) (n + 1) ∧ + y.Bounded (w + 2) ∧ + (∀ β < numBands L m N, ∀ τ, + y.mem (catalogueAddress L m N false β (Encoding.leafIndex τ).val) = + encode (encodingL (bandArrayL (MaskTable.layout hm) X β) τ)) ∧ + (∀ β < numBands L m N, ∀ τ, + y.mem (catalogueAddress L m N true β (Encoding.leafIndex τ).val) = + encode (encodingR (bandArrayR (MaskTable.layout hm) Y β) τ)) ∧ + (∀ a < 2 ^ L, y.mem (inverseBase L m N + a) = MaskTable.lookup L m a) ∧ + (∀ gh : Fin (K0 L m) × Fin (K0 L m), + y.mem (tableBase L m N + (MaskTable.pairRank gh).val) = + (MaskTable.selectedMask gh).val) ∧ + (∀ b, ¬(inputEnd m N ≤ b ∧ b < memoryEnd L m N) → + y.mem b = (State.init (input L m N X Y) : State 96).data.mem b) := by + obtain ⟨hi, y, he, hb, hp, hv, ht, hf, hframe⟩ := + exec_input hw hm X Y A hA hX hY hmag hfit + obtain ⟨n, hn, he⟩ := he + refine ⟨hi, y, n, hn, ?_, ?_, hb, ?_, ?_, ht, hf, hframe⟩ + · simpa only [Nat.zero_add, State.init] using he.runs (pc := 0) + (p := program.program catalogue) (Program.Contains.self _).append_left + · simpa only [hp.catalogue, State.init] using he.executes catalogue + · intro β hβ τ + rw [hv false β hβ _ (Encoding.leafIndex τ).isLt, + catalogueValue_left hm X _ (by intro i k; simp [matrixValues])] + · intro β hβ τ + rw [hv true β hβ _ (Encoding.leafIndex τ).isLt, + catalogueValue_right hm Y _ (by intro k i; simp [matrixValues])] + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandHeader.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandHeader.lean new file mode 100644 index 0000000000..abc4fe7e9c --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandHeader.lean @@ -0,0 +1,278 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram + +/-! +# Actual header loads, native powers and concrete band addresses + +## Main results + +* `BandPreparation.exec_powerCode`: exact reuse of the five-register native power routine. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Renaming the actual power routine has no instruction overhead and changes only its +accumulator, counter and constant-one register. -/ +theorem exec_powerCode {w : ℕ} (v b n : Fin 96) + (hinj : Function.Injective (powerRegisters v b n)) (B M : ℕ) (hB : 1 ≤ B) + (x : Data 96) (hb : x.regs b = B) (hn : x.regs n = M) + (hfit : Fits w (B ^ M + M + 2)) : + ∃ y, ExecWithin w (powerCode v b n) x y (4 * M + 5) ∧ + y.regs v = B ^ M ∧ y.regs one = 1 ∧ y.mem = x.mem ∧ + (∀ i, i ≠ v → i ≠ counter → i ≠ one → y.regs i = x.regs i) := by + obtain ⟨z, hz, hv, _, ho, hm, _, _⟩ := + Power.exec B M hB (x.comap (powerRegisters v b n)) hb hn hfit + have he := hz.rename hinj + refine ⟨z.extend (powerRegisters v b n) x, he, ?_, ?_, hm, ?_⟩ + · have h := congrArg (fun d : Data 5 => d.regs Power.value) (Data.comap_extend hinj z x) + exact h.trans hv + · have h := congrArg (fun d : Data 5 => d.regs Power.one) (Data.comap_extend hinj z x) + exact h.trans ho + · intro i hv hc ho + apply he.regs_eq + simp [Power.program, Power.body, Code.rename, Code.Preserves, Op.rename, Op.dest, Op.bin, + powerRegisters, Power.value, Power.counter, Power.base, Power.exponent, Power.one, + Ne.symm hv, Ne.symm hc, Ne.symm ho] + +/-- Scalar parameters loaded from the three actual input headers and small constants. -/ +structure BasicParameters (L m N : ℕ) (x : Data 96) : Prop where + /-- Input depth. -/ + levels : x.regs levels = L + /-- Input inner cardinality. -/ + inner : x.regs inner = m + /-- Input matrix dimension. -/ + rows : x.regs rows = N + /-- Constant zero. -/ + zero : x.regs zero = 0 + /-- Constant one. -/ + one : x.regs one = 1 + /-- Constant two. -/ + two : x.regs two = 2 + /-- Constant three. -/ + three : x.regs three = 3 + /-- Constant four. -/ + four : x.regs four = 4 + /-- Constant seven. -/ + seven : x.regs seven = 7 + /-- Constant ten. -/ + ten : x.regs ten = 10 + +private def basicRegisters : Finset (Fin 96) := + {levels, inner, rows, zero, two, three, four, seven, ten} + +private theorem BasicParameters.of_frame {L m N : ℕ} {x y : Data 96} + (hx : BasicParameters L m N x) (hone : y.regs BandPreparation.one = 1) + (h : ∀ i ∈ basicRegisters, y.regs i = x.regs i) : BasicParameters L m N y := + ⟨(h _ (by decide)).trans hx.levels, (h _ (by decide)).trans hx.inner, + (h _ (by decide)).trans hx.rows, (h _ (by decide)).trans hx.zero, hone, + (h _ (by decide)).trans hx.two, (h _ (by decide)).trans hx.three, + (h _ (by decide)).trans hx.four, (h _ (by decide)).trans hx.seven, + (h _ (by decide)).trans hx.ten⟩ + +private theorem BasicParameters.after_power {L m N : ℕ} {v : Fin 96} {x y : Data 96} + (hx : BasicParameters L m N x) (hv : v ∉ basicRegisters) + (hone : y.regs BandPreparation.one = 1) + (hf : ∀ i, i ≠ v → i ≠ counter → i ≠ BandPreparation.one → y.regs i = x.regs i) : + BasicParameters L m N y := hx.of_frame hone (by + intro i hi + have hstatic : ∀ i ∈ basicRegisters, i ≠ counter ∧ i ≠ BandPreparation.one := by decide + exact hf i (fun heq => hv (heq ▸ hi)) (hstatic i hi).1 (hstatic i hi).2) + +/-- The global numerical budget bounds every final header power and its exponent. -/ +theorem power_fits {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) + (b n : ℕ) (hb : b ≤ 10) (hn : n ≤ L) : Fits w (b ^ n + n + 2) := by + have hp : b ^ n ≤ 10 ^ L := + (Nat.pow_le_pow_left hb _).trans (Nat.pow_le_pow_right (by decide) hn) + have hpos : 1 ≤ 10 ^ L := Nat.one_le_pow L 10 (by decide) + apply hfit.of_le + unfold wordBudget + nlinarith only [hp, hn, hpos, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +/-- All header powers are computed by five actual native loops. Their register embeddings +are injective and their complete cost is linear in the original depth. -/ +theorem exec_headerPowers {w L m N : ℕ} (hm : m ≤ L) (x : Data 96) + (hx : BasicParameters L m N x) (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w headerPowers x y (16 * L + 26) ∧ BasicParameters L m N y ∧ + y.regs columns = D m ∧ y.regs outerSize = N0 L m ∧ y.regs masks = 2 ^ L ∧ + y.regs cells = 7 ^ L ∧ y.regs leaves = 10 ^ L ∧ y.mem = x.mem := by + obtain ⟨x1, h1, hD, ho1, hm1, hf1⟩ := exec_powerCode columns four inner (by decide) + 4 m (by decide) x hx.four hx.inner (power_fits hfit 4 m (by decide) hm) + have hp1 : BasicParameters L m N x1 := hx.after_power (by decide) ho1 hf1 + let x2 := x1.set temporary (L - m) + have hsub : ExecWithin w (ram[ temporary ← levels - inner ]) x1 x2 1 := by + have hL := power_fits hfit 1 L (by decide) le_rfl + have hs : Fits w (L - m) := hL.of_le (by simp only [one_pow]; omega) + ram_eval [x2, hp1.levels, hp1.inner, hs, temporary, levels, inner] + have ho2 : x2.regs one = 1 := by simpa [x2, temporary, one] using hp1.one + have hp2 : BasicParameters L m N x2 := hp1.of_frame ho2 (by + intro i hi + apply hsub.regs_eq + have hf : ∀ i ∈ basicRegisters, (ram[ temporary ← levels - inner ]).Preserves i := by + simp only [Code.Preserves, Op.dest_bin] + decide + exact hf i hi) + have htemp : x2.regs temporary = L - m := by simp [x2] + obtain ⟨x3, h3, hN0, ho3, hm3, hf3⟩ := exec_powerCode outerSize three temporary (by decide) + 3 (L - m) (by decide) x2 hp2.three htemp + (power_fits hfit 3 (L - m) (by decide) (Nat.sub_le _ _)) + have hp3 : BasicParameters L m N x3 := hp2.after_power (by decide) ho3 hf3 + obtain ⟨x4, h4, hS2, ho4, hm4, hf4⟩ := exec_powerCode masks two levels (by decide) + 2 L (by decide) x3 hp3.two hp3.levels (power_fits hfit 2 L (by decide) le_rfl) + have hp4 : BasicParameters L m N x4 := hp3.after_power (by decide) ho4 hf4 + obtain ⟨x5, h5, hS7, ho5, hm5, hf5⟩ := exec_powerCode cells seven levels (by decide) + 7 L (by decide) x4 hp4.seven hp4.levels (power_fits hfit 7 L (by decide) le_rfl) + have hp5 : BasicParameters L m N x5 := hp4.after_power (by decide) ho5 hf5 + obtain ⟨y, h6, hS10, ho6, hm6, hf6⟩ := exec_powerCode leaves ten levels (by decide) + 10 L (by decide) x5 hp5.ten hp5.levels (power_fits hfit 10 L (by decide) le_rfl) + refine ⟨y, ?_, hp5.after_power (by decide) ho6 hf6, ?_, ?_, ?_, ?_, hS10, ?_⟩ + · apply (h1.seq (hsub.seq (h3.seq (h4.seq (h5.seq h6))))).mono + omega + · have hD2 : x2.regs columns = 4 ^ m := by simpa [x2, temporary, columns] using hD + exact (hf6 columns (by decide) (by decide) (by decide)).trans + ((hf5 columns (by decide) (by decide) (by decide)).trans + ((hf4 columns (by decide) (by decide) (by decide)).trans + ((hf3 columns (by decide) (by decide) (by decide)).trans hD2))) + · exact (hf6 outerSize (by decide) (by decide) (by decide)).trans + ((hf5 outerSize (by decide) (by decide) (by decide)).trans + ((hf4 outerSize (by decide) (by decide) (by decide)).trans hN0)) + · exact (hf6 masks (by decide) (by decide) (by decide)).trans + ((hf5 masks (by decide) (by decide) (by decide)).trans hS2) + · exact (hf6 cells (by decide) (by decide) (by decide)).trans hS7 + · exact hm6.trans (hm5.trans (hm4.trans (hm3.trans hm1))) + +/-- Concrete computed powers and addresses after native header initialization. -/ +structure HeaderPost (L m N : ℕ) (x : Data 96) : Prop where + /-- Original scalar headers and constants. -/ + fixed : BasicParameters L m N x + /-- Thin dimension. -/ + columns : x.regs columns = D m + /-- Outer block dimension. -/ + outerSize : x.regs outerSize = N0 L m + /-- Complete mask scan range. -/ + masks : x.regs masks = 2 ^ L + /-- Flat input length of each band array. -/ + cells : x.regs cells = 7 ^ L + /-- Flat output length of each band encoding. -/ + leaves : x.regs leaves = 10 ^ L + /-- Original initialized prefix length. -/ + prefixSize : x.regs prefixSize = inputEnd m N + /-- Protected original-input copy address. -/ + snapshot : x.regs snapshot = snapshotBase L m N + /-- Dense inverse-table address. -/ + inverse : x.regs inverse = inverseBase L m N + /-- Selected-mask forward-table address. -/ + table : x.regs table = tableBase L m N + /-- Shared band-encoding output address. -/ + catalogue : x.regs catalogue = catalogueBase L m N + +/-- Header parameters are loaded from real memory, all five powers are computed natively, +and every concrete address calculation is charged. Header zero and all unrelated cells may +contain arbitrary words; initialization reads only addresses one through three. -/ +theorem exec_header {w L m N : ℕ} (hm : m ≤ L) (x : Data 96) + (hL : x.mem 1 = L) (hmemory : x.mem 2 = m) (hN : x.mem 3 = N) + (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w header x y (16 * L + 49) ∧ HeaderPost L m N y ∧ y.mem = x.mem := by + have h256 : 256 ≤ wordBudget L m N := Nat.le_mul_of_pos_right _ (Nat.succ_pos _) + have hsmall (v : ℕ) (hv : v ≤ 256) : Fits w v := hfit.of_le (hv.trans h256) + let x0 := (((((((((((((x.set address 1).set levels L).set address 2).set inner m).set + address 3).set rows N).set zero 0).set one 1).set two 2).set three 3).set four 4).set + seven 7).set ten 10) + have hs : ExecWithin w headerLoads x x0 13 := by + ram_eval [headerLoads, x0, hL, hmemory, hN, hsmall 0 (by decide), hsmall 1 (by decide), + hsmall 2 (by decide), hsmall 3 (by decide), hsmall 4 (by decide), + hsmall 7 (by decide), hsmall 10 (by decide), + address, levels, inner, rows, zero, one, two, three, four, seven, ten, Fin.mk.injEq] + have hp : BasicParameters L m N x0 := by + constructor <;> simp [x0, address, levels, inner, rows, zero, one, two, three, four, seven, ten] + obtain ⟨z, hz, hbasic, hD, hN0, hS2, hS7, hS10, hzmem⟩ := exec_headerPowers hm x0 hp hfit + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + have h2 : 2 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + have hlocal : Fits w (4 * (N * D m) + 8 * 10 ^ L + 16) := by + apply hfit.of_le + unfold wordBudget + nlinarith only [Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + let B := N * D m * 2 + 4 + let S := 10 ^ L * 2 + 7 ^ L + 2 + B + let I := S + B + let T := I + 2 ^ L + let C := T + 2 ^ L + have hnd : Fits w (N * D m) := hlocal.of_le (by omega) + have hnd2 : Fits w (N * D m * 2) := hlocal.of_le (by omega) + have hB : Fits w B := hlocal.of_le (by dsimp [B]; omega) + have ht2 : Fits w (10 ^ L * 2) := hlocal.of_le (by omega) + have hts : Fits w (10 ^ L * 2 + 7 ^ L) := hlocal.of_le (by omega) + have hts2 : Fits w (10 ^ L * 2 + 7 ^ L + 2) := hlocal.of_le (by omega) + have hS : Fits w S := hlocal.of_le (by dsimp [S, B]; omega) + have hI : Fits w I := hlocal.of_le (by dsimp [I, S, B]; omega) + have hT : Fits w T := hlocal.of_le (by dsimp [T, I, S, B]; omega) + have hC : Fits w C := hlocal.of_le (by dsimp [C, T, I, S, B]; omega) + let y := ((((((((((z.set prefixSize (N * D m)).set prefixSize (N * D m * 2)).set + prefixSize B).set snapshot (10 ^ L * 2)).set snapshot (10 ^ L * 2 + 7 ^ L)).set + snapshot (10 ^ L * 2 + 7 ^ L + 2)).set snapshot S).set inverse I).set table T).set catalogue C) + have hy : ExecWithin w headerLayout z y 10 := by + ram_eval [headerLayout, y, B, S, I, T, C, hbasic.rows, hbasic.two, hbasic.four, + hS2, hS7, hS10, hD, hnd, hnd2, hB, ht2, hts, hts2, hS, hI, hT, hC, + prefixSize, snapshot, inverse, table, catalogue, rows, columns, leaves, cells, + masks, two, four, Fin.mk.injEq, Function.update_idem] + have hfixed : BasicParameters L m N y := hbasic.of_frame (by + simpa [y, prefixSize, snapshot, inverse, table, catalogue, one] using hbasic.one) (by + intro i hi + apply hy.regs_eq + have hf : ∀ i ∈ basicRegisters, headerLayout.Preserves i := by + simp only [headerLayout, Code.Preserves, Op.dest_bin] + decide + exact hf i hi) + refine ⟨y, ?_, ?_, hzmem⟩ + · exact (hs.seq (hz.seq hy)).mono (by omega) + · constructor + · exact hfixed + · simpa [y, prefixSize, snapshot, inverse, table, catalogue, columns] using hD + · simpa [y, prefixSize, snapshot, inverse, table, catalogue, outerSize] using hN0 + · simpa [y, prefixSize, snapshot, inverse, table, catalogue, masks] using hS2 + · simpa [y, prefixSize, snapshot, inverse, table, catalogue, cells] using hS7 + · simpa [y, prefixSize, snapshot, inverse, table, catalogue, leaves] using hS10 + · have hb : y.regs prefixSize = B := by + simp +decide [y, prefixSize, snapshot, inverse, table, catalogue] + rw [hb] + dsimp only [B, inputEnd] + omega + · have hs : y.regs snapshot = S := by + simp +decide [y, prefixSize, snapshot, inverse, table, catalogue] + rw [hs] + dsimp only [S, B, snapshotBase, inputEnd, Encoding.scratchBase, Encoding.outputBase] + omega + · have hi : y.regs inverse = I := by + simp +decide [y, prefixSize, snapshot, inverse, table, catalogue] + rw [hi] + dsimp only [I, S, B, inverseBase, snapshotBase, inputEnd, + Encoding.scratchBase, Encoding.outputBase] + omega + · have ht : y.regs table = T := by + simp +decide [y, prefixSize, snapshot, inverse, table, catalogue] + rw [ht] + dsimp only [T, I, S, B, tableBase, inverseBase, snapshotBase, inputEnd, + Encoding.scratchBase, Encoding.outputBase] + omega + · have hc : y.regs catalogue = C := by simp [y] + rw [hc] + dsimp only [C, T, I, S, B, catalogueBase, tableBase, inverseBase, snapshotBase, + inputEnd, Encoding.scratchBase, Encoding.outputBase] + omega + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandInput.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandInput.lean new file mode 100644 index 0000000000..5a60b0da5f --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandInput.lean @@ -0,0 +1,149 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandExec +public import Complexity.Foundations.Data.List.RowMajor +import Mathlib.Data.List.GetD + +/-! +# Actual initialized header and signed-matrix inputs + +## Main definitions + +* `BandPreparation.matrixValues`: zero extension of the two original finite matrices for + stating the already proved native guarded-load specification. + +## Main results + +* `BandPreparation.input_fits`: the actual initialized list fits from original magnitude + and parameter bounds. +* `BandPreparation.exec_input`: complete charged native execution from `State.init`. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Original matrix entries, zero extended to natural indices. Right-side `r` denotes a +column of `Y`. This is only the mathematical memory specification; the program loads the +concrete row-major input words, after its actual bounds guard. -/ +def matrixValues {m N : ℕ} (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (side : Bool) (r k : ℕ) : ℤ := + if hr : r < N then if hk : k < D m then + if side then Y ⟨k, hk⟩ ⟨r, hr⟩ else X ⟨r, hr⟩ ⟨k, hk⟩ else 0 else 0 + +/-- The input prefix consists of the native length header, three scalar words and both +original matrices, with no powers, tables or precomputed outputs. -/ +theorem input_length {m N : ℕ} (L : ℕ) (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) : (input L m N X Y).length + 1 = inputEnd m N := by + simp only [input, List.flatMap_ofFn_prod, List.length_append, List.length_ofFn, + List.length_cons, List.length_nil, inputEnd] + ring + +private theorem init_matrix {L m N : ℕ} (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (side : Bool) (r : ℕ) (hr : r < N) + (k : ℕ) (hk : k < D m) : + (State.init (input L m N X Y) : State 96).data.mem (matrixOffset m N side r k + 4) = + encode (matrixValues X Y side r k) := by + let xs := (List.finRange N).flatMap (fun i => List.ofFn fun k => encode (X i k)) + let ys := (List.finRange (D m)).flatMap (fun k => List.ofFn fun j => encode (Y k j)) + have hxs : xs.length = N * D m := by + simp only [xs, List.flatMap_ofFn_prod, List.length_ofFn] + change (Data.init (input L m N X Y) : Data 96).mem _ = _ + rw [show matrixOffset m N side r k + 4 = (matrixOffset m N side r k + 3) + 1 by omega, + Data.init_mem_succ] + change ((L :: m :: N :: (xs ++ ys)).getD (matrixOffset m N side r k + 3) 0) = _ + rw [show matrixOffset m N side r k + 3 = + ((matrixOffset m N side r k + 1) + 1) + 1 by omega] + simp only [List.getD_cons_succ] + cases side + · have hi : r * D m + k < N * D m := by + have h := (finProdFinEquiv (⟨r, hr⟩, ⟨k, hk⟩)).isLt + simpa only [finProdFinEquiv, Equiv.coe_fn_mk, Nat.add_comm, Nat.mul_comm] using h + simp only [matrixOffset, Bool.false_eq_true, ↓reduceIte] + rw [List.getD_append _ _ _ _ (by simpa only [hxs] using hi)] + simpa only [xs, matrixValues, dite_eq_left hr, dite_eq_left hk, + Bool.false_eq_true, ↓reduceIte] using + (List.getD_flatMap_ofFn (fun i k => encode (X i k)) 0 ⟨r, hr⟩ ⟨k, hk⟩) + · simp only [matrixOffset, ↓reduceIte] + rw [List.getD_append_right _ _ _ _ (by rw [hxs]; omega), hxs, + Nat.add_sub_cancel_right] + simpa only [ys, matrixValues, dite_eq_left hr, dite_eq_left hk, ↓reduceIte] using + (List.getD_flatMap_ofFn (fun k j => encode (Y k j)) 0 ⟨k, hk⟩ ⟨r, hr⟩) + +/-- The sign/magnitude input encoding uses at most `w+1` bits; the actual native program +uses `w+2` bits for its signed arithmetic kernels. Every input word is covered. -/ +theorem input_fits {w L m N : ℕ} (hw : 0 < w) (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (A : ℤ) (hA : 0 ≤ A) (hX : ∀ i k, |X i k| ≤ A) (hY : ∀ k j, |Y k j| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : InputFits (w + 2) (input L m N X Y) := by + have hscalar : Fits (w + 2) (L + N) := (local_fits hfit).of_le (by + nlinarith only [Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)]) + have hpow : (1 : ℤ) ≤ (7 ^ L : ℕ) := by + exact_mod_cast Nat.one_le_pow L 7 (by decide) + have hAfit : A < (2 : ℤ) ^ (w - 1) := by nlinarith only [hpow, hA, hmag] + have hvalue (z : ℤ) (hz : |z| ≤ A) : Fits (w + 2) (encode z) := by + have hb := abs_lt.mp (hz.trans_lt hAfit) + exact (encode_fits (show GenericRAM.Signed.Fits w z from ⟨hw, hb.1.le, hb.2⟩)).mono + (by omega) + refine ⟨?_, ?_⟩ + · apply ((memory_fits hfit).of_le (input_le_memoryEnd L m N)).of_le + have := input_length L X Y + omega + · intro v hv + simp only [input, List.mem_append, List.mem_cons, List.not_mem_nil, or_false] at hv + rcases hv with (hv | hv) | hv + · rcases hv with rfl | rfl | rfl <;> exact hscalar.of_le (by omega) + · obtain ⟨i, _, hv⟩ := List.mem_flatMap.mp hv + obtain ⟨k, rfl⟩ := List.mem_ofFn.mp hv + exact hvalue _ (hX i k) + · obtain ⟨k, _, hv⟩ := List.mem_flatMap.mp hv + obtain ⟨j, rfl⟩ := List.mem_ofFn.mp hv + exact hvalue _ (hY k j) + +/-- Complete actual `State.init` execution, with both encoded catalogues, computed tables, +preserved input and outside frame, bounded final state and the complete charged time bound. -/ +theorem exec_input {w L m N : ℕ} (hw : 0 < w) (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (A : ℤ) (hA : 0 ≤ A) (hX : ∀ i k, |X i k| ≤ A) (hY : ∀ k j, |Y k j| ≤ A) + (hmag : (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget L m N)) : + InputFits (w + 2) (input L m N X Y) ∧ + ∃ y, ExecWithin (w + 2) program (State.init (input L m N X Y)).data y (timeBound L m N) ∧ + y.Bounded (w + 2) ∧ BandParameters L m N y ∧ + (∀ side β, β < numBands L m N → ∀ p, p < 10 ^ L → + y.mem (catalogueAddress L m N side β p) = + encode (catalogueValue L m N side (matrixValues X Y side) β p)) ∧ + (∀ a < 2 ^ L, y.mem (inverseBase L m N + a) = MaskTable.lookup L m a) ∧ + (∀ gh : Fin (K0 L m) × Fin (K0 L m), + y.mem (tableBase L m N + (MaskTable.pairRank gh).val) = + (MaskTable.selectedMask gh).val) ∧ + (∀ b, ¬(inputEnd m N ≤ b ∧ b < memoryEnd L m N) → + y.mem b = (State.init (input L m N X Y) : State 96).data.mem b) := by + have hi := input_fits hw hm X Y A hA hX hY hmag hfit + let x : Data 96 := (State.init (input L m N X Y)).data + have hL : x.mem 1 = L := by simp [x, State.init, Data.init, input] + have hm' : x.mem 2 = m := by simp [x, State.init, Data.init, input] + have hN : x.mem 3 = N := by simp [x, State.init, Data.init, input] + have hvalues (side : Bool) (r : ℕ) (hr : r < N) (k : ℕ) (hk : k < D m) : + |matrixValues X Y side r k| ≤ A := by + cases side <;> simp only [matrixValues, dite_eq_left hr, dite_eq_left hk, + Bool.false_eq_true, ↓reduceIte] + · exact hX _ _ + · exact hY _ _ + obtain ⟨y, he, hp, hv, ht, hf, hframe⟩ := exec hw hm x hL hm' hN + (matrixValues X Y) (init_matrix X Y) A hA hvalues hmag hfit + exact ⟨hi, y, he, he.bounded (Data.Bounded.init hi), hp, hv, ht, hf, hframe⟩ + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMeaning.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMeaning.lean new file mode 100644 index 0000000000..67daa9a4e0 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMeaning.lean @@ -0,0 +1,351 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.SplitDigits + +/-! +# Meaning of the native band indices + +## Main results + +* `BandPreparation.leftAt_isInner`, `BandPreparation.rightAt_isInner`: the original + encoder numberings classify the same inner positions as the native decoder. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open Encoding BandDigits + +/-- Digits three through six are exactly the inner left variables. -/ +theorem leftAt_isInner (d : Fin 7) : (leftAt d).IsInner ↔ ¬d.val < 3 := by decide +revert + +/-- The same four digits classify the inner right variables. -/ +theorem rightAt_isInner (d : Fin 7) : (rightAt d).IsInner ↔ ¬d.val < 3 := by decide +revert + +/-- An outer left digit is its ordinary base-three value. -/ +theorem leftAt_outer (d : Fin 7) (h : d.val < 3) : + (leftAt d).outerIndex.val = d.val := by decide +revert + +/-- An outer right digit is its ordinary base-three value. -/ +theorem rightAt_outer (d : Fin 7) (h : d.val < 3) : + (rightAt d).outerIndex.val = d.val := by decide +revert + +/-- The four inner left digits use the existing row-major pair numbering. -/ +theorem leftAt_inner (d : Fin 7) (h : ¬d.val < 3) : + (finProdFinEquiv (leftAt d).innerIndex).val = d.val - 3 := by decide +revert + +/-- The four inner right digits use the same pair numbering. -/ +theorem rightAt_inner (d : Fin 7) (h : ¬d.val < 3) : + (finProdFinEquiv (rightAt d).innerIndex).val = d.val - 3 := by decide +revert + +private theorem filter_outer {L : ℕ} (u : Fin L → Fin 7) + (Q : Finset (Fin L)) (hQ : ∀ ℓ, ℓ ∈ Q ↔ (u ℓ).val < 3) + (n : ℕ) (hn : Q.card = n) (f : Fin L → ℕ) + (hf : ∀ ℓ ∈ Q, f ℓ = (u ℓ).val) : + (List.ofFn fun ℓ => (u ℓ).val).filter (fun d => d < 3) = + List.ofFn (fun k => f (Q.orderEmbOfFin hn k)) := by + rw [Finset.ofFn_orderEmbOfFin, Finset.sort_eq_filter_finRange, + List.ofFn_eq_map, List.filter_map] + have hp : (fun ℓ : Fin L => decide ((u ℓ).val < 3)) = (fun ℓ => decide (ℓ ∈ Q)) := by + funext ℓ + simp only [hQ] + simp only [Function.comp_def, hp] + apply List.map_congr_left + intro ℓ hℓ + exact (hf ℓ (by simpa only [List.mem_filter, List.mem_finRange, true_and, + decide_eq_true_eq] using hℓ)).symm + +/-- Base-three compression follows the increasing outer positions used by `outerPartL`. -/ +theorem outer_digits_left {L m : ℕ} (u : Fin L → Fin 7) + (h : (innerSetL (fun ℓ => leftAt (u ℓ))).card = m) : + (List.ofFn fun ℓ => (u ℓ).val).filter (fun d => d < 3) = + List.ofFn (fun k => (outerPartL (fun ℓ => leftAt (u ℓ)) h k).val) := by + have hQ : ∀ ℓ, ℓ ∈ (innerSetL (fun ℓ => leftAt (u ℓ)))ᶜ ↔ (u ℓ).val < 3 := by + intro ℓ + simp [innerSetL, leftAt_isInner] + apply filter_outer u (innerSetL (fun ℓ => leftAt (u ℓ)))ᶜ hQ (L - m) + (Finset.card_compl_of_card_eq _ h) (fun ℓ => (leftAt (u ℓ)).outerIndex.val) + intro ℓ hℓ + exact leftAt_outer _ ((hQ ℓ).mp hℓ) + +/-- Base-three compression follows the increasing outer positions used by `outerPartR`. -/ +theorem outer_digits_right {L m : ℕ} (u : Fin L → Fin 7) + (h : (innerSetR (fun ℓ => rightAt (u ℓ))).card = m) : + (List.ofFn fun ℓ => (u ℓ).val).filter (fun d => d < 3) = + List.ofFn (fun k => (outerPartR (fun ℓ => rightAt (u ℓ)) h k).val) := by + have hQ : ∀ ℓ, ℓ ∈ (innerSetR (fun ℓ => rightAt (u ℓ)))ᶜ ↔ (u ℓ).val < 3 := by + intro ℓ + simp [innerSetR, rightAt_isInner] + apply filter_outer u (innerSetR (fun ℓ => rightAt (u ℓ)))ᶜ hQ (L - m) + (Finset.card_compl_of_card_eq _ h) (fun ℓ => (rightAt (u ℓ)).outerIndex.val) + intro ℓ hℓ + exact rightAt_outer _ ((hQ ℓ).mp hℓ) + +private theorem filter_inner {L : ℕ} (u : Fin L → Fin 7) + (Q : Finset (Fin L)) (hQ : ∀ ℓ, ℓ ∈ Q ↔ ¬(u ℓ).val < 3) + (n : ℕ) (hn : Q.card = n) (f : Fin L → ℕ) + (hf : ∀ ℓ ∈ Q, f ℓ = (u ℓ).val - 3) : + (((List.ofFn fun ℓ => (u ℓ).val).filter (fun d => ¬d < 3)).map fun d => d - 3) = + List.ofFn (fun k => f (Q.orderEmbOfFin hn k)) := by + rw [Finset.ofFn_orderEmbOfFin, Finset.sort_eq_filter_finRange, + List.ofFn_eq_map, List.filter_map, List.map_map] + have hp : (fun ℓ : Fin L => decide (¬(u ℓ).val < 3)) = (fun ℓ => decide (ℓ ∈ Q)) := by + funext ℓ + simp only [hQ] + simp only [Function.comp_def, hp] + apply List.map_congr_left + intro ℓ hℓ + exact (hf ℓ (by simpa only [List.mem_filter, List.mem_finRange, true_and, + decide_eq_true_eq] using hℓ)).symm + +/-- Base-four compression follows the increasing inner positions used by `innerPartL`. -/ +theorem inner_digits_left {L m : ℕ} (u : Fin L → Fin 7) + (h : (innerSetL (fun ℓ => leftAt (u ℓ))).card = m) : + (((List.ofFn fun ℓ => (u ℓ).val).filter (fun d => ¬d < 3)).map fun d => d - 3) = + List.ofFn (fun k => (finProdFinEquiv (innerPartL (fun ℓ => leftAt (u ℓ)) h k)).val) := by + have hQ : ∀ ℓ, ℓ ∈ innerSetL (fun ℓ => leftAt (u ℓ)) ↔ ¬(u ℓ).val < 3 := by + intro ℓ + simp [innerSetL, leftAt_isInner] + apply filter_inner u (innerSetL (fun ℓ => leftAt (u ℓ))) hQ m h + (fun ℓ => (finProdFinEquiv (leftAt (u ℓ)).innerIndex).val) + intro ℓ hℓ + exact leftAt_inner _ ((hQ ℓ).mp hℓ) + +/-- The right inner index uses the identical four pair digits. -/ +theorem inner_digits_right {L m : ℕ} (u : Fin L → Fin 7) + (h : (innerSetR (fun ℓ => rightAt (u ℓ))).card = m) : + (((List.ofFn fun ℓ => (u ℓ).val).filter (fun d => ¬d < 3)).map fun d => d - 3) = + List.ofFn (fun k => (finProdFinEquiv (innerPartR (fun ℓ => rightAt (u ℓ)) h k)).val) := by + have hQ : ∀ ℓ, ℓ ∈ innerSetR (fun ℓ => rightAt (u ℓ)) ↔ ¬(u ℓ).val < 3 := by + intro ℓ + simp [innerSetR, rightAt_isInner] + apply filter_inner u (innerSetR (fun ℓ => rightAt (u ℓ))) hQ m h + (fun ℓ => (finProdFinEquiv (rightAt (u ℓ)).innerIndex).val) + intro ℓ hℓ + exact rightAt_inner _ ((hQ ℓ).mp hℓ) + +/-- Native mask compression agrees with the exact inner level set of the original array. -/ +theorem parts_mask_left {L : ℕ} (a : Fin (7 ^ L)) : + (parts L a.val).mask = + (MaskTable.maskCode (innerSetL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ)))).val := by + rw [(parts_lists L a.val).1, sevenDigits_eq_ofFn] + change _ = (finFunctionFinEquiv (fun ℓ => + if ℓ ∈ innerSetL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ)) + then (1 : Fin 2) else 0)).val + rw [← Nat.ofDigits_finFunctionFinEquiv, List.map_ofFn] + congr 1 + apply congrArg List.ofFn + funext ℓ + change (if (finFunctionFinEquiv.symm a ℓ).val < 3 then 0 else 1) = _ + by_cases h : (finFunctionFinEquiv.symm a ℓ).val < 3 <;> + simp only [innerSetL, Finset.mem_filter, Finset.mem_univ, true_and, + leftAt_isInner, h, not_true_eq_false, not_false_eq_true, ↓reduceIte, + Fin.val_zero, Fin.val_one] + +/-- The right mask is classified identically, without changing the encoder numbering. -/ +theorem parts_mask_right {L : ℕ} (a : Fin (7 ^ L)) : + (parts L a.val).mask = + (MaskTable.maskCode (innerSetR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ)))).val := by + rw [(parts_lists L a.val).1, sevenDigits_eq_ofFn] + change _ = (finFunctionFinEquiv (fun ℓ => + if ℓ ∈ innerSetR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ)) + then (1 : Fin 2) else 0)).val + rw [← Nat.ofDigits_finFunctionFinEquiv, List.map_ofFn] + congr 1 + apply congrArg List.ofFn + funext ℓ + change (if (finFunctionFinEquiv.symm a ℓ).val < 3 then 0 else 1) = _ + by_cases h : (finFunctionFinEquiv.symm a ℓ).val < 3 <;> + simp only [innerSetR, Finset.mem_filter, Finset.mem_univ, true_and, + rightAt_isInner, h, not_true_eq_false, not_false_eq_true, ↓reduceIte, + Fin.val_zero, Fin.val_one] + +/-- Native left outer offset, with the exact public mathematical string interface. -/ +theorem parts_outer_left {L m : ℕ} (a : Fin (7 ^ L)) + (h : (innerSetL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ))).card = m) : + (parts L a.val).outer = + (finFunctionFinEquiv + (outerPartL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ)) h)).val := by + rw [(parts_lists L a.val).2.1, sevenDigits_eq_ofFn] + change Nat.ofDigits 3 ((List.ofFn fun ℓ => (finFunctionFinEquiv.symm a ℓ).val).filter + (fun d => d < 3)) = _ + rw [outer_digits_left _ h, Nat.ofDigits_finFunctionFinEquiv] + +/-- Native right outer offset in the same least-significant-first convention. -/ +theorem parts_outer_right {L m : ℕ} (a : Fin (7 ^ L)) + (h : (innerSetR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ))).card = m) : + (parts L a.val).outer = + (finFunctionFinEquiv + (outerPartR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ)) h)).val := by + rw [(parts_lists L a.val).2.1, sevenDigits_eq_ofFn] + change Nat.ofDigits 3 ((List.ofFn fun ℓ => (finFunctionFinEquiv.symm a ℓ).val).filter + (fun d => d < 3)) = _ + rw [outer_digits_right _ h, Nat.ofDigits_finFunctionFinEquiv] + +/-- Native left thin-dimension offset agrees with the row-major pair digits of `Layout`. -/ +theorem parts_inner_left {L m : ℕ} (a : Fin (7 ^ L)) + (h : (innerSetL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ))).card = m) : + (parts L a.val).inner = (finFunctionFinEquiv fun k => + finProdFinEquiv (innerPartL (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ)) h k)).val := by + rw [(parts_lists L a.val).2.2, sevenDigits_eq_ofFn] + change Nat.ofDigits 4 (((List.ofFn fun ℓ => (finFunctionFinEquiv.symm a ℓ).val).filter + (fun d => ¬d < 3)).map fun d => d - 3) = _ + rw [inner_digits_left _ h, Nat.ofDigits_finFunctionFinEquiv] + +/-- Native right thin-dimension offset uses the same `Layout` equivalence. -/ +theorem parts_inner_right {L m : ℕ} (a : Fin (7 ^ L)) + (h : (innerSetR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ))).card = m) : + (parts L a.val).inner = (finFunctionFinEquiv fun k => + finProdFinEquiv (innerPartR (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ)) h k)).val := by + rw [(parts_lists L a.val).2.2, sevenDigits_eq_ofFn] + change Nat.ofDigits 4 (((List.ofFn fun ℓ => (finFunctionFinEquiv.symm a ℓ).val).filter + (fun d => ¬d < 3)).map fun d => d - 3) = _ + rw [inner_digits_right _ h, Nat.ofDigits_finFunctionFinEquiv] + +/-- Numeric specification of one guarded band-array read. `matrix r k` denotes an original +matrix cell (or the transposed indexing of the right matrix); the native implementation +computes its actual row-major address and loads it. -/ +def bandValue (L m N β : ℕ) (rightSide : Bool) (matrix : ℕ → ℕ → ℤ) (a : ℕ) : ℤ := + let p := parts L a + let q := MaskTable.lookup L m p.mask + if q = 0 then 0 else + let b := if rightSide then (q - 1) % K0 L m else (q - 1) / K0 L m + let r := (β * K0 L m + b) * N0 L m + p.outer + if r < N then matrix r p.inner else 0 + +/-- The numeric left-band specification is the unchanged mathematical `bandArrayL` for +the computed deterministic layout, including all unselected strings and padding. -/ +theorem bandValue_left {L m N : ℕ} (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (matrix : ℕ → ℕ → ℤ) + (hX : ∀ i k, matrix i.val k.val = X i k) (β : ℕ) (a : Fin (7 ^ L)) : + bandValue L m N β false matrix a.val = + bandArrayL (MaskTable.layout hm) X β (fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ)) := by + let u : LeftStr L := fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ) + let lay := MaskTable.layout hm + have hmask : (parts L a.val).mask = (MaskTable.maskCode (innerSetL u)).val := + parts_mask_left a + by_cases hz : MaskTable.lookup L m (parts L a.val).mask = 0 + · have hmissing : ∀ gh, lay.table gh ≠ innerSetL u := by + intro gh heq + have hc := congrArg (fun Q => (MaskTable.maskCode Q).val) heq + change (MaskTable.maskCode (MaskTable.maskSet (MaskTable.selectedMask gh))).val = _ at hc + rw [MaskTable.maskCode_maskSet] at hc + have ha : (MaskTable.selectedMask gh).val = (parts L a.val).mask := hc.trans hmask.symm + have hne : MaskTable.lookup L m (MaskTable.selectedMask gh).val ≠ 0 := by + rw [MaskTable.lookup_selectedMask] + omega + exact hne (ha ▸ hz) + have hf := bandFamilyL_eq_zero lay X β hmissing + simp only [bandValue, ↓reduceIte, hz] + change 0 = (arrayL m (bandFamilyL lay X β)) u + unfold arrayL + split <;> simp only [hf, Matrix.zero_apply] + · obtain ⟨gh, hs⟩ := MaskTable.exists_selectedMask_of_lookup hz + have hcode : MaskTable.selectedMask gh = MaskTable.maskCode (innerSetL u) := + Fin.ext (hs.trans hmask) + have hset : lay.table gh = innerSetL u := by + change MaskTable.maskSet (MaskTable.selectedMask gh) = _ + rw [hcode, MaskTable.maskSet_maskCode] + have hcard : (innerSetL u).card = m := hset ▸ lay.table_card gh + have hf : bandFamilyL lay X β (innerSetL u) = rowBlock lay X (β * K0 L m + gh.1) := by + rw [← hset] + exact bandFamilyL_table lay X β gh.1 gh.2 + have ho : (parts L a.val).outer = (lay.rowIdx.symm (outerPartL u hcard)).val := + parts_outer_left a hcard + have hi : (parts L a.val).inner = (lay.innerIdx.symm (innerPartL u hcard)).val := by + change (parts L a.val).inner = + (finFunctionFinEquiv fun k => finProdFinEquiv (innerPartL u hcard k)).val + exact parts_inner_left a hcard + have hg : (MaskTable.lookup L m (parts L a.val).mask - 1) / K0 L m = gh.1.val := by + rw [← hs] + exact (MaskTable.lookup_blockPair gh).1 + change bandValue L m N β false matrix a.val = (arrayL m (bandFamilyL lay X β)) u + rw [arrayL, dite_eq_left hcard, hf] + simp only [bandValue, hz, ↓reduceIte, Bool.false_eq_true, hg, ho, hi, + rowBlock, padRows] + split_ifs with hr + · exact hX ⟨_, hr⟩ _ + · rfl + +/-- The numeric right-band specification is the unchanged mathematical `bandArrayR`. +Its natural row parameter refers to a column of the original right matrix. -/ +theorem bandValue_right {L m N : ℕ} (hm : m ≤ L) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (matrix : ℕ → ℕ → ℤ) + (hY : ∀ i k, matrix i.val k.val = Y k i) (β : ℕ) (a : Fin (7 ^ L)) : + bandValue L m N β true matrix a.val = + bandArrayR (MaskTable.layout hm) Y β (fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ)) := by + let u : RightStr L := fun ℓ => rightAt (finFunctionFinEquiv.symm a ℓ) + let lay := MaskTable.layout hm + have hmask : (parts L a.val).mask = (MaskTable.maskCode (innerSetR u)).val := + parts_mask_right a + by_cases hz : MaskTable.lookup L m (parts L a.val).mask = 0 + · have hmissing : ∀ gh, lay.table gh ≠ innerSetR u := by + intro gh heq + have hc := congrArg (fun Q => (MaskTable.maskCode Q).val) heq + change (MaskTable.maskCode (MaskTable.maskSet (MaskTable.selectedMask gh))).val = _ at hc + rw [MaskTable.maskCode_maskSet] at hc + have ha : (MaskTable.selectedMask gh).val = (parts L a.val).mask := hc.trans hmask.symm + have hne : MaskTable.lookup L m (MaskTable.selectedMask gh).val ≠ 0 := by + rw [MaskTable.lookup_selectedMask] + omega + exact hne (ha ▸ hz) + have hf := bandFamilyR_eq_zero lay Y β hmissing + simp only [bandValue, ↓reduceIte, hz] + change 0 = (arrayR m (bandFamilyR lay Y β)) u + unfold arrayR + split <;> simp only [hf, Matrix.zero_apply] + · obtain ⟨gh, hs⟩ := MaskTable.exists_selectedMask_of_lookup hz + have hcode : MaskTable.selectedMask gh = MaskTable.maskCode (innerSetR u) := + Fin.ext (hs.trans hmask) + have hset : lay.table gh = innerSetR u := by + change MaskTable.maskSet (MaskTable.selectedMask gh) = _ + rw [hcode, MaskTable.maskSet_maskCode] + have hcard : (innerSetR u).card = m := hset ▸ lay.table_card gh + have hf : bandFamilyR lay Y β (innerSetR u) = colBlock lay Y (β * K0 L m + gh.2) := by + rw [← hset] + exact bandFamilyR_table lay Y β gh.1 gh.2 + have ho : (parts L a.val).outer = (lay.colIdx.symm (outerPartR u hcard)).val := + parts_outer_right a hcard + have hi : (parts L a.val).inner = (lay.innerIdx.symm (innerPartR u hcard)).val := by + change (parts L a.val).inner = + (finFunctionFinEquiv fun k => finProdFinEquiv (innerPartR u hcard k)).val + exact parts_inner_right a hcard + have hg : (MaskTable.lookup L m (parts L a.val).mask - 1) % K0 L m = gh.2.val := by + rw [← hs] + exact (MaskTable.lookup_blockPair gh).2 + change bandValue L m N β true matrix a.val = (arrayR m (bandFamilyR lay Y β)) u + rw [arrayR, dite_eq_left hcard, hf] + simp only [bandValue, hz, ↓reduceIte, hg, ho, hi, colBlock, padCols] + split_ifs with hr + · exact hY ⟨_, hr⟩ _ + · rfl + +/-- A selected native mask ensures that its compressed offsets fit the actual matrix +dimensions, rather than merely the coarser depth-based bounds. -/ +theorem selected_parts_bounds {L m : ℕ} (a : Fin (7 ^ L)) + (hz : MaskTable.lookup L m (parts L a.val).mask ≠ 0) : + (parts L a.val).outer < N0 L m ∧ (parts L a.val).inner < D m := by + obtain ⟨gh, hs⟩ := MaskTable.exists_selectedMask_of_lookup hz + let u : LeftStr L := fun ℓ => leftAt (finFunctionFinEquiv.symm a ℓ) + have hc : MaskTable.selectedMask gh = MaskTable.maskCode (innerSetL u) := + Fin.ext (hs.trans (parts_mask_left a)) + have hset : MaskTable.maskSet (MaskTable.selectedMask gh) = innerSetL u := by + rw [hc, MaskTable.maskSet_maskCode] + have hcard : (innerSetL u).card = m := hset ▸ MaskTable.card_selectedMask gh + constructor + · rw [parts_outer_left a hcard] + exact (finFunctionFinEquiv (outerPartL u hcard)).isLt + · rw [parts_inner_left a hcard] + exact (finFunctionFinEquiv (fun k => finProdFinEquiv (innerPartL u hcard k))).isLt + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMemory.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMemory.lean new file mode 100644 index 0000000000..44f7f18870 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandMemory.lean @@ -0,0 +1,124 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram + +/-! +# Concrete disjoint memory regions and derived native fits + +## Main results + +* `BandPreparation.workspace_le_snapshot`: protected data are above the complete encoder. +* `BandPreparation.memory_fits`: the actual allocation end fits the explicit input budget. +* `BandPreparation.tables_fit`, `BandPreparation.split_fits`, `BandPreparation.encoder_fits`: + component inequalities derived from the same budget, without an intermediate-state oracle. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The native input prefix is disjoint from its protected copy. -/ +theorem input_le_snapshot (L m N : ℕ) : inputEnd m N ≤ snapshotBase L m N := by + unfold snapshotBase + exact (Nat.le_add_right (inputEnd m N) (Encoding.scratchBase L)).trans + (Nat.le_add_right _ (10 ^ L)) + +/-- The encoder's complete low workspace is below the protected original input. -/ +theorem workspace_le_snapshot (L m N : ℕ) : + Encoding.scratchBase L + 10 ^ L ≤ snapshotBase L m N := by + unfold snapshotBase + omega + +/-- The protected original-input copy ends exactly at the inverse table. -/ +theorem snapshot_end (L m N : ℕ) : + snapshotBase L m N + inputEnd m N = inverseBase L m N := rfl + +/-- The inverse table ends exactly at the selected-mask table. -/ +theorem inverse_end (L m N : ℕ) : inverseBase L m N + 2 ^ L = tableBase L m N := rfl + +/-- The selected-mask table ends exactly at the shared output catalogue. -/ +theorem table_end (L m N : ℕ) : tableBase L m N + 2 ^ L = catalogueBase L m N := rfl + +/-- The catalogue base is below the complete allocation end, including an empty catalogue. -/ +theorem catalogue_le_memoryEnd (L m N : ℕ) : catalogueBase L m N ≤ memoryEnd L m N := + Nat.le_add_right _ _ + +/-- Each table block precedes the complete allocation end. -/ +theorem table_le_memoryEnd (L m N : ℕ) : tableBase L m N ≤ memoryEnd L m N := + (Nat.le_add_right _ _ : tableBase L m N ≤ tableBase L m N + 2 ^ L).trans + (catalogue_le_memoryEnd L m N) + +/-- The protected prefix ends below every table/output allocation. -/ +theorem inverse_le_memoryEnd (L m N : ℕ) : inverseBase L m N ≤ memoryEnd L m N := + (Nat.le_add_right _ _ : inverseBase L m N ≤ inverseBase L m N + 2 ^ L).trans + (table_le_memoryEnd L m N) + +/-- The snapshot base is within the complete allocation. -/ +theorem snapshot_le_memoryEnd (L m N : ℕ) : snapshotBase L m N ≤ memoryEnd L m N := + (Nat.le_add_right _ _ : snapshotBase L m N ≤ snapshotBase L m N + inputEnd m N).trans + (inverse_le_memoryEnd L m N) + +/-- Every original input address also fits the complete allocation. -/ +theorem input_le_memoryEnd (L m N : ℕ) : inputEnd m N ≤ memoryEnd L m N := + (input_le_snapshot L m N).trans (snapshot_le_memoryEnd L m N) + +/-- The constant margin in the original budget covers all fixed program constants. -/ +theorem small_fits {w L m N v : ℕ} (hfit : Fits w (wordBudget L m N)) (hv : v ≤ 256) : + Fits w v := + hfit.of_le (hv.trans (Nat.le_mul_of_pos_right _ (Nat.succ_pos _))) + +/-- The entire real allocation, including both catalogue halves, fits the original numerical +word-budget inequality. The sizes of unreachable padded row-major products are irrelevant. -/ +theorem memory_fits {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) : + Fits w (memoryEnd L m N) := by + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + apply hfit.of_le + unfold memoryEnd catalogueBase tableBase inverseBase snapshotBase inputEnd + Encoding.scratchBase Encoding.outputBase wordBudget + nlinarith only [h7, Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +/-- Both table-store addresses and all population/count arithmetic fit before table execution. -/ +theorem tables_fit {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) : + Fits w (inverseBase L m N + tableBase L m N + 2 ^ L + L + 4) := by + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + apply hfit.of_le + unfold tableBase inverseBase snapshotBase inputEnd Encoding.scratchBase Encoding.outputBase + wordBudget + nlinarith only [h7, Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +/-- Every original band index admits the already proved native digit-loop fits. -/ +theorem split_fits {w L m N a : ℕ} (hfit : Fits w (wordBudget L m N)) (ha : a < 7 ^ L) : + Fits w (a + L + 16 * 10 ^ L + 16) := by + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + apply hfit.of_le + unfold wordBudget + nlinarith only [ha, h7, Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +/-- The unchanged native encoder's arithmetic/address premise follows from this client's +explicit original-parameter bound. Its signed magnitude premise remains separate. -/ +theorem encoder_fits {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) : + Fits w (Encoding.addressBound L) := by + apply hfit.of_le + unfold Encoding.addressBound wordBudget + have hpos : 1 ≤ 10 ^ L := Nat.one_le_pow L 10 (by decide) + nlinarith only [hpos, Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandParameters.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandParameters.lean new file mode 100644 index 0000000000..519cb3b84a --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandParameters.lean @@ -0,0 +1,152 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandHeader +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCount +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandRoot +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.BandBounds + +/-! +# Complete computed scalar parameters for native bands + +## Main results + +* `BandPreparation.exec_countCode`: compute binomial count, block count and band count + from actual scalar parameters, retaining the previously computed header layout. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Registers containing the persistent header and concrete layout. -/ +def headerRegisters : Finset (Fin 96) := + {levels, inner, rows, columns, outerSize, masks, cells, leaves, prefixSize, + snapshot, inverse, table, catalogue, zero, one, two, three, four, seven, ten} + +/-- The concrete layout survives any native execution that preserves its scalar registers. -/ +theorem headerPost_reframe {L m N : ℕ} {x y : Data 96} (hx : HeaderPost L m N x) + (h : ∀ i ∈ headerRegisters, y.regs i = x.regs i) : HeaderPost L m N y := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩ + · exact ⟨(h _ (by decide)).trans hx.fixed.levels, (h _ (by decide)).trans hx.fixed.inner, + (h _ (by decide)).trans hx.fixed.rows, (h _ (by decide)).trans hx.fixed.zero, + (h _ (by decide)).trans hx.fixed.one, (h _ (by decide)).trans hx.fixed.two, + (h _ (by decide)).trans hx.fixed.three, (h _ (by decide)).trans hx.fixed.four, + (h _ (by decide)).trans hx.fixed.seven, (h _ (by decide)).trans hx.fixed.ten⟩ + all_goals first + | exact (h _ (by decide)).trans hx.columns + | exact (h _ (by decide)).trans hx.outerSize + | exact (h _ (by decide)).trans hx.masks + | exact (h _ (by decide)).trans hx.cells + | exact (h _ (by decide)).trans hx.leaves + | exact (h _ (by decide)).trans hx.prefixSize + | exact (h _ (by decide)).trans hx.snapshot + | exact (h _ (by decide)).trans hx.inverse + | exact (h _ (by decide)).trans hx.table + | exact (h _ (by decide)).trans hx.catalogue + +/-- Persistent concrete layout and the additional scalar results actually computed by +the mask-count/root/band-count code. No fit or execution bound is bundled here. -/ +structure BandParameters (L m N : ℕ) (x : Data 96) : Prop extends HeaderPost L m N x where + /-- Computed floor square root of the binomial count. -/ + root : x.regs root = K0 L m + /-- Number of selected subset masks. -/ + square : x.regs square = K0 L m * K0 L m + /-- Computed number of shared bands per side. -/ + bands : x.regs bands = numBands L m N + +/-- All scalar parameters retained throughout preparation and encoder calls. -/ +def parameterRegisters : Finset (Fin 96) := headerRegisters ∪ {root, square, bands} + +/-- Preserve the computed parameters by preserving their actual registers. -/ +theorem parameters_reframe {L m N : ℕ} {x y : Data 96} (hx : BandParameters L m N x) + (h : ∀ i ∈ parameterRegisters, y.regs i = x.regs i) : BandParameters L m N y := by + refine ⟨headerPost_reframe hx.toHeaderPost (fun i hi => h i (Finset.mem_union_left _ hi)), + (h _ (by decide)).trans hx.root, (h _ (by decide)).trans hx.square, + (h _ (by decide)).trans hx.bands⟩ + +/-- A convenient local numerical bound, derived from the original global word budget. -/ +theorem local_fits {w L m N : ℕ} (hfit : Fits w (wordBudget L m N)) : + Fits w (16 * (L + N + N * D m + 2 ^ L + 10 ^ L + 1)) := by + apply hfit.of_le + unfold wordBudget + nlinarith only [Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L), + Nat.zero_le (2 * numBands L m N * 10 ^ L)] + +/-- Compute the root and band counts with the actual native loops. -/ +theorem exec_countCode {w L m N : ℕ} (x : Data 96) + (hx : HeaderPost L m N x) (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w countCode x y ((6 * L + 16) * 2 ^ L + 14) ∧ + BandParameters L m N y ∧ y.mem = x.mem := by + have hl := local_fits hfit + have hK : K L m ≤ 2 ^ L := Nat.choose_le_two_pow _ _ + obtain ⟨c, hc, hcount, _, hcmem, _⟩ := exec_countMasks L m x hx.fixed.levels hx.fixed.inner + hx.masks hx.fixed.one hx.fixed.two (hl.of_le (by + nlinarith only [Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)])) + have hcp : HeaderPost L m N c := headerPost_reframe hx (by + intro i hi + apply hc.regs_eq + have hp : ∀ i ∈ headerRegisters, countMasks.Preserves i := by + simp only [countMasks, countBody, populationCode, Code.Preserves, Op.dest, Op.bin] + decide + exact hp i hi) + obtain ⟨r, hr, hroot, hsquare, hrmem, _⟩ := exec_root (K L m) c hcount hcp.fixed.one + (hl.of_le (by + nlinarith only [hK, Nat.zero_le L, Nat.zero_le N, + Nat.zero_le (N * D m), Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)])) + have hrp : HeaderPost L m N r := headerPost_reframe hcp (by + intro i hi + apply hr.regs_eq + have hp : ∀ i ∈ headerRegisters, rootCode.Preserves i := by + simp only [rootCode, rootBody, Code.Preserves, Op.dest, Op.bin] + decide + exact hp i hi) + have hB : bandSize L m ≤ 7 ^ L := bandSize_le_seven_pow L m + have h7 : 7 ^ L ≤ 10 ^ L := Nat.pow_le_pow_left (by decide) _ + have hbfit : Fits w (bandSize L m) := hl.of_le (by + nlinarith only [hB, h7, Nat.zero_le L, Nat.zero_le N, + Nat.zero_le (N * D m), Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)]) + have hnfit : Fits w (N + bandSize L m) := hl.of_le (by + nlinarith only [hB, h7, Nat.zero_le L, Nat.zero_le N, + Nat.zero_le (N * D m), Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)]) + have hnsubfit : Fits w (N + bandSize L m - 1) := + hnfit.of_le (Nat.sub_le _ _) + have hdivfit : Fits w ((N + bandSize L m - 1) / bandSize L m) := + hnsubfit.of_le (Nat.div_le_self _ _) + let y := (((r.set temporary (bandSize L m)).set bands (N + bandSize L m)).set + bands (N + bandSize L m - 1)).set bands (numBands L m N) + have hroot0 : r.regs root = K0 L m := by simpa only [K0] using hroot + simp only [bandSize] at hbfit hnfit hnsubfit hdivfit + have hy : ExecWithin w bandCount r y 4 := by + ram_eval [bandCount, y, hroot0, hrp.outerSize, hrp.fixed.rows, hrp.fixed.one, + bandSize, numBands, hbfit, hnfit, hnsubfit, hdivfit, + temporary, root, outerSize, bands, rows, one, Fin.mk.injEq] + have he : ExecWithin w countCode x y ((6 * L + 16) * 2 ^ L + 14) := by + apply (hc.seq (hr.seq hy)).mono + have hsqrt : Nat.sqrt (K L m) ≤ 2 ^ L := (Nat.sqrt_le_self _).trans hK + nlinarith + have hyp : HeaderPost L m N y := headerPost_reframe hx (by + intro i hi + apply he.regs_eq + have hp : ∀ i ∈ headerRegisters, countCode.Preserves i := by + simp only [countCode, countMasks, countBody, populationCode, rootCode, rootBody, + bandCount, Code.Preserves, Op.dest, Op.bin] + decide + exact hp i hi) + refine ⟨y, he, ⟨hyp, ?_, ?_, ?_⟩, hrmem.trans hcmem⟩ + · simpa [y, temporary, bands, root] using hroot0 + · simpa [y, temporary, bands, square, K0] using hsquare + · simp [y] + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPopulation.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPopulation.lean new file mode 100644 index 0000000000..7e6315663c --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPopulation.lean @@ -0,0 +1,114 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.MaskTable + +/-! +# Charged native binary-mask population scans + +## Main results + +* `BandPreparation.exec_population`: the actual remainder/division loop computes the + population in `6L+5` instructions, preserves memory and changes only its four local registers. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open Cslib.Algorithms.MatrixMultiplication.Schonhage.MaskTable + +/-- Popcount is a native loop, with a fit derived for each arithmetic operation from a +bound on its original mask and depth. No memory cell is touched. -/ +theorem exec_population {w : ℕ} (L a : ℕ) (x : Data 96) + (hl : x.regs levels = L) (ha : x.regs mask = a) + (ho : x.regs one = 1) (ht : x.regs two = 2) (hfit : Fits w (a + L + 2)) : + ∃ y, ExecWithin w populationCode x y (6 * L + 5) ∧ + y.regs population = MaskTable.population L a ∧ y.mem = x.mem ∧ + (∀ i, i ≠ residue → i ≠ population → i ≠ counter → i ≠ digit → + y.regs i = x.regs i) := by + have hzero : Fits w 0 := hfit.of_le (by omega) + let start := ((x.set residue a).set population 0).set counter 0 + let I (y : Data 96) : Prop := + y.mem = x.mem ∧ y.regs levels = L ∧ y.regs one = 1 ∧ y.regs two = 2 ∧ + y.regs counter ≤ L ∧ y.regs residue ≤ a ∧ + y.regs population + MaskTable.population (L - y.regs counter) (y.regs residue) = + MaskTable.population L a + have hi : I start := by + simp [I, start, hl, ho, ht, counter, residue, population, levels, one, two] + have hb (y : Data 96) (hy : I y) (hc : (Cond.lt counter levels).Holds y) : + ∃ z, ExecWithin w (ram[ + digit ← residue % two; population ← population + digit; + residue ← residue / two; counter ← counter + one ]) y z 4 ∧ + I z ∧ L - z.regs counter < L - y.regs counter := by + obtain ⟨hm, hL, h1, h2, hj, hr, hp⟩ := hy + have hjL : y.regs counter < L := by simpa only [Cond.holds_lt, hL] using hc + have hrec : MaskTable.population (L - y.regs counter) (y.regs residue) = + y.regs residue % 2 + + MaskTable.population (L - (y.regs counter + 1)) (y.regs residue / 2) := by + rw [show L - y.regs counter = (L - (y.regs counter + 1)) + 1 by omega, + population_succ] + have hpop := population_le L a + have hfpop : Fits w (y.regs population + y.regs residue % 2) := + hfit.of_le (by omega) + have hfmod : Fits w (y.regs residue % 2) := hfit.of_le (by omega) + have hfdiv : Fits w (y.regs residue / 2) := + hfit.of_le ((Nat.div_le_self _ _).trans (by omega)) + have hfnext : Fits w (y.regs counter + 1) := hfit.of_le (by omega) + let z := (((y.set digit (y.regs residue % 2)).set population + (y.regs population + y.regs residue % 2)).set residue (y.regs residue / 2)).set + counter (y.regs counter + 1) + have he : ExecWithin w (ram[ + digit ← residue % two; population ← population + digit; + residue ← residue / two; counter ← counter + one ]) y z 4 := by + ram_eval [z, hfpop, hfmod, hfdiv, hfnext, h1, h2, + digit, residue, population, counter, one, two, Fin.mk.injEq] + refine ⟨z, he, ?_, ?_⟩ + · have hsum : (y.regs population + y.regs residue % 2) + + MaskTable.population (L - (y.regs counter + 1)) (y.regs residue / 2) = + MaskTable.population L a := by omega + simpa +decide [I, z, levels, one, two, residue, population, counter, digit, + Function.update_apply] using + (show y.mem = x.mem ∧ y.regs levels = L ∧ y.regs one = 1 ∧ y.regs two = 2 ∧ + y.regs counter + 1 ≤ L ∧ y.regs residue / 2 ≤ a ∧ + (y.regs population + y.regs residue % 2) + + MaskTable.population (L - (y.regs counter + 1)) (y.regs residue / 2) = + MaskTable.population L a from + ⟨hm, hL, h1, h2, by omega, (Nat.div_le_self _ _).trans hr, hsum⟩) + · simp only [z, Data.set_regs, Function.update_self] + omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => L - y.regs counter) hb hi + have hj : y.regs counter = L := by + simp only [Cond.holds_lt, hy.2.1] at hc + have := hy.2.2.2.2.1 + omega + have hwhole : ExecWithin w populationCode x y (6 * L + 5) := by + have hs : ExecWithin (B := .standard) w (.op (.copy residue mask)) + x (x.set residue a) 1 := by + simpa only [ha] using (ExecWithin.copy (B := .standard) (w := w) + (x := x) (d := residue) (s := mask)) + have hwhole := hs.seq ((ExecWithin.const (d := population) hzero).seq + ((ExecWithin.const (d := counter) hzero).seq he)) + apply hwhole.mono + simp only [start, Data.set_regs, Function.update_self, Nat.sub_zero] + omega + refine ⟨y, hwhole, ?_, hy.1, ?_⟩ + · have h := hy.2.2.2.2.2.2 + simpa only [hj, Nat.sub_self, population_zero, Nat.add_zero] using h + · intro i hr hp hc hd + apply hwhole.regs_eq + have hframe : ∀ i : Fin 96, i ≠ residue → i ≠ population → i ≠ counter → i ≠ digit → + populationCode.Preserves i := by decide + exact hframe i hr hp hc hd + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPrepare.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPrepare.lean new file mode 100644 index 0000000000..1364489d73 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandPrepare.lean @@ -0,0 +1,182 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandRead + +/-! +# Complete native padded-band preparation + +## Main results + +* `BandPreparation.exec_prepareBand`: initialize and fill every band-array cell by actual + native digit decoding, inverse-table lookup and guarded signed matrix loads. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +-- The static frame checks reduce finite programs over 96 registers and their parameter sets. +set_option maxRecDepth 4096 + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Every prepared-band cell is below the protected original input. -/ +theorem band_input_le_snapshot (L m N : ℕ) : 2 + 7 ^ L ≤ snapshotBase L m N := by + have h := workspace_le_snapshot L m N + unfold Encoding.scratchBase Encoding.outputBase at h + have := Nat.zero_le (10 ^ L) + omega + +private theorem split_preserves_parameters : + ∀ i ∈ parameterRegisters, splitCode.Preserves i := by + simp only [parameterRegisters, headerRegisters, splitCode, splitHead, splitBody, + Code.Preserves, Op.dest, Op.bin] + decide + +private theorem prepare_preserves_parameters (side : Bool) : + ∀ i ∈ parameterRegisters, (prepareBody side).Preserves i := by + cases side <;> simp +decide only [] + +/-- The entire input band is materialized with charged initialization, digit/index work, +padding and loop control. No cell outside the depth header and the band input is changed. -/ +theorem exec_prepareBand {w L m N : ℕ} (side : Bool) (β : ℕ) (x : Data 96) + (hx : BandParameters L m N x) (hβ : β < numBands L m N) (hb : x.regs band = β) + (htable : ∀ t < 2 ^ L, x.mem (inverseBase L m N + t) = MaskTable.lookup L m t) + (matrix : ℕ → ℕ → ℤ) + (hmat : ∀ r < N, ∀ k < D m, + x.mem (matrixAddress L m N side r k) = encode (matrix r k)) + (hfit : Fits w (wordBudget L m N)) : + ∃ y, ExecWithin w (prepareBand side) x y ((13 * L + 39) * 7 ^ L + 5) ∧ + BandParameters L m N y ∧ y.regs band = β ∧ y.mem 1 = L ∧ + (∀ a < 7 ^ L, y.mem (2 + a) = encode (bandValue L m N β side matrix a)) ∧ + (∀ b, ¬(1 ≤ b ∧ b < 2 + 7 ^ L) → y.mem b = x.mem b) := by + have hL : Fits w L := (local_fits hfit).of_le (by + nlinarith only [Nat.zero_le L, Nat.zero_le N, Nat.zero_le (N * D m), + Nat.zero_le (2 ^ L), Nat.zero_le (10 ^ L)]) + let start := ((x.set address 1).write 1 L).set index 0 + have hs : ExecWithin w prepareHead x start 3 := by + ram_eval [prepareHead, start, hx.fixed.levels, hL, + small_fits hfit (by decide : 1 ≤ 256), address, index, levels, Fin.mk.injEq] + have hsp : BandParameters L m N start := parameters_reframe hx (by + intro i hi + apply hs.regs_eq + have hp : ∀ i ∈ parameterRegisters, prepareHead.Preserves i := by + simp only [parameterRegisters, headerRegisters, prepareHead, Code.Preserves, Op.dest] + decide + exact hp i hi) + let I (y : Data 96) : Prop := BandParameters L m N y ∧ y.regs band = β ∧ + y.regs index ≤ 7 ^ L ∧ y.mem 1 = L ∧ + (∀ a < y.regs index, y.mem (2 + a) = encode (bandValue L m N β side matrix a)) ∧ + (∀ b, ¬(1 ≤ b ∧ b < 2 + 7 ^ L) → y.mem b = x.mem b) + have hstart : I start := by + refine ⟨hsp, ?_, ?_, ?_, ?_, ?_⟩ + · simpa [start, band, address, index] using hb + · simp [start] + · simp [start] + · simp [start] + · intro b hbb + have hne : b ≠ 1 := by + have := Nat.one_le_pow L 7 (by decide) + omega + simp [start, hne] + have hbody (y : Data 96) (hy : I y) (hc : (Cond.lt index cells).Holds y) : + ∃ z, ExecWithin w (prepareBody side) y z (13 * L + 37) ∧ I z ∧ + 7 ^ L - z.regs index < 7 ^ L - y.regs index := by + let i := y.regs index + have hi : i < 7 ^ L := by simpa only [Cond.holds_lt, hy.1.cells] using hc + have hsplit : SplitParameters L y := ⟨hy.1.fixed.levels, hy.1.fixed.one, + hy.1.fixed.two, hy.1.fixed.three, hy.1.fixed.four, hy.1.fixed.seven⟩ + obtain ⟨s, hes, _, hmask, houter, hinner, hsmem⟩ := + exec_split L i y hsplit rfl (split_fits hfit hi) + have hsp : BandParameters L m N s := parameters_reframe hy.1 + (fun j hj => hes.regs_eq (split_preserves_parameters j hj)) + have hsband : s.regs band = β := (hes.regs_eq (by + simp only [splitCode, splitHead, splitBody, Code.Preserves, Op.dest, Op.bin] + decide)).trans hy.2.1 + have hsindex : s.regs index = i := hes.regs_eq (by + simp only [splitCode, splitHead, splitBody, Code.Preserves, Op.dest, Op.bin] + decide) + have htable' (t : ℕ) (ht : t < 2 ^ L) : + s.mem (inverseBase L m N + t) = MaskTable.lookup L m t := by + rw [hsmem, hy.2.2.2.2.2 _ (by + have hlow := band_input_le_snapshot L m N + unfold inverseBase + omega)] + exact htable t ht + have hmat' (r : ℕ) (hr : r < N) (k : ℕ) (hk : k < D m) : + s.mem (matrixAddress L m N side r k) = encode (matrix r k) := by + rw [hsmem, hy.2.2.2.2.2 _ (by + have hlow := band_input_le_snapshot L m N + unfold matrixAddress + omega)] + exact hmat r hr k hk + obtain ⟨t, het, hv, htmem⟩ := exec_readCell side ⟨i, hi⟩ β s hsp hβ hsband + hmask houter hinner htable' matrix hmat' hfit + have htindex : t.regs index = i := (het.regs_eq (by + cases side <;> simp only [readCell, Code.Preserves, Op.dest, Op.bin] <;> decide)).trans + hsindex + have htone : t.regs one = 1 := (het.regs_eq (by + cases side <;> simp only [readCell, Code.Preserves, Op.dest, Op.bin] <;> decide)).trans + hsp.fixed.one + have httwo : t.regs two = 2 := (het.regs_eq (by + cases side <;> simp only [readCell, Code.Preserves, Op.dest, Op.bin] <;> decide)).trans + hsp.fixed.two + have hadd : Fits w (i + 2) := + ((memory_fits hfit).of_le (snapshot_le_memoryEnd L m N)).of_le (by + have hlow := band_input_le_snapshot L m N + omega) + have hinc : Fits w (i + 1) := hadd.of_le (by omega) + let z := ((t.set address (i + 2)).write (i + 2) + (encode (bandValue L m N β side matrix i))).set index (i + 1) + have hez : ExecWithin w (ram[ + address ← index + two; *address ← value; index ← index + one ]) t z 3 := by + ram_eval [z, htindex, htone, httwo, hv, hadd, hinc, + address, index, two, one, value, Fin.mk.injEq] + have he : ExecWithin w (prepareBody side) y z (13 * L + 37) := by + simpa only [prepareBody, Nat.add_assoc] using hes.seq (het.seq hez) + have hzindex : z.regs index = i + 1 := by simp [z] + have hzmem : z.mem = Function.update y.mem (i + 2) + (encode (bandValue L m N β side matrix i)) := by + simp only [z, Data.set_mem, Data.write_mem, htmem, hsmem] + refine ⟨z, he, ?_, ?_⟩ + · refine ⟨parameters_reframe hy.1 (fun j hj => + he.regs_eq (prepare_preserves_parameters side j hj)), ?_, ?_, ?_, ?_, ?_⟩ + · exact (he.regs_eq (by + cases side <;> + simp +decide only [])).trans hy.2.1 + · rw [hzindex]; omega + · rw [hzmem, Function.update_of_ne (by omega : 1 ≠ i + 2)] + exact hy.2.2.2.1 + · intro a ha + rw [hzindex] at ha + rw [hzmem] + by_cases hai : a = i + · subst a; rw [Nat.add_comm 2 i, Function.update_self] + · rw [Function.update_of_ne (by omega : 2 + a ≠ i + 2)] + exact hy.2.2.2.2.1 a (by omega) + · intro b hbb + rw [hzmem, Function.update_of_ne (by omega : b ≠ i + 2)] + exact hy.2.2.2.2.2 b hbb + · rw [hzindex]; change 7 ^ L - (i + 1) < 7 ^ L - i; omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => 7 ^ L - y.regs index) hbody hstart + have hindex : y.regs index = 7 ^ L := by + simp only [Cond.holds_lt, hy.1.cells] at hc + have := hy.2.2.1 + omega + refine ⟨y, ?_, hy.1, hy.2.1, hy.2.2.2.1, ?_, hy.2.2.2.2.2⟩ + · apply (hs.seq he).mono + have hzero : start.regs index = 0 := by simp [start] + rw [hzero, Nat.sub_zero] + nlinarith only [Nat.zero_le (7 ^ L)] + · simpa only [hindex] using hy.2.2.2.2.1 + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandProgram.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandProgram.lean new file mode 100644 index 0000000000..127663f3c0 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandProgram.lean @@ -0,0 +1,322 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.Flat +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.ThinProduct +public import Complexity.Computability.WordRAM.Power + +/-! +# Native preparation and encoding of shared matrix bands + +The input is `L :: m :: N :: X ++ Y`, with both matrices in row-major order and their +signed entries encoded by `Signed.encode`. The native length header occupies address zero. +The program snapshots the complete initialized input before using the encoder's low workspace, +computes its subset tables, prepares and encodes each band once, and restores the input. + +The dense inverse table uses zero for an unselected mask and rank plus one for a selected mask. +The forward table lists the first `K0²` masks of cardinality `m` in increasing numeric order. +The shared output catalogue consists of all left bands followed by all right bands, with +`10^L` consecutive cells per band. A tile refers to its two band entries in this catalogue. + +`LeftVar.equiv` and `RightVar.equiv` identify the mathematical alphabets with +`Fin 3 ⊕ (Fin 2 × Fin 2)`, and `Term.equiv` identifies terms with `Option (Fin 3 × Fin 3)`. +The existing native `Encoding.leftDigits`, `Encoding.rightDigits` and `Encoding.termDigits` +also fix the executable flat `Fin 7`/`Fin 10` order: the nine product terms are row-major, +with digit `3*i+j`, and `P0` has digit 9. This preparation retains that numbering unchanged. + +## Main definitions + +* `BandPreparation.program`: fixed native code, including the actual shared encoders. +* `BandPreparation.input`: flat header and matrix input for `State.init`. +* `BandPreparation.catalogueAddress`: the location of a band's encoding leaf. + +## References + +Selective adaptation of the band preparation in [AlmanWilliams2026, Section 2.4.1] and +Anthropic's formal-math revision `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`, +[band-array preparation](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Programs/Sec2/Theorem5/Encode/BandArray.lean) +and [all-band encoding](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Programs/Sec2/Theorem5/Encode/AllBands.lean). +The mask scans, dense inverse table, native control and memory layout here are new; +the source's Light machine and execution semantics are not used. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Original initialized input prefix, including the native length header and three parameters. -/ +def inputEnd (m N : ℕ) : ℕ := 4 + 2 * (N * D m) + +/-- The encoder's low workspace and the original input are both below the snapshot. -/ +def snapshotBase (L m N : ℕ) : ℕ := inputEnd m N + Encoding.scratchBase L + 10 ^ L + +/-- Dense mask-to-rank inverse table. -/ +def inverseBase (L m N : ℕ) : ℕ := snapshotBase L m N + inputEnd m N + +/-- Forward selected-mask table; `2^L` cells are reserved. -/ +def tableBase (L m N : ℕ) : ℕ := inverseBase L m N + 2 ^ L + +/-- Shared encoding catalogue, after the snapshot and both tables. -/ +def catalogueBase (L m N : ℕ) : ℕ := tableBase L m N + 2 ^ L + +/-- Every left band, then every right band, occupies one block of `10^L` words. -/ +def catalogueAddress (L m N : ℕ) (rightSide : Bool) (β leaf : ℕ) : ℕ := + catalogueBase L m N + ((if rightSide then numBands L m N else 0) + β) * 10 ^ L + leaf + +/-- One past all allocated memory. -/ +def memoryEnd (L m N : ℕ) : ℕ := catalogueBase L m N + 2 * numBands L m N * 10 ^ L + +/-- Explicit input-parameter bound used to derive every unsigned arithmetic/address fit. +The signed magnitude condition is separate. This is an ordinary numerical inequality, +without quantification over executions or intermediate states. -/ +def wordBudget (L m N : ℕ) : ℕ := + 256 * (L + N + N * D m + 2 ^ L + 10 ^ L + 2 * numBands L m N * 10 ^ L + 1) + +/-- Row-major signed matrix input; no table, power or encoded band is supplied by the caller. -/ +def input (L m N : ℕ) (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) : List ℕ := + [L, m, N] ++ + (List.finRange N).flatMap (fun i => List.ofFn fun k => Signed.encode (X i k)) ++ + (List.finRange (D m)).flatMap (fun k => List.ofFn fun j => Signed.encode (Y k j)) + +/-- Low registers are reserved for the existing encoder and copy routines. Persistent +parameters and preparation counters occupy registers 32 through 84 of a fixed bank of 96. -/ +abbrev levels : Fin 96 := 32 +@[inherit_doc levels] abbrev copySrc : Fin 96 := 0 +@[inherit_doc levels] abbrev copyDst : Fin 96 := 1 +@[inherit_doc levels] abbrev copyLen : Fin 96 := 2 +@[inherit_doc levels] abbrev inner : Fin 96 := 33 +@[inherit_doc levels] abbrev rows : Fin 96 := 34 +@[inherit_doc levels] abbrev columns : Fin 96 := 35 +@[inherit_doc levels] abbrev outerSize : Fin 96 := 36 +@[inherit_doc levels] abbrev masks : Fin 96 := 37 +@[inherit_doc levels] abbrev cells : Fin 96 := 38 +@[inherit_doc levels] abbrev leaves : Fin 96 := 39 +@[inherit_doc levels] abbrev prefixSize : Fin 96 := 40 +@[inherit_doc levels] abbrev snapshot : Fin 96 := 41 +@[inherit_doc levels] abbrev inverse : Fin 96 := 42 +@[inherit_doc levels] abbrev table : Fin 96 := 43 +@[inherit_doc levels] abbrev catalogue : Fin 96 := 44 +@[inherit_doc levels] abbrev bands : Fin 96 := 45 +@[inherit_doc levels] abbrev count : Fin 96 := 46 +@[inherit_doc levels] abbrev root : Fin 96 := 47 +@[inherit_doc levels] abbrev square : Fin 96 := 48 +@[inherit_doc levels] abbrev temporary : Fin 96 := 49 +@[inherit_doc levels] abbrev one : Fin 96 := 50 +@[inherit_doc levels] abbrev two : Fin 96 := 51 +@[inherit_doc levels] abbrev three : Fin 96 := 52 +@[inherit_doc levels] abbrev four : Fin 96 := 53 +@[inherit_doc levels] abbrev seven : Fin 96 := 54 +@[inherit_doc levels] abbrev ten : Fin 96 := 55 +@[inherit_doc levels] abbrev counter : Fin 96 := 56 +@[inherit_doc levels] abbrev mask : Fin 96 := 57 +@[inherit_doc levels] abbrev population : Fin 96 := 58 +@[inherit_doc levels] abbrev residue : Fin 96 := 59 +@[inherit_doc levels] abbrev digit : Fin 96 := 60 +@[inherit_doc levels] abbrev limit : Fin 96 := 61 +@[inherit_doc levels] abbrev rank : Fin 96 := 62 +@[inherit_doc levels] abbrev address : Fin 96 := 63 +@[inherit_doc levels] abbrev band : Fin 96 := 64 +@[inherit_doc levels] abbrev index : Fin 96 := 66 +@[inherit_doc levels] abbrev maskAcc : Fin 96 := 67 +@[inherit_doc levels] abbrev outerAcc : Fin 96 := 68 +@[inherit_doc levels] abbrev innerAcc : Fin 96 := 69 +@[inherit_doc levels] abbrev maskPower : Fin 96 := 70 +@[inherit_doc levels] abbrev outerPower : Fin 96 := 71 +@[inherit_doc levels] abbrev innerPower : Fin 96 := 72 +@[inherit_doc levels] abbrev level : Fin 96 := 73 +@[inherit_doc levels] abbrev block : Fin 96 := 75 +@[inherit_doc levels] abbrev row : Fin 96 := 77 +@[inherit_doc levels] abbrev source : Fin 96 := 79 +@[inherit_doc levels] abbrev value : Fin 96 := 80 +@[inherit_doc levels] abbrev zero : Fin 96 := 84 + +/-- Register embedding for the unchanged shared encoders. -/ +def encoderRegisters (i : Fin 32) : Fin 96 := ⟨i.val, by omega⟩ + +/-- Register embedding for the unchanged seven-register copy program. -/ +def copyRegisters (i : Fin 7) : Fin 96 := ⟨i.val, by omega⟩ + +/-- Embed the five-register native power routine in the caller's chosen scalar registers. -/ +def powerRegisters (v b n : Fin 96) : Fin 5 → Fin 96 := ![v, counter, b, n, one] + +/-- The actual native power routine, with register renaming only. -/ +def powerCode (v b n : Fin 96) : Code .standard 96 := + Power.program.rename (powerRegisters v b n) + +/-- Popcount is computed by an explicit fixed-depth loop of divisions and remainders. -/ +def populationCode : Code .standard 96 := ram[ + residue ← mask; population ← #0; counter ← #0; + while counter < levels do [ + digit ← residue % two; population ← population + digit; + residue ← residue / two; counter ← counter + one ] ] + +/-- Load the three actual header parameters and the small constants. -/ +def headerLoads : Code .standard 96 := ram[ + address ← #1; levels ← *address; address ← #2; inner ← *address; + address ← #3; rows ← *address; + zero ← #0; one ← #1; two ← #2; three ← #3; four ← #4; seven ← #7; ten ← #10 ] + +/-- Compute all five powers with actual renamed native routines. -/ +def headerPowers : Code .standard 96 := + .seq (powerCode columns four inner) (.seq (ram[ temporary ← levels - inner ]) + (.seq (powerCode outerSize three temporary) (.seq (powerCode masks two levels) + (.seq (powerCode cells seven levels) (powerCode leaves ten levels))))) + +/-- Compute concrete snapshot, table and catalogue addresses. -/ +def headerLayout : Code .standard 96 := ram[ + prefixSize ← rows * columns; prefixSize ← prefixSize * two; + prefixSize ← prefixSize + four; + snapshot ← leaves * two; snapshot ← snapshot + cells; snapshot ← snapshot + two; + snapshot ← snapshot + prefixSize; + inverse ← snapshot + prefixSize; table ← inverse + masks; catalogue ← table + masks ] + +/-- Compute the three header parameters, powers and addresses with charged native code. -/ +def header : Code .standard 96 := .seq headerLoads (.seq headerPowers headerLayout) + +/-- Save every original input word, including the actual length header. -/ +def saveInput : Code .standard 96 := ram[ + copySrc ← #0; copyDst ← snapshot; copyLen ← prefixSize; + ⟪ArrayCopy.program.rename copyRegisters⟫ ] + +/-- Visit one mask and update the number of accepted masks. -/ +def countBody : Code .standard 96 := + .seq populationCode (ram[ + if population = inner then [ count ← count + one ] else [ count ← count ]; + mask ← mask + one ]) + +/-- Count every size-`m` mask. -/ +def countMasks : Code .standard 96 := + .seq (ram[ mask ← #0; count ← #0 ]) (.loop (.lt mask masks) countBody) + +/-- One explicitly squared trial in the integer square-root search. -/ +def rootBody : Code .standard 96 := ram[ + root ← root + one; temporary ← root + one; square ← temporary * temporary ] + +/-- Compute the integer square root by squared comparisons. No square-root instruction +is added to the native basis. -/ +def rootCode : Code .standard 96 := + .seq (ram[ + root ← #0; square ← #1; limit ← count + one + ]) (.seq (.loop (.lt square limit) rootBody) (ram[ square ← root * root ])) + +/-- Compute the ceiling quotient giving the number of bands. -/ +def bandCount : Code .standard 96 := ram[ + temporary ← root * outerSize; + bands ← rows + temporary; bands ← bands - one; bands ← bands / temporary ] + +/-- Complete scalar band parameters from native mask counting and root computation. -/ +def countCode : Code .standard 96 := .seq countMasks (.seq rootCode bandCount) + +/-- Store the two table entries belonging to the current mask and advance both counters. -/ +def tableWrite : Code .standard 96 := ram[ + address ← inverse + mask; *address ← zero; + if population = inner then [ + if rank < square then [ + temporary ← rank + one; *address ← temporary; + address ← table + rank; *address ← mask ] else [ rank ← rank ]; + rank ← rank + one ] else [ rank ← rank ]; + mask ← mask + one ] + +/-- One complete table iteration, including the actual population scan. -/ +def tableBody : Code .standard 96 := .seq populationCode tableWrite + +/-- Materialize both tables, explicitly initializing every inverse-table cell. -/ +def tablesCode : Code .standard 96 := + .seq (ram[ mask ← #0; rank ← #0 ]) (.loop (.lt mask masks) tableBody) + +/-- Initialize the three digit accumulators and their positional weights. -/ +def splitHead : Code .standard 96 := ram[ + residue ← index; maskAcc ← #0; outerAcc ← #0; innerAcc ← #0; + maskPower ← #1; outerPower ← #1; innerPower ← #1; level ← #0 ] + +/-- One base-seven digit, with only ordinary native arithmetic and control. -/ +def splitBody : Code .standard 96 := ram[ + digit ← residue % seven; residue ← residue / seven; + if digit < three then [ + temporary ← digit * outerPower; outerAcc ← outerAcc + temporary; + outerPower ← outerPower * three ] else [ + digit ← digit - three; temporary ← digit * innerPower; + innerAcc ← innerAcc + temporary; innerPower ← innerPower * four; + maskAcc ← maskAcc + maskPower ]; + maskPower ← maskPower * two; level ← level + one ] + +/-- Decode a base-seven string into its inner mask and compressed base-three/base-four +offsets, in increasing level order. Every division, remainder, product and increment is native. -/ +def splitCode : Code .standard 96 := .seq splitHead (.loop (.lt level levels) splitBody) + +/-- Read a real matrix cell only after checking the padded row/column against `N`. +In particular, the potentially oversized row-major product is inside the guard. -/ +def readCell (rightSide : Bool) : Code .standard 96 := + .seq (ram[ + address ← inverse + maskAcc; rank ← *address; value ← #0 ]) + (.branch (.lt zero rank) + (.seq (ram[ rank ← rank - one ]) + (.seq (if rightSide then ram[ block ← rank % root ] else ram[ block ← rank / root ]) + (ram[ + temporary ← band * root; block ← temporary + block; + row ← block * outerSize; row ← row + outerAcc; + if row < rows then [ + ⟪if rightSide then ram[ + address ← innerAcc * rows; address ← address + row; + temporary ← rows * columns; address ← address + temporary ] + else ram[ address ← row * columns; address ← address + innerAcc ]⟫; + address ← address + four; address ← snapshot + address; value ← *address + ] else [ value ← value ] ]))) + (ram[ value ← value ])) + +/-- Actual header and counter initialization for each band. -/ +def prepareHead : Code .standard 96 := ram[ + address ← #1; *address ← levels; index ← #0 ] + +/-- Decode, guard, load and store one band-array entry. -/ +def prepareBody (rightSide : Bool) : Code .standard 96 := ram[ + ⟪splitCode⟫; ⟪readCell rightSide⟫; + address ← index + two; *address ← value; index ← index + one ] + +/-- Prepare every entry of one low-workspace band array, including all zero padding. -/ +def prepareBand (rightSide : Bool) : Code .standard 96 := + .seq prepareHead (.loop (.lt index cells) (prepareBody rightSide)) + +/-- Set up the actual encoder-output copy, including its catalogue address arithmetic. -/ +def collectHead (rightSide : Bool) : Code .standard 96 := ram[ + copySrc ← cells + two; + ⟪if rightSide then ram[ temporary ← bands + band ] else ram[ temporary ← band ]⟫; + temporary ← temporary * leaves; temporary ← catalogue + temporary; + copyDst ← temporary; copyLen ← leaves ] + +/-- Copy the actual encoder output to this band's shared catalogue entry. -/ +def collectBand (rightSide : Bool) : Code .standard 96 := + .seq (collectHead rightSide) (ArrayCopy.program.rename copyRegisters) + +/-- Prepare, encode and collect one shared band, then advance its counter. -/ +def bandBody (rightSide : Bool) : Code .standard 96 := + .seq (.seq (prepareBand rightSide) + (.seq ((Encoding.program rightSide).rename encoderRegisters) (collectBand rightSide))) + (ram[ band ← band + one ]) + +/-- Each row/column band calls the actual unchanged finite native encoder exactly once. -/ +def encodeBands (rightSide : Bool) : Code .standard 96 := + .seq (ram[ band ← #0 ]) (.loop (.lt band bands) (bandBody rightSide)) + +/-- Restore the complete original initialized input after all low-workspace calls. -/ +def restoreInput : Code .standard 96 := ram[ + copySrc ← snapshot; copyDst ← #0; copyLen ← prefixSize; + ⟪ArrayCopy.program.rename copyRegisters⟫ ] + +/-- One fixed finite native program for the entire band-preparation component. Halt returns +`catalogue`, the base address of the shared encodings; it is appended by `Code.program`. -/ +def program : Code .standard 96 := + .seq header (.seq saveInput (.seq countCode (.seq tablesCode + (.seq (encodeBands false) (.seq (encodeBands true) restoreInput))))) + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRead.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRead.lean new file mode 100644 index 0000000000..cc88b44cc8 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRead.lean @@ -0,0 +1,196 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandParameters +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandMemory +public import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandSplit + +/-! +# Guarded native matrix reads for band arrays + +## Main definitions + +* `BandPreparation.matrixAddress`: real row-major matrix addresses in the protected input. + +## Main results + +* `BandPreparation.exec_readCell`: the actual guarded read returns exactly one encoded + band value; padded rows are rejected before row-major multiplication. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open BandDigits + +/-- Offset in the concatenated original matrix inputs. The right matrix has `D` rows and +`N` columns, while its band specification labels a column by `r`. -/ +def matrixOffset (m N : ℕ) (rightSide : Bool) (r k : ℕ) : ℕ := + if rightSide then k * N + r + N * D m else r * D m + k + +/-- Address of an original signed input word in its protected snapshot. -/ +def matrixAddress (L m N : ℕ) (rightSide : Bool) (r k : ℕ) : ℕ := + snapshotBase L m N + (matrixOffset m N rightSide r k + 4) + +private theorem matrix_indices_lt {m N r k : ℕ} (hr : r < N) (hk : k < D m) : + r * D m + k < N * D m ∧ k * N + r < N * D m := by + have hleft := Nat.mul_le_mul_right (D m) (Nat.succ_le_of_lt hr) + have hright := Nat.mul_le_mul_right N (Nat.succ_le_of_lt hk) + rw [Nat.succ_mul] at hleft hright + rw [Nat.mul_comm (D m) N] at hright + constructor <;> omega + +/-- A real matrix cell lies strictly inside its protected snapshot. This uses the +real row/column guard, before constructing its flat address. -/ +theorem matrixAddress_lt_inverse {L m N r k : ℕ} (side : Bool) + (hr : r < N) (hk : k < D m) : + matrixAddress L m N side r k < inverseBase L m N := by + obtain ⟨hl, hr'⟩ := matrix_indices_lt hr hk + unfold matrixAddress matrixOffset inverseBase inputEnd + cases side <;> simp only [Bool.false_eq_true, ↓reduceIte] <;> omega + +/-- The complete guarded load uses at most twenty-four actual target instructions. +All intermediate products are bounded using the original parameters and the real `r + simp only [Bool.false_eq_true, ↓reduceIte] at hbfit hblockfit hrowfit hrfit hr <;> + simp only [Bool.false_eq_true, ↓reduceIte] at hleftMul hleftAdd hrightAdd + <;> simp only [Bool.false_eq_true, ↓reduceIte] at hrightOff hleftFour hrightFour + <;> simp only [Bool.false_eq_true, ↓reduceIte] at hleftAddress hrightAddress hcell <;> + ram_eval [readCell, y, v, u, b, r, hx.inverse, hx.snapshot, hx.root, hx.outerSize, + hx.columns, hx.fixed.zero, hx.fixed.one, hx.fixed.four, hx.fixed.rows, + hm, ho, hi, hb, hlookup, haddr, hqpos, hr, + small_fits hfit (by decide : 0 ≤ 256), hqfit, hbfit, hβfit, hblockfit, + hrowfit, hrfit, hnd, hleftMul, hleftAdd, hrightMul, hrightAdd, + hrightOff, hleftFour, hrightFour, hleftAddress, hrightAddress, + hcell, matrixAddress, matrixOffset, + inverse, snapshot, root, outerSize, columns, zero, one, four, rows, + maskAcc, outerAcc, innerAcc, band, address, rank, value, block, temporary, + row, Fin.mk.injEq] + · change y.regs value = encode (if q = 0 then 0 else + if r < N then matrix r p.inner else 0) + rw [ite_eq_right hq, ite_eq_left hr] + cases side <;> simp [y] + · cases side <;> rfl + · let y := v.set value 0 + dsimp only [b, r] at hbfit hblockfit hrowfit hrfit hr + refine ⟨y, ?_, ?_, rfl⟩ + · cases side <;> + simp only [Bool.false_eq_true, ↓reduceIte] at hbfit hblockfit hrowfit hrfit hr <;> + ram_eval [readCell, y, v, u, b, r, hx.inverse, hx.root, hx.outerSize, + hx.fixed.zero, hx.fixed.one, hx.fixed.rows, hm, ho, hi, hb, + hlookup, haddr, hqpos, hr, small_fits hfit (by decide : 0 ≤ 256), + hqfit, hbfit, hβfit, hblockfit, hrowfit, hrfit, + inverse, root, outerSize, zero, one, rows, maskAcc, outerAcc, band, + address, rank, value, block, temporary, row, Fin.mk.injEq] + · change y.regs value = encode (if q = 0 then 0 else + if r < N then matrix r p.inner else 0) + rw [ite_eq_right hq, ite_eq_right hr] + simp only [y, Data.set_regs, Function.update_self] + rfl + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRoot.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRoot.lean new file mode 100644 index 0000000000..48ca528194 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandRoot.lean @@ -0,0 +1,101 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandProgram + +/-! +# Native integer square root for band parameters + +## Main results + +* `BandPreparation.exec_root`: successive squared comparisons compute the ordinary natural + square root within `5 * sqrt k + 6` instructions, with every trial square proved to fit. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM + +/-- The root loop uses only increment, multiplication and unsigned comparison. The bound on +trial squares is derived from the original count, including the first failing trial. -/ +theorem exec_root {w : ℕ} (k : ℕ) (x : Data 96) + (hk : x.regs count = k) (ho : x.regs one = 1) (hfit : Fits w (4 * k + 4)) : + ∃ y, ExecWithin w rootCode x y (5 * Nat.sqrt k + 6) ∧ + y.regs root = Nat.sqrt k ∧ y.regs square = Nat.sqrt k * Nat.sqrt k ∧ + y.mem = x.mem ∧ + (∀ i, i ≠ root → i ≠ square → i ≠ temporary → i ≠ limit → y.regs i = x.regs i) := by + have h0 : Fits w 0 := hfit.of_le (by omega) + have h1 : Fits w 1 := hfit.of_le (by omega) + have hnext : Fits w (k + 1) := hfit.of_le (by omega) + let start := ((x.set root 0).set square 1).set limit (k + 1) + have hs : ExecWithin w (ram[ root ← #0; square ← #1; limit ← count + one ]) + x start 3 := by + ram_eval [start, hk, ho, h0, h1, hnext, root, square, limit, count, one, Fin.mk.injEq] + let I (y : Data 96) : Prop := + y.mem = x.mem ∧ y.regs one = 1 ∧ y.regs limit = k + 1 ∧ + y.regs root ≤ Nat.sqrt k ∧ + y.regs square = (y.regs root + 1) * (y.regs root + 1) + have hi : I start := by simp [I, start, ho, root, square, limit, one] + have hb (y : Data 96) (hy : I y) (hc : (Cond.lt square limit).Holds y) : + ∃ z, ExecWithin w rootBody y z 3 ∧ I z ∧ + Nat.sqrt k - z.regs root < Nat.sqrt k - y.regs root := by + obtain ⟨hm, h1, hlimit, hr, hsq⟩ := hy + have htrial : (y.regs root + 1) * (y.regs root + 1) ≤ k := by + simpa only [Cond.holds_lt, hlimit, hsq, Nat.lt_add_one_iff] using hc + have hrnext : y.regs root + 1 ≤ Nat.sqrt k := Nat.le_sqrt.mpr htrial + have hlarge : (y.regs root + 2) * (y.regs root + 2) ≤ 4 * k + 4 := by + nlinarith only [htrial] + have hrfit : Fits w (y.regs root + 1) := + hfit.of_le (by have := Nat.sqrt_le_self k; omega) + have htfit : Fits w (y.regs root + 2) := + hfit.of_le (by have := Nat.sqrt_le_self k; omega) + have hsfit : Fits w ((y.regs root + 2) * (y.regs root + 2)) := hfit.of_le hlarge + let z := ((y.set root (y.regs root + 1)).set temporary (y.regs root + 2)).set + square ((y.regs root + 2) * (y.regs root + 2)) + have he : ExecWithin w rootBody y z 3 := by + ram_eval [rootBody, z, h1, hrfit, htfit, hsfit, + root, temporary, square, one, Fin.mk.injEq, Nat.add_assoc] + refine ⟨z, he, ?_, ?_⟩ + · simpa +decide [I, z, root, temporary, square, one, limit, Function.update_apply, + Nat.add_assoc] using + (show y.mem = x.mem ∧ y.regs one = 1 ∧ y.regs limit = k + 1 ∧ + y.regs root + 1 ≤ Nat.sqrt k ∧ + (y.regs root + 2) * (y.regs root + 2) = + (y.regs root + 1 + 1) * (y.regs root + 1 + 1) from + ⟨hm, h1, hlimit, hrnext, by ring⟩) + · have hz : z.regs root = y.regs root + 1 := by simp [z, root, square, temporary] + rw [hz] + omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => Nat.sqrt k - y.regs root) hb hi + have hr : y.regs root = Nat.sqrt k := by + have hstop : k < (y.regs root + 1) * (y.regs root + 1) := by + simp only [Cond.holds_lt, hy.2.2.1, hy.2.2.2.2] at hc + omega + exact Nat.eq_sqrt.mpr ⟨(Nat.le_sqrt.mp hy.2.2.2.1), hstop⟩ + let z := y.set square (Nat.sqrt k * Nat.sqrt k) + have ht : ExecWithin w (ram[ square ← root * root ]) y z 1 := by + ram_eval [z, hr, hfit.of_le (Nat.sqrt_le k |>.trans (by omega)), square, root] + have hwhole : ExecWithin w rootCode x z (5 * Nat.sqrt k + 6) := by + apply (hs.seq (he.seq ht)).mono + have hsroot : start.regs root = 0 := by simp [start, root, square, limit] + rw [hsroot, Nat.sub_zero] + omega + refine ⟨z, hwhole, ?_, ?_, hy.1, ?_⟩ + · simpa [z, square, root] using hr + · simp [z] + · intro i hr hs ht hl + apply hwhole.regs_eq + have hf : ∀ i : Fin 96, i ≠ root → i ≠ square → i ≠ temporary → i ≠ limit → + rootCode.Preserves i := by decide + exact hf i hr hs ht hl + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandSplit.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandSplit.lean new file mode 100644 index 0000000000..30463b6899 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandSplit.lean @@ -0,0 +1,247 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandMeaning + +/-! +# Charged native splitting of tensor-string digits + +## Main results + +* `BandPreparation.exec_split`: the actual base-seven digit loop computes the binary mask + and both compressed offsets with every native arithmetic operation charged. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open BandDigits + +/-- Fixed scalar inputs to a digit loop, without execution or fit assumptions. -/ +structure SplitParameters (L : ℕ) (x : Data 96) : Prop where + /-- Original depth. -/ + levels : x.regs levels = L + /-- Constant one. -/ + one : x.regs one = 1 + /-- Constant two. -/ + two : x.regs two = 2 + /-- Constant three. -/ + three : x.regs three = 3 + /-- Constant four. -/ + four : x.regs four = 4 + /-- Constant seven. -/ + seven : x.regs seven = 7 + +private def splitParameterRegisters : Finset (Fin 96) := {levels, one, two, three, four, seven} + +private theorem SplitParameters.of_eq {L : ℕ} {x y : Data 96} (hx : SplitParameters L x) + (h : ∀ i ∈ splitParameterRegisters, y.regs i = x.regs i) : SplitParameters L y := + ⟨(h _ (by decide)).trans hx.levels, (h _ (by decide)).trans hx.one, + (h _ (by decide)).trans hx.two, (h _ (by decide)).trans hx.three, + (h _ (by decide)).trans hx.four, (h _ (by decide)).trans hx.seven⟩ + +/-- Remaining-digit equations and elementary counter bounds during the native scan. -/ +structure SplitProgress (L a : ℕ) (x : Data 96) : Prop where + /-- Parameters are retained. -/ + fixed : SplitParameters L x + /-- Number of consumed levels. -/ + level_le : x.regs level ≤ L + /-- Division only decreases the residual index. -/ + residue_le : x.regs residue ≤ a + /-- Binary place weight. -/ + mask_power : x.regs maskPower = 2 ^ x.regs level + /-- At most one outer place is consumed per level. -/ + outer_power : x.regs outerPower ≤ 3 ^ x.regs level + /-- At most one inner place is consumed per level. -/ + inner_power : x.regs innerPower ≤ 4 ^ x.regs level + /-- Partial binary value is below its next place. -/ + mask_acc : x.regs maskAcc < x.regs maskPower + /-- Partial base-three value is below its next place. -/ + outer_acc : x.regs outerAcc < x.regs outerPower + /-- Partial base-four value is below its next place. -/ + inner_acc : x.regs innerAcc < x.regs innerPower + /-- The binary accumulator and remaining digits give the complete mask. -/ + mask_value : x.regs maskAcc + x.regs maskPower * + (parts (L - x.regs level) (x.regs residue)).mask = (parts L a).mask + /-- The outer accumulator and remaining digits give the complete outer offset. -/ + outer_value : x.regs outerAcc + x.regs outerPower * + (parts (L - x.regs level) (x.regs residue)).outer = (parts L a).outer + /-- The inner accumulator and remaining digits give the complete inner offset. -/ + inner_value : x.regs innerAcc + x.regs innerPower * + (parts (L - x.regs level) (x.regs residue)).inner = (parts L a).inner + +/-- One native digit iteration, using only bounds derived from the original index/depth. -/ +theorem exec_splitBody {w : ℕ} (L a : ℕ) (x : Data 96) (hx : SplitProgress L a x) + (hc : (Cond.lt level levels).Holds x) (hfit : Fits w (a + L + 16 * 10 ^ L + 16)) : + ∃ y, ExecWithin w splitBody x y 11 ∧ SplitProgress L a y ∧ y.mem = x.mem ∧ + L - y.regs level < L - x.regs level := by + have hj : x.regs level < L := by simpa only [Cond.holds_lt, hx.fixed.levels] using hc + have hd : x.regs residue % 7 < 7 := Nat.mod_lt _ (by decide) + have hpow (b : ℕ) (hb : b ≤ 10) : b ^ x.regs level ≤ 10 ^ L := + (Nat.pow_le_pow_left hb _).trans (Nat.pow_le_pow_right (by decide) hx.level_le) + have hmp : x.regs maskPower ≤ 10 ^ L := by rw [hx.mask_power]; exact hpow 2 (by decide) + have hop : x.regs outerPower ≤ 10 ^ L := hx.outer_power.trans (hpow 3 (by decide)) + have hip : x.regs innerPower ≤ 10 ^ L := hx.inner_power.trans (hpow 4 (by decide)) + have hdFit : Fits w (x.regs residue % 7) := hfit.of_le (by omega) + have hrFit : Fits w (x.regs residue / 7) := + hfit.of_le ((Nat.div_le_self _ _).trans (by have := hx.residue_le; omega)) + have hmFit : Fits w (x.regs maskPower * 2) := hfit.of_le (by omega) + have hjFit : Fits w (x.regs level + 1) := hfit.of_le (by omega) + have hrest : L - x.regs level = L - (x.regs level + 1) + 1 := by omega + by_cases houter : x.regs residue % 7 < 3 + · have hmul : Fits w (x.regs residue % 7 * x.regs outerPower) := + hfit.of_le (by nlinarith) + have hadd : Fits w (x.regs outerAcc + x.regs residue % 7 * x.regs outerPower) := + hfit.of_le (by have := hx.outer_acc; nlinarith) + have hpower : Fits w (x.regs outerPower * 3) := hfit.of_le (by omega) + let y := (((((((x.set digit (x.regs residue % 7)).set residue (x.regs residue / 7)).set + temporary (x.regs residue % 7 * x.regs outerPower)).set outerAcc + (x.regs outerAcc + x.regs residue % 7 * x.regs outerPower)).set outerPower + (x.regs outerPower * 3)).set maskPower (x.regs maskPower * 2)).set + level (x.regs level + 1)) + have he : ExecWithin w splitBody x y 11 := by + ram_eval [splitBody, y, hx.fixed.one, hx.fixed.two, hx.fixed.three, hx.fixed.seven, + hdFit, hrFit, hmFit, hjFit, hmul, hadd, hpower, houter, + digit, residue, temporary, outerAcc, outerPower, maskPower, level, one, two, + three, seven, Fin.mk.injEq] + have hfixed : SplitParameters L y := hx.fixed.of_eq (by + intro i hi + apply he.regs_eq + have hf : ∀ i ∈ splitParameterRegisters, splitBody.Preserves i := by + simp only [splitBody, Code.Preserves, Op.dest, Op.bin] + decide + exact hf i hi) + have hmval := hx.mask_value + have hoval := hx.outer_value + have hival := hx.inner_value + rw [hrest, parts, ite_eq_left houter] at hmval hoval hival + refine ⟨y, he, ?_, rfl, ?_⟩ + · constructor + · exact hfixed + all_goals simp +decide only [y, Data.set_regs, Function.update_apply, ↓reduceIte] + · omega + · exact (Nat.div_le_self _ _).trans hx.residue_le + · simpa only [pow_succ] using congrArg (fun n => n * 2) hx.mask_power + · rw [pow_succ] + exact Nat.mul_le_mul_right 3 hx.outer_power + · exact hx.inner_power.trans (Nat.pow_le_pow_right (by decide) (by omega)) + · have := hx.mask_acc; omega + · have := hx.outer_acc; nlinarith only [this, houter] + · exact hx.inner_acc + · nlinarith only [hmval] + · nlinarith only [hoval] + · exact hival + · have hyl : y.regs level = x.regs level + 1 := by simp [y] + rw [hyl] + omega + · have hsub : Fits w (x.regs residue % 7 - 3) := hfit.of_le (by omega) + have hsuble : x.regs residue % 7 - 3 < 4 := by omega + have hmul : Fits w ((x.regs residue % 7 - 3) * x.regs innerPower) := + hfit.of_le (by nlinarith only [hip, hsuble]) + have hadd : Fits w (x.regs innerAcc + (x.regs residue % 7 - 3) * x.regs innerPower) := + hfit.of_le (by nlinarith only [hx.inner_acc, hip, hsuble]) + have hpower : Fits w (x.regs innerPower * 4) := hfit.of_le (by omega) + have hmask : Fits w (x.regs maskAcc + x.regs maskPower) := + hfit.of_le (by have := hx.mask_acc; omega) + let y := (((((((((x.set digit (x.regs residue % 7)).set residue (x.regs residue / 7)).set + digit (x.regs residue % 7 - 3)).set temporary + ((x.regs residue % 7 - 3) * x.regs innerPower)).set innerAcc + (x.regs innerAcc + (x.regs residue % 7 - 3) * x.regs innerPower)).set innerPower + (x.regs innerPower * 4)).set maskAcc (x.regs maskAcc + x.regs maskPower)).set + maskPower (x.regs maskPower * 2)).set level (x.regs level + 1)) + have he : ExecWithin w splitBody x y 11 := by + ram_eval [splitBody, y, hx.fixed.one, hx.fixed.two, hx.fixed.three, hx.fixed.four, + hx.fixed.seven, hdFit, hrFit, hmFit, hjFit, hsub, hmul, hadd, hpower, hmask, houter, + digit, residue, temporary, innerAcc, innerPower, maskAcc, maskPower, level, + one, two, three, four, seven, Fin.mk.injEq] + have hfixed : SplitParameters L y := hx.fixed.of_eq (by + intro i hi + apply he.regs_eq + have hf : ∀ i ∈ splitParameterRegisters, splitBody.Preserves i := by + simp only [splitBody, Code.Preserves, Op.dest, Op.bin] + decide + exact hf i hi) + have hmval := hx.mask_value + have hoval := hx.outer_value + have hival := hx.inner_value + rw [hrest, parts, ite_eq_right houter] at hmval hoval hival + refine ⟨y, he, ?_, rfl, ?_⟩ + · constructor + · exact hfixed + all_goals simp +decide only [y, Data.set_regs, Function.update_apply, ↓reduceIte] + · omega + · exact (Nat.div_le_self _ _).trans hx.residue_le + · simpa only [pow_succ] using congrArg (fun n => n * 2) hx.mask_power + · exact hx.outer_power.trans (Nat.pow_le_pow_right (by decide) (by omega)) + · rw [pow_succ] + exact Nat.mul_le_mul_right 4 hx.inner_power + · have := hx.mask_acc; omega + · exact hx.outer_acc + · nlinarith only [hx.inner_acc, hsuble] + · nlinarith only [hmval] + · exact hoval + · nlinarith only [hival] + · have hyl : y.regs level = x.regs level + 1 := by simp [y] + rw [hyl] + omega + +/-- All digits are split by the fixed native loop. Scratch need not start at zero: every +accumulator and place weight is explicitly initialized and its initialization is charged. -/ +theorem exec_split {w : ℕ} (L a : ℕ) (x : Data 96) (hx : SplitParameters L x) + (ha : x.regs index = a) (hfit : Fits w (a + L + 16 * 10 ^ L + 16)) : + ∃ y, ExecWithin w splitCode x y (13 * L + 10) ∧ SplitParameters L y ∧ + y.regs maskAcc = (parts L a).mask ∧ y.regs outerAcc = (parts L a).outer ∧ + y.regs innerAcc = (parts L a).inner ∧ y.mem = x.mem := by + have h0 : Fits w 0 := hfit.of_le (Nat.zero_le _) + have h1 : Fits w 1 := hfit.of_le (by omega) + let start := ((((((((x.set residue a).set maskAcc 0).set outerAcc 0).set innerAcc 0).set + maskPower 1).set outerPower 1).set innerPower 1).set level 0) + have hs : ExecWithin w splitHead x start 8 := by + ram_eval [splitHead, start, ha, h0, h1, residue, index, maskAcc, outerAcc, innerAcc, + maskPower, outerPower, innerPower, level, Fin.mk.injEq] + have hp : SplitParameters L start := hx.of_eq (by + intro i hi + apply hs.regs_eq + have hf : ∀ i ∈ splitParameterRegisters, splitHead.Preserves i := by + simp only [splitHead, Code.Preserves, Op.dest] + decide + exact hf i hi) + have hi : SplitProgress L a start := by + constructor + · exact hp + all_goals simp [start, residue, maskAcc, outerAcc, innerAcc, + maskPower, outerPower, innerPower, level] + let I (y : Data 96) : Prop := SplitProgress L a y ∧ y.mem = x.mem + have hb (y : Data 96) (hy : I y) (hc : (Cond.lt level levels).Holds y) : + ∃ z, ExecWithin w splitBody y z 11 ∧ I z ∧ L - z.regs level < L - y.regs level := by + obtain ⟨z, he, hz, hm, ht⟩ := exec_splitBody L a y hy.1 hc hfit + exact ⟨z, he, ⟨hz, hm.trans hy.2⟩, ht⟩ + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => L - y.regs level) hb ⟨hi, rfl⟩ + have hlevel : y.regs level = L := by + simp only [Cond.holds_lt, hy.1.fixed.levels] at hc + have := hy.1.level_le + omega + have hwhole : ExecWithin w splitCode x y (13 * L + 10) := by + apply (hs.seq he).mono + have hs0 : start.regs level = 0 := by simp [start] + rw [hs0, Nat.sub_zero] + omega + refine ⟨y, hwhole, hy.1.fixed, ?_, ?_, ?_, hy.2⟩ + · simpa only [hlevel, Nat.sub_self, parts, Nat.mul_zero, Nat.add_zero] + using hy.1.mask_value + · simpa only [hlevel, Nat.sub_self, parts, Nat.mul_zero, Nat.add_zero] + using hy.1.outer_value + · simpa only [hlevel, Nat.sub_self, parts, Nat.mul_zero, Nat.add_zero] + using hy.1.inner_value + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandTables.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandTables.lean new file mode 100644 index 0000000000..41bf96539a --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandTables.lean @@ -0,0 +1,280 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandPopulation + +/-! +# Charged materialization of the subset tables + +## Main definitions + +* `BandPreparation.tableStep`: the exact memory writes made by one table iteration. + +## Main results + +* `BandPreparation.exec_tableWrite`: actual inverse/forward table stores and counter updates. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open MaskTable + +/-- The inverse word is initialized even for an unselected mask. Selected masks also write +one forward-table cell. This is only a mathematical description of the two native stores. -/ +def tableStep (L m inv tab t : ℕ) (μ : ℕ → ℕ) : ℕ → ℕ := + let μ' := Function.update μ (inv + t) (lookup L m t) + if lookup L m t = 0 then μ' else Function.update μ' (tab + prefixCount L m t) t + +/-- The part of both tables computed after visiting `t` masks. The outside frame permits +arbitrary caller memory, including future query inputs and scratch regions. -/ +structure TablesPrefix (L m inv tab t : ℕ) (initial current : ℕ → ℕ) : Prop where + /-- Every visited inverse word has been initialized. -/ + inverse : ∀ a < t, current (inv + a) = lookup L m a + /-- Every selected visited mask has been placed at its prefix rank. -/ + forward : ∀ a < t, lookup L m a ≠ 0 → current (tab + prefixCount L m a) = a + /-- Only the two explicitly reserved table blocks may change. -/ + frame : ∀ a, ¬(inv ≤ a ∧ a < inv + 2 ^ L) → ¬(tab ≤ a ∧ a < tab + 2 ^ L) → + current a = initial a + +/-- Empty scan prefix, with no initialized-memory assumption. -/ +theorem TablesPrefix.zero (L m inv tab : ℕ) (μ : ℕ → ℕ) : + TablesPrefix L m inv tab 0 μ μ := + ⟨fun _ h => (Nat.not_lt_zero _ h).elim, + fun _ h => (Nat.not_lt_zero _ h).elim, fun _ _ _ => rfl⟩ + +/-- The actual two stores extend the initialized table prefix by one mask. -/ +theorem TablesPrefix.step {L m inv tab t : ℕ} {μ μ' : ℕ → ℕ} + (h : TablesPrefix L m inv tab t μ μ') (ht : t < 2 ^ L) + (hsep : inv + 2 ^ L ≤ tab) : + TablesPrefix L m inv tab (t + 1) μ (tableStep L m inv tab t μ') := by + have hr : prefixCount L m t < 2 ^ L := (prefixCount_le L m t).trans_lt ht + constructor + · intro a ha + by_cases hat : a = t + · subst a + unfold tableStep + split + · exact Function.update_self _ _ _ + · rw [Function.update_of_ne (by omega), Function.update_self] + · have haold : a < t := by omega + unfold tableStep + split + · rw [Function.update_of_ne (by omega)] + exact h.inverse a haold + · rw [Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + exact h.inverse a haold + · intro a ha hselected + by_cases hat : a = t + · subst a + simp only [tableStep, hselected, ↓reduceIte, Function.update_self] + · have haold : a < t := by omega + have hpop := (lookup_ne_zero_iff L m a).mp hselected |>.2.1 + have hless := prefixCount_lt hpop haold + unfold tableStep + split + · rw [Function.update_of_ne (by omega)] + exact h.forward a haold hselected + · rw [Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + exact h.forward a haold hselected + · intro a hi hf + have hia : a ≠ inv + t := by omega + have hta : a ≠ tab + prefixCount L m t := by omega + unfold tableStep + split + · rw [Function.update_of_ne hia] + exact h.frame a hi hf + · rw [Function.update_of_ne hta, Function.update_of_ne hia] + exact h.frame a hi hf + +/-- Execute the stores after the population loop. The code neither searches a supplied +table nor receives a table-validity certificate. -/ +theorem exec_tableWrite {w : ℕ} (L m inv tab t : ℕ) (x : Data 96) + (ht : t < 2 ^ L) (hinv : x.regs inverse = inv) (htab : x.regs table = tab) + (hmask : x.regs mask = t) (hrank : x.regs rank = prefixCount L m t) + (hpop : x.regs population = MaskTable.population L t) (hm : x.regs inner = m) + (hsq : x.regs square = K0 L m * K0 L m) + (h0 : x.regs zero = 0) (h1 : x.regs one = 1) + (hfit : Fits w (inv + tab + 2 ^ L + 4)) : + ∃ y, ExecWithin w tableWrite x y 10 ∧ y.mem = tableStep L m inv tab t x.mem ∧ + y.regs mask = t + 1 ∧ y.regs rank = prefixCount L m (t + 1) ∧ + (∀ i, i ≠ address → i ≠ temporary → i ≠ mask → i ≠ rank → y.regs i = x.regs i) := by + have hbound := prefixCount_le L m t + have ha : Fits w (inv + t) := hfit.of_le (by omega) + have htaddr : Fits w (tab + prefixCount L m t) := hfit.of_le (by omega) + have hr : Fits w (prefixCount L m t + 1) := hfit.of_le (by omega) + have ht1 : Fits w (t + 1) := hfit.of_le (by omega) + by_cases hp : MaskTable.population L t = m + · by_cases hselected : prefixCount L m t < K0 L m * K0 L m + · let y := (((((((x.set address (inv + t)).write (inv + t) 0).set temporary + (prefixCount L m t + 1)).write (inv + t) (prefixCount L m t + 1)).set + address (tab + prefixCount L m t)).write (tab + prefixCount L m t) t).set + rank (prefixCount L m t + 1)).set mask (t + 1) + refine ⟨y, ?_, ?_, ?_, ?_, ?_⟩ + · ram_eval [tableWrite, y, hinv, htab, hmask, hrank, hpop, hm, hsq, h0, h1, + ha, htaddr, hr, ht1, hp, hselected, + inverse, table, mask, rank, population, inner, square, zero, one, address, + temporary, Fin.mk.injEq] + · simp [y, tableStep, lookup, ht, hp, hselected, Data.write, Data.set, + Function.update_idem] + · simp [y, Data.set, address, temporary, mask, rank] + · simp [y, Data.set, address, temporary, mask, rank, prefixCount_succ, hp] + · intro i ha ht hm hr + simp only [y, Data.set_regs, Data.write_regs, + Function.update_of_ne ha, Function.update_of_ne ht, + Function.update_of_ne hm, Function.update_of_ne hr] + · let y := ((((x.set address (inv + t)).write (inv + t) 0).set rank + (x.regs rank)).set rank (prefixCount L m t + 1)).set mask (t + 1) + refine ⟨y, ?_, ?_, ?_, ?_, ?_⟩ + · ram_eval [tableWrite, y, hinv, htab, hmask, hrank, hpop, hm, hsq, h0, h1, + ha, ht1, hr, hp, hselected, inverse, table, mask, rank, population, inner, + square, zero, one, address, temporary, Fin.mk.injEq] + · simp [y, tableStep, lookup, ht, hp, hselected, Data.write, Data.set] + · simp [y, Data.set, address, mask, rank] + · simp [y, Data.set, address, mask, rank, prefixCount_succ, hp] + · intro i ha _ hm hr + simp only [y, Data.set_regs, Data.write_regs, Function.update_of_ne ha, + Function.update_of_ne hm, Function.update_of_ne hr] + · let y := (((x.set address (inv + t)).write (inv + t) 0).set rank + (x.regs rank)).set mask (t + 1) + refine ⟨y, ?_, ?_, ?_, ?_, ?_⟩ + · ram_eval [tableWrite, y, hinv, htab, hmask, hrank, hpop, hm, hsq, h0, h1, + ha, ht1, hp, inverse, mask, rank, population, inner, zero, one, + address, Fin.mk.injEq] + · simp [y, tableStep, lookup, ht, hp, Data.write, Data.set] + · simp [y, Data.set, address, mask, rank] + · simp [y, Data.set, address, mask, rank, prefixCount_succ, hp, hrank] + · intro i ha _ hm hr + simp only [y, Data.set_regs, Data.write_regs, Function.update_of_ne ha, + Function.update_of_ne hm, Function.update_of_ne hr] + +/-- Persistent scalar parameters read by the table loops. These are register contents only; +there is no execution, fit or cost assumption in this postcondition. -/ +structure TableParameters (L m inv tab : ℕ) (x : Data 96) : Prop where + /-- Original depth. -/ + levels : x.regs levels = L + /-- Original inner cardinality. -/ + inner : x.regs inner = m + /-- Computed complete binary-mask range. -/ + masks : x.regs masks = 2 ^ L + /-- Inverse-table base. -/ + inverse : x.regs inverse = inv + /-- Forward-table base. -/ + table : x.regs table = tab + /-- Computed number of selected masks. -/ + square : x.regs square = K0 L m * K0 L m + /-- Constant zero. -/ + zero : x.regs zero = 0 + /-- Constant one. -/ + one : x.regs one = 1 + /-- Constant two. -/ + two : x.regs two = 2 + +private def tableParameterRegisters : Finset (Fin 96) := + {levels, inner, masks, inverse, table, square, zero, one, two} + +private theorem TableParameters.of_eq {L m inv tab : ℕ} {x y : Data 96} + (hx : TableParameters L m inv tab x) + (h : ∀ i ∈ tableParameterRegisters, y.regs i = x.regs i) : + TableParameters L m inv tab y := by + exact ⟨(h _ (by decide)).trans hx.levels, (h _ (by decide)).trans hx.inner, + (h _ (by decide)).trans hx.masks, (h _ (by decide)).trans hx.inverse, + (h _ (by decide)).trans hx.table, (h _ (by decide)).trans hx.square, + (h _ (by decide)).trans hx.zero, (h _ (by decide)).trans hx.one, + (h _ (by decide)).trans hx.two⟩ + +/-- The complete native scan computes both tables, with all inner population scans and +outer controls charged. Outside memory may contain arbitrary caller data. -/ +theorem exec_tables {w : ℕ} (L m inv tab : ℕ) (x : Data 96) + (hx : TableParameters L m inv tab x) (hsep : inv + 2 ^ L ≤ tab) + (hfit : Fits w (inv + tab + 2 ^ L + L + 4)) : + ∃ y, ExecWithin w tablesCode x y ((6 * L + 17) * 2 ^ L + 4) ∧ + TableParameters L m inv tab y ∧ + (∀ a < 2 ^ L, y.mem (inv + a) = lookup L m a) ∧ + (∀ gh : Fin (K0 L m) × Fin (K0 L m), + y.mem (tab + (pairRank gh).val) = (selectedMask gh).val) ∧ + (∀ a, ¬(inv ≤ a ∧ a < inv + 2 ^ L) → ¬(tab ≤ a ∧ a < tab + 2 ^ L) → + y.mem a = x.mem a) := by + have h0 : Fits w 0 := hfit.of_le (Nat.zero_le _) + let start := (x.set mask 0).set rank 0 + have hs : ExecWithin (B := .standard) w (ram[ mask ← #0; rank ← #0 ]) x start 2 := by + ram_eval [start, h0, mask, rank, Fin.mk.injEq] + have hstart : TableParameters L m inv tab start := hx.of_eq (by + intro i hi + have hf : ∀ i ∈ tableParameterRegisters, i ≠ mask ∧ i ≠ rank := by decide + simp only [start, Data.set_regs, Function.update_of_ne (hf i hi).1, + Function.update_of_ne (hf i hi).2]) + let I (y : Data 96) : Prop := + TableParameters L m inv tab y ∧ y.regs mask ≤ 2 ^ L ∧ + y.regs rank = prefixCount L m (y.regs mask) ∧ + TablesPrefix L m inv tab (y.regs mask) x.mem y.mem + have hi : I start := by + refine ⟨hstart, ?_, ?_, ?_⟩ + · simp [start, mask, rank] + · simp [start, mask, rank, prefixCount_zero] + · simpa [start, mask, rank] using TablesPrefix.zero L m inv tab x.mem + have hb (y : Data 96) (hy : I y) (hc : (Cond.lt mask masks).Holds y) : + ∃ z, ExecWithin w tableBody y z (6 * L + 15) ∧ I z ∧ + 2 ^ L - z.regs mask < 2 ^ L - y.regs mask := by + obtain ⟨hp, ht, hrank, hmem⟩ := hy + have htlt : y.regs mask < 2 ^ L := by simpa only [Cond.holds_lt, hp.masks] using hc + obtain ⟨u, hu, hpop, humem, huf⟩ := exec_population L (y.regs mask) y hp.levels rfl + hp.one hp.two (hfit.of_le (by omega)) + have hup : TableParameters L m inv tab u := hp.of_eq (by + intro i hi + have hf : ∀ i ∈ tableParameterRegisters, + i ≠ residue ∧ i ≠ population ∧ i ≠ counter ∧ i ≠ digit := by decide + exact huf i (hf i hi).1 (hf i hi).2.1 (hf i hi).2.2.1 (hf i hi).2.2.2) + have hut : u.regs mask = y.regs mask := + huf mask (by decide) (by decide) (by decide) (by decide) + have hur : u.regs rank = prefixCount L m (y.regs mask) := + (huf rank (by decide) (by decide) (by decide) (by decide)).trans hrank + obtain ⟨z, hz, hzmem, hzt, hzr, _⟩ := exec_tableWrite L m inv tab (y.regs mask) u + htlt hup.inverse hup.table hut hur hpop hup.inner hup.square hup.zero hup.one + (hfit.of_le (by omega)) + have he : ExecWithin w tableBody y z (6 * L + 15) := hu.seq hz + have hzp : TableParameters L m inv tab z := hp.of_eq (by + intro i hi + apply he.regs_eq + have hf : ∀ i ∈ tableParameterRegisters, tableBody.Preserves i := by + simp only [tableBody, tableWrite, populationCode, Code.Preserves, Op.dest, Op.bin] + decide + exact hf i hi) + refine ⟨z, he, ⟨hzp, ?_, ?_, ?_⟩, ?_⟩ + · rw [hzt] + omega + · simpa only [hzt] using hzr + · rw [hzt, hzmem, humem] + exact hmem.step htlt hsep + · rw [hzt] + omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => 2 ^ L - y.regs mask) hb hi + have ht : y.regs mask = 2 ^ L := by + simp only [Cond.holds_lt, hy.1.masks] at hc + have := hy.2.1 + omega + have hwhole : ExecWithin w tablesCode x y ((6 * L + 17) * 2 ^ L + 4) := by + apply (hs.seq he).mono + have hs0 : start.regs mask = 0 := by simp [start, mask, rank] + rw [hs0, Nat.sub_zero] + nlinarith + refine ⟨y, hwhole, hy.1, ?_, ?_, hy.2.2.2.frame⟩ + · intro a ha + exact hy.2.2.2.inverse a (ht ▸ ha) + · intro gh + have hlookup : lookup L m (selectedMask gh).val ≠ 0 := by rw [lookup_selectedMask]; omega + have hb : (selectedMask gh).val < y.regs mask := by rw [ht]; exact (selectedMask gh).isLt + have h := hy.2.2.2.forward (selectedMask gh).val hb hlookup + simpa only [prefixCount_selectedMask] using h + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandThin.lean b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandThin.lean new file mode 100644 index 0000000000..eb7c622005 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/BandThin.lean @@ -0,0 +1,76 @@ +/- +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.WordRAM.MatrixMultiplication.Schonhage.BandFlat + +/-! +# Thin-regime native instruction bound + +## Main results + +* `BandPreparation.runs_thin`: actual initialized native execution, including the halt, + within `81920*(L+1)*N²/(4^m)^(1/18)` at `L=19m`, `N≥(4^m)^18`. +* `BandPreparation.width_exists`: original magnitude and parameter inequalities are jointly + satisfiable for every finite input bound, including positive thin parameters. + +The result concerns shared-band preparation and encoding only. Pruned queries, products and +the whole ordinary 3SUM endpoint remain separate obligations. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation + +open Cslib.WordRAM +open Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- All charged preparation and encoder instructions, including the final halt, satisfy +the finite thin-regime saving. Both the coefficient and polynomial exponent are fixed. -/ +theorem runs_thin {w m N : ℕ} (hw : 0 < w) (hN : D m ^ 18 ≤ N) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (A : ℤ) (hA : 0 ≤ A) (hX : ∀ i k, |X i k| ≤ A) (hY : ∀ k j, |Y k j| ≤ A) + (hmag : (7 ^ (19 * m) : ℕ) * A < (2 : ℤ) ^ (w - 1)) + (hfit : Fits (w + 2) (wordBudget (19 * m) m N)) : + ∃ n, Executes (w + 2) (program.program catalogue) + (State.init (input (19 * m) m N X Y)) (catalogueBase (19 * m) m N) n ∧ + (n : ℝ) ≤ 81920 * ((19 * m + 1 : ℕ) : ℝ) * (N : ℝ) ^ 2 / + (D m : ℝ) ^ (1 / 18 : ℝ) := by + obtain ⟨_, _, n, hn, _, he, _⟩ := runs hw (show m ≤ 19 * m by omega) + X Y A hA hX hY hmag hfit + refine ⟨n + 1, he, ?_⟩ + have hpos : 1 ≤ timeBound (19 * m) m N := by + apply Nat.succ_le_of_lt + unfold timeBound + positivity + have hc : n + 1 ≤ 2 * timeBound (19 * m) m N := by omega + have hcR : ((n + 1 : ℕ) : ℝ) ≤ 2 * (timeBound (19 * m) m N : ℝ) := by + exact_mod_cast hc + have ht := thin_timeBound_div m N hN + calc + ((n + 1 : ℕ) : ℝ) ≤ 2 * (timeBound (19 * m) m N : ℝ) := hcR + _ ≤ 2 * (40960 * ((19 * m + 1 : ℕ) : ℝ) * (N : ℝ) ^ 2 / + (D m : ℝ) ^ (1 / 18 : ℝ)) := mul_le_mul_of_nonneg_left ht (by positivity) + _ = _ := by ring + +/-- A concrete, deliberately generous width witnesses simultaneous signed-magnitude and +unsigned-address fits. This proves nonvacuity without enumerating a thin matrix. -/ +theorem width_exists (L m N : ℕ) (A : ℤ) (hA : 0 ≤ A) : + ∃ w, 0 < w ∧ (7 ^ L : ℕ) * A < (2 : ℤ) ^ (w - 1) ∧ + Fits (w + 2) (wordBudget L m N) := by + let q := wordBudget L m N + 7 ^ L * A.natAbs + 3 + refine ⟨q + 1, by omega, ?_, ?_⟩ + · have hnat : 7 ^ L * A.natAbs < 2 ^ q := + (Nat.lt_two_pow_self (n := 7 ^ L * A.natAbs)).trans_le + (Nat.pow_le_pow_right (by decide) (by dsimp [q]; omega)) + have hz : ((7 ^ L * A.natAbs : ℕ) : ℤ) < (2 : ℤ) ^ q := by exact_mod_cast hnat + simpa only [Nat.add_sub_cancel, Nat.cast_mul, Int.natCast_natAbs, + abs_of_nonneg hA] using hz + · exact (fits_self (wordBudget L m N)).mono (by dsimp [q]; omega) + +end Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandPreparation diff --git a/Complexity/Computability/WordRAM/Power.lean b/Complexity/Computability/WordRAM/Power.lean new file mode 100644 index 0000000000..a49be46a73 --- /dev/null +++ b/Complexity/Computability/WordRAM/Power.lean @@ -0,0 +1,93 @@ +/- +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.Computability.WordRAM.Notation + +/-! +# Native natural powers + +## Main definitions + +* `Power.program`: repeated native multiplication, with base and exponent in registers. + +## Main results + +* `Power.exec`: `b^n` within `4n+5` instructions, with input-derived intermediate fits. +-/ + +@[expose] public section + +namespace Cslib.WordRAM.Power + +/-- Result accumulator, initially one. -/ +abbrev value : Fin 5 := 0 +/-- Number of completed multiplications. -/ +abbrev counter : Fin 5 := 1 +/-- Input base. -/ +abbrev base : Fin 5 := 2 +/-- Input exponent. -/ +abbrev exponent : Fin 5 := 3 +/-- Constant one, initialized by the routine. -/ +abbrev one : Fin 5 := 4 + +/-- One multiplication and one counter increment. -/ +def body : Code .standard 5 := ram[ value ← value * base; counter ← counter + one ] + +/-- A fixed program for natural exponentiation by repeated multiplication. -/ +def program : Code .standard 5 := + .seq (ram[ value ← #1; counter ← #0; one ← #1 ]) (.loop (.lt counter exponent) body) + +/-- The routine is uniform in its arguments and preserves memory and both input registers. +Every multiplication fit follows from the final power and exponent bound. -/ +theorem exec {w : ℕ} (b n : ℕ) (hb : 1 ≤ b) (x : Data 5) + (hx : x.regs base = b) (hn : x.regs exponent = n) (hfit : Fits w (b ^ n + n + 2)) : + ∃ y, ExecWithin w program x y (4 * n + 5) ∧ y.regs value = b ^ n ∧ + y.regs counter = n ∧ y.regs one = 1 ∧ y.mem = x.mem ∧ + y.regs base = b ∧ y.regs exponent = n := by + have h0 : Fits w 0 := hfit.of_le (Nat.zero_le _) + have h1 : Fits w 1 := hfit.of_le (by omega) + let start := ((x.set value 1).set counter 0).set one 1 + have hs : ExecWithin (B := .standard) w + (ram[ value ← #1; counter ← #0; one ← #1 ]) x start 3 := by + ram_eval [start, h0, h1, value, counter, one, Fin.mk.injEq] + let I (y : Data 5) : Prop := + y.regs base = b ∧ y.regs exponent = n ∧ y.regs one = 1 ∧ y.regs counter ≤ n ∧ + y.regs value = b ^ y.regs counter ∧ y.mem = x.mem + have hi : I start := by simp [I, start, hx, hn, base, exponent, one, counter, value] + have hbody (y : Data 5) (hy : I y) (hc : (Cond.lt counter exponent).Holds y) : + ∃ z, ExecWithin w body y z 2 ∧ I z ∧ n - z.regs counter < n - y.regs counter := by + obtain ⟨hbv, hnv, ho, hc_le, hv, hm⟩ := hy + have hc_lt : y.regs counter < n := by simpa only [Cond.holds_lt, hnv] using hc + have hpn : b ^ (y.regs counter + 1) ≤ b ^ n := Nat.pow_le_pow_right hb (by omega) + have hp : Fits w (b ^ (y.regs counter + 1)) := hfit.of_le (hpn.trans (by omega)) + have hj : Fits w (y.regs counter + 1) := hfit.of_le (by omega) + let z := (y.set value (b ^ (y.regs counter + 1))).set counter (y.regs counter + 1) + have he : ExecWithin w body y z 2 := by + ram_eval [body, z, hbv, ho, hv, ← pow_succ, hp, hj, + value, counter, base, one, Fin.mk.injEq] + refine ⟨z, he, ?_, ?_⟩ + · simpa +decide [I, z, base, exponent, one, counter, value, Function.update_apply] using + (show y.regs base = b ∧ y.regs exponent = n ∧ y.regs one = 1 ∧ + y.regs counter + 1 ≤ n ∧ y.mem = x.mem from ⟨hbv, hnv, ho, by omega, hm⟩) + · simp only [z, Data.set_regs, Function.update_self] + omega + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => n - y.regs counter) hbody hi + have hj : y.regs counter = n := by + simp only [Cond.holds_lt, hy.2.1] at hc + have := hy.2.2.2.1 + omega + refine ⟨y, ?_, ?_, hj, hy.2.2.1, hy.2.2.2.2.2, hy.1, hy.2.1⟩ + · apply (hs.seq he).mono + have hs0 : start.regs counter = 0 := by simp [start, counter, one] + rw [hs0, Nat.sub_zero] + omega + · simpa only [hj] using hy.2.2.2.2.1 + +end Cslib.WordRAM.Power diff --git a/Complexity/Foundations/Data/Finset/OrderedFilter.lean b/Complexity/Foundations/Data/Finset/OrderedFilter.lean new file mode 100644 index 0000000000..81ed3ec1f9 --- /dev/null +++ b/Complexity/Foundations/Data/Finset/OrderedFilter.lean @@ -0,0 +1,42 @@ +/- +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.Init +public import Mathlib.Data.Finset.Sort +public import Mathlib.Data.List.OfFn + +/-! +# Ordered finite-set enumeration and filtering + +## Main results + +* `Finset.sort_eq_filter_finRange`: a level set in increasing order is the filtered level list. +* `Finset.ofFn_orderEmbOfFin`: values at its ordered positions form its sorted mapped list. +-/ + +@[expose] public section + +namespace Finset + +/-- Increasing level-set enumeration is ordinary filtering of all levels. -/ +theorem sort_eq_filter_finRange {L : ℕ} (s : Finset (Fin L)) : + s.sort = (List.finRange L).filter fun i => decide (i ∈ s) := by + apply s.sortedLT_sort.eq_of_mem_iff + ((List.sortedLT_finRange L).pairwise.filter _ |>.sortedLT) + intro a + simp + +/-- Applying a function to the ordered finite-set positions agrees with its sorted list. -/ +theorem ofFn_orderEmbOfFin {α β : Type*} [LinearOrder α] (s : Finset α) + {k : ℕ} (h : s.card = k) (f : α → β) : + List.ofFn (fun i => f (s.orderEmbOfFin h i)) = s.sort.map f := by + rw [← s.listMap_orderEmbOfFin_finRange h, List.map_map, List.ofFn_eq_map] + rfl + +end Finset diff --git a/Complexity/Foundations/Data/List/RowMajor.lean b/Complexity/Foundations/Data/List/RowMajor.lean new file mode 100644 index 0000000000..c03ab40239 --- /dev/null +++ b/Complexity/Foundations/Data/List/RowMajor.lean @@ -0,0 +1,64 @@ +/- +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.Init +public import Mathlib.Data.List.OfFn +public import Mathlib.Data.List.FinRange +public import Mathlib.Logic.Equiv.Fin.Basic +import Mathlib.Data.List.Basic +import Mathlib.Data.List.GetD + +/-! +# Row-major finite arrays as lists + +## Main results + +* `List.flatMap_ofFn_prod`: nested row-major enumeration agrees with the existing product + equivalence on finite indices. +* `List.getD_flatMap_ofFn`: read a cell at its ordinary flat row-major index. +-/ + +@[expose] public section + +namespace List + +/-- Flattening all rows enumerates the standard finite product equivalence. -/ +theorem flatMap_ofFn_prod {α : Type*} {n d : ℕ} (f : Fin n → Fin d → α) : + (finRange n).flatMap (fun i => ofFn (f i)) = + ofFn (fun p : Fin (n * d) => f (finProdFinEquiv.symm p).1 + (finProdFinEquiv.symm p).2) := by + rw [ofFn_mul] + simp only [ofFn_eq_map, ← flatMap_def] + apply flatMap_congr + intro i _ + apply map_congr_left + intro j _ + have heq : (⟨i.val * d + j.val, by + calc + i.val * d + j.val < i.val * d + d := Nat.add_lt_add_left j.isLt _ + _ = (i.val + 1) * d := by rw [Nat.add_mul, Nat.one_mul] + _ ≤ n * d := Nat.mul_le_mul_right d i.isLt⟩ : Fin (n * d)) = finProdFinEquiv (i, j) := by + apply Fin.ext + simp [finProdFinEquiv, Nat.mul_comm, Nat.add_comm] + simp only [heq, Equiv.symm_apply_apply] + +/-- Ordinary flat indexing reads the corresponding original row and column. -/ +theorem getD_flatMap_ofFn {α : Type*} {n d : ℕ} (f : Fin n → Fin d → α) + (default : α) (i : Fin n) (j : Fin d) : + ((finRange n).flatMap (fun i => ofFn (f i))).getD (i.val * d + j.val) default = + f i j := by + rw [flatMap_ofFn_prod] + have heq : i.val * d + j.val = (finProdFinEquiv (i, j)).val := by + simp [finProdFinEquiv, Nat.mul_comm, Nat.add_comm] + rw [heq, getD_eq_getElem _ _ (by + rw [length_ofFn] + exact (finProdFinEquiv (i, j)).isLt), List.getElem_ofFn] + exact congrArg (fun p => f p.1 p.2) (finProdFinEquiv.symm_apply_apply (i, j)) + +end List diff --git a/Complexity/Foundations/Data/Nat/OfDigits.lean b/Complexity/Foundations/Data/Nat/OfDigits.lean new file mode 100644 index 0000000000..9fd7ecce81 --- /dev/null +++ b/Complexity/Foundations/Data/Nat/OfDigits.lean @@ -0,0 +1,32 @@ +/- +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.Foundations.Data.Nat.Radix + +/-! +# Finite-function digit numbering + +## Main results + +* `Nat.ofDigits_finFunctionFinEquiv`: agreement with Mathlib's finite-function numbering. + +The weighted-digit identity is the existing `Nat.ofDigits_ofFn_eq_sum` from fixed-length radix +arithmetic; this module supplies only its bridge to Mathlib's finite-function enumeration. +-/ + +@[expose] public section + +namespace Nat + +/-- `finFunctionFinEquiv` uses exactly `ofDigits`'s least-significant-first convention. -/ +theorem ofDigits_finFunctionFinEquiv {b L : ℕ} (f : Fin L → Fin b) : + ofDigits b (List.ofFn fun i => (f i).val) = (finFunctionFinEquiv f).val := + ofDigits_ofFn_eq_sum b _ + +end Nat diff --git a/ComplexityTests/SchonhageBandPreparation.lean b/ComplexityTests/SchonhageBandPreparation.lean new file mode 100644 index 0000000000..bd2d3af56b --- /dev/null +++ b/ComplexityTests/SchonhageBandPreparation.lean @@ -0,0 +1,167 @@ +/- +Copyright (c) 2026 Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Erik Demaine +-/ + +import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandThin +import Complexity.Algorithms.WordRAM.MatrixMultiplication.Schonhage.BandCatalogue +import Complexity.Algorithms.MatrixMultiplication.Schonhage.MaskTable + +-- These regressions evaluate actual native code; general execution proofs are checked separately. +set_option linter.hashCommand false + +namespace ComplexityTests.SchonhageBandPreparation + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.MatrixMultiplication.Schonhage +open Cslib.Algorithms.WordRAM.MatrixMultiplication.Schonhage +open BandPreparation + +private def leafAt {L : ℕ} (i : Fin (10 ^ L)) : Leaf L := + fun ℓ => Encoding.termAt (finFunctionFinEquiv.symm i ℓ) + +private def expected {L m N : ℕ} (hm : m ≤ L) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) : List ℤ := + ((List.range (numBands L m N)).flatMap fun β => + List.ofFn fun i => encodingL (bandArrayL (MaskTable.layout hm) X β) (leafAt i)) ++ + ((List.range (numBands L m N)).flatMap fun β => + List.ofFn fun i => encodingR (bandArrayR (MaskTable.layout hm) Y β) (leafAt i)) + +private def observe (w L m N : ℕ) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) : + Option ((List ℤ × Bool) × List ℕ) := + let initial := ((State.init (input L m N X Y) : State 96).data.write + (memoryEnd L m N) 91) + (program.eval w 1024 initial).map fun (y, _) => + (((List.range (2 * numBands L m N * 10 ^ L)).map fun i => + decode (y.mem (catalogueBase L m N + i)), + decide ((∀ a < inputEnd m N, y.mem a = initial.mem a) ∧ + y.mem (memoryEnd L m N) = 91)), + (List.range (2 ^ L)).map fun i => y.mem (inverseBase L m N + i)) + +-- Distinct tiles share the existing row/column band blocks, with no per-tile duplication. +#guard tileInputBases 1 0 4 (0, 1) == (57, 87) +#guard tileInputBases 1 0 4 (1, 1) == (67, 87) + +-- A zero-dimensional matrix computes its tables and has no bands. Its input is restored. +#guard observe 8 0 0 0 (fun _ _ => 0) (fun _ _ => 0) == some (([], true), [1]) + +-- A genuine signed zero-depth band on each side calls both actual encoder programs. +#guard observe 8 0 0 1 (fun _ _ => -3) (fun _ _ => 2) == some (([-3, 2], true), [1]) + +-- All leaves of four distinct shared bands agree with the unchanged mathematical arrays. +#guard (observe 12 1 1 2 (fun i k => (i.val : ℤ) - k.val) + (fun k j => (k.val : ℤ) - 2 * j.val)).map (fun z => z.1) == + some (expected (N := 2) (by decide : 1 ≤ 1) (fun i k => (i.val : ℤ) - k.val) + (fun k j => (k.val : ℤ) - 2 * j.val), true) + +-- Two row and column bands with padding: the last band has only one real row/column. +#guard (observe 12 1 0 4 (fun i _ => (i.val : ℤ) - 2) + (fun _ j => 1 - (j.val : ℤ))).map (fun z => z.1) == + some (expected (N := 4) (by decide : 0 ≤ 1) (fun i _ => (i.val : ℤ) - 2) + (fun _ j => 1 - (j.val : ℤ)), true) + +-- The second size-one mask is unselected when K0² = 1; the first is mask one. +#guard (observe 12 2 1 1 (fun _ _ => -1) (fun _ _ => 1)).map Prod.snd == + some [0, 1, 0, 0] + +-- Address-width boundary for a padded, two-band call. The nonzero sentinel also survives. +#guard observe 6 1 0 4 (fun _ _ => 0) (fun _ _ => 0) == none +#guard (observe 7 1 0 4 (fun _ _ => 0) (fun _ _ => 0)).map (fun z => z.1.2) == some true + +/-- The original numerical width premise has a checked strict boundary. -/ +example : ¬Fits 13 (wordBudget 1 0 4) ∧ Fits 14 (wordBudget 1 0 4) := by decide +kernel + +/-- A negative zero-depth initialized client. The actual final memory belongs to the +charged prefix run, which is followed by its charged halt. -/ +example : ∃ y : Data 96, ∃ n, n ≤ timeBound 0 0 1 ∧ + Runs 11 (program.program catalogue) + (State.init (input 0 0 1 (fun _ _ => -3) (fun _ _ => 2))) ⟨program.size, y⟩ n ∧ + Executes 11 (program.program catalogue) + (State.init (input 0 0 1 (fun _ _ => -3) (fun _ _ => 2))) + (catalogueBase 0 0 1) (n + 1) ∧ + y.mem (catalogueAddress 0 0 1 false 0 0) = 7 ∧ + y.mem (catalogueAddress 0 0 1 true 0 0) = 4 := by + obtain ⟨_, y, n, hn, hr, he, _, hl, hright, _⟩ := runs (w := 9) (N := 1) (by decide) + (by decide : 0 ≤ 0) (fun _ _ => -3) (fun _ _ => 2) 3 (by decide) + (by intro i k; decide) (by intro k j; decide) (by decide) (by decide +kernel) + refine ⟨y, n, hn, hr, he, ?_, ?_⟩ + · have h := hl 0 (by decide +kernel) (Fin.elim0 : Leaf 0) + have hi : (Encoding.leafIndex (Fin.elim0 : Leaf 0)).val = 0 := by + have := (Encoding.leafIndex (Fin.elim0 : Leaf 0)).isLt + norm_num at this ⊢ + have hv : encodingL (bandArrayL (N := 1) (MaskTable.layout (by decide : 0 ≤ 0)) + (fun _ _ => -3) 0) Fin.elim0 = -3 := by decide +kernel + simpa only [hi, hv, show encode (-3) = 7 from rfl] using h + · have h := hright 0 (by decide +kernel) (Fin.elim0 : Leaf 0) + have hi : (Encoding.leafIndex (Fin.elim0 : Leaf 0)).val = 0 := by + have := (Encoding.leafIndex (Fin.elim0 : Leaf 0)).isLt + norm_num at this ⊢ + have hv : encodingR (bandArrayR (N := 1) (MaskTable.layout (by decide : 0 ≤ 0)) + (fun _ _ => 2) 0) Fin.elim0 = 2 := by decide +kernel + simpa only [hi, hv, show encode 2 = 4 from rfl] using h + +/-- A positive-depth, signed padded multi-band client of the general kernel-checked run. +The prefix is tied to its actual final memory, including restored input and outside frame. -/ +example : ∃ y : Data 96, ∃ n, n ≤ timeBound 1 0 4 ∧ + Runs 14 (program.program catalogue) + (State.init (input 1 0 4 (fun i _ => (i.val : ℤ) - 2) + (fun _ j => 1 - (j.val : ℤ)))) ⟨program.size, y⟩ n ∧ + Executes 14 (program.program catalogue) + (State.init (input 1 0 4 (fun i _ => (i.val : ℤ) - 2) + (fun _ j => 1 - (j.val : ℤ)))) (catalogueBase 1 0 4) (n + 1) ∧ + y.mem (catalogueAddress 1 0 4 false 0 + (Encoding.leafIndex (fun _ : Fin 1 => Term.P 0 0)).val) = encode (-2) ∧ + (∀ b, ¬(inputEnd 0 4 ≤ b ∧ b < memoryEnd 1 0 4) → + y.mem b = (State.init (input 1 0 4 (fun i _ => (i.val : ℤ) - 2) + (fun _ j => 1 - (j.val : ℤ))) : State 96).data.mem b) := by + obtain ⟨_, y, n, hn, hr, he, _, hl, _, _, _, hf⟩ := runs (w := 12) (N := 4) (by decide) + (by decide : 0 ≤ 1) (fun i _ => (i.val : ℤ) - 2) (fun _ j => 1 - (j.val : ℤ)) + 2 (by decide) (by decide +kernel) (by decide +kernel) (by decide) (by decide +kernel) + refine ⟨y, n, hn, hr, he, ?_, hf⟩ + have hv : encodingL (bandArrayL (N := 4) (MaskTable.layout (by decide : 0 ≤ 1)) + (fun i _ => (i.val : ℤ) - 2) 0) (fun _ : Fin 1 => Term.P 0 0) = -2 := by decide +kernel + simpa only [hv] using hl 0 (by decide +kernel) (fun _ : Fin 1 => Term.P 0 0) + +private def framedEmpty : Data 96 := + ((State.init (input 0 0 0 (fun _ _ => 0) (fun _ _ => 0)) : State 96).data.write + 0 77).write (memoryEnd 0 0 0) 91 + +/-- General-call framing allows an unrelated native header zero and a nonzero adjacent +query/scratch sentinel. The full original prefix, including header zero, is restored. -/ +example : ∃ y : Data 96, ∃ n, n ≤ timeBound 0 0 0 ∧ + Runs 10 (program.program catalogue) ⟨0, framedEmpty⟩ ⟨program.size, y⟩ n ∧ + y.mem 0 = 77 ∧ y.mem (memoryEnd 0 0 0) = 91 := by + obtain ⟨y, ⟨n, hn, he⟩, _, _, _, _, hf⟩ := exec (w := 8) (L := 0) (m := 0) + (N := 0) (by decide) (by decide) framedEmpty (by decide +kernel) + (by decide +kernel) (by decide +kernel) (fun _ _ _ => 0) + (by intro side r hr; omega) 0 (by decide) (by intro side r hr; omega) + (by decide) (by decide +kernel) + refine ⟨y, n, hn, ?_, ?_, ?_⟩ + · simpa only [Nat.zero_add] using he.runs (pc := 0) + (p := program.program catalogue) (Program.Contains.self _).append_left + · rw [hf 0 (by decide +kernel)] + decide +kernel + · rw [hf (memoryEnd 0 0 0) (by omega)] + rfl + +/-- Symbolic positive-depth thin parameters and nonzero signed matrices simultaneously +satisfy all width premises and the genuine charged native saving, without enumerating them. -/ +example (m : ℕ) : let M := m + 1 + ∃ w n, 0 < w ∧ 0 < M ∧ Fits (w + 2) (wordBudget (19 * M) M (D M ^ 18)) ∧ + (7 ^ (19 * M) : ℕ) * (1 : ℤ) < (2 : ℤ) ^ (w - 1) ∧ + Executes (w + 2) (program.program catalogue) + (State.init (input (19 * M) M (D M ^ 18) (fun _ _ => 1) (fun _ _ => -1))) + (catalogueBase (19 * M) M (D M ^ 18)) n ∧ + (n : ℝ) ≤ 81920 * ((19 * M + 1 : ℕ) : ℝ) * (D M ^ 18 : ℕ) ^ (2 : ℕ) / + (D M : ℝ) ^ (1 / 18 : ℝ) := by + dsimp only + obtain ⟨w, hw, hmag, hfit⟩ := width_exists (19 * (m + 1)) (m + 1) + (D (m + 1) ^ 18) 1 (by decide) + obtain ⟨n, he, ht⟩ := runs_thin (m := m + 1) (N := D (m + 1) ^ 18) hw (Nat.le_refl _) + (fun _ _ => 1) (fun _ _ => -1) 1 (by decide) + (by intro i k; decide) (by intro k j; decide) hmag hfit + exact ⟨w, n, hw, Nat.succ_pos _, hfit, hmag, he, ht⟩ + +end ComplexityTests.SchonhageBandPreparation diff --git a/NOTICE b/NOTICE index 468bfa04bc..15f6c9aef7 100644 --- a/NOTICE +++ b/NOTICE @@ -31,3 +31,7 @@ under the Apache License, Version 2.0. The shared tensor encoder in Complexity/Algorithms/WordRAM/MatrixMultiplication/Schonhage/ is also part of the Anthropic adaptation described above. + +The native shared-band preparation selectively adapts the band-array and all-band +organization from that same pinned source. Anthropic, PBC credit is retained in +BandProgram.lean; its native mask scans, dense inverse table and execution proofs are new. diff --git a/ROADMAP.md b/ROADMAP.md index c0f55bba8a..8af6cf8ba6 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -441,6 +441,11 @@ each statement against the numbered result. Organize modules by subject, not cha - **Shared native matrix encoders** ([#1493](https://github.com/cslib-community/complexitylib/issues/1493)): compute all `10^L` leaf values from `7^L` signed inputs within `240*(L+1)*10^L` standard WordRAM instructions, plus the final halt, with fixed code/registers and concrete initialized inputs, frames and explicit width bounds. +- **Native shared-band preparation** ([#1523](https://github.com/cslib-community/complexitylib/issues/1523)): computes deterministic subset tables and a valid + layout, prepares padded signed matrix bands, invokes both shared encoders for every band, and + restores the original initialized input. Fixed 96-register code has a charged bound, including + halt, of `4096*(L+1)*(1+N*4^m+2^L+2*numBands*10^L)+1`; at `L=19m`, `N≥(4^m)^18`, it is at most + `81920*(L+1)*N²/(4^m)^(1/18)`. - **Word RAM:** an AC⁰ basis; wider-fragment and Turing-machine simulations; transfer of Anthropic's modular machine and 3SUM/APSP algorithms; direct 3SUM adaptation ([#1394]) has the Convolution-3SUM/Exact Triangle correctness reduction, a charged native single-instance diff --git a/UPSTREAM.md b/UPSTREAM.md index 7803bd2dda..004abdb9d9 100644 --- a/UPSTREAM.md +++ b/UPSTREAM.md @@ -170,6 +170,21 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he blocks in either address order or identical blocks. This is the encoder component only; band/tile preparation and full native composition remain in [#1394](https://github.com/cslib-community/complexitylib/issues/1394). +* **Native shared-band preparation.** The native shared-band preparation component + ([#1523](https://github.com/cslib-community/complexitylib/issues/1523)) computes a deterministic + `Layout` and both subset tables by native mask/popcount scans, computes the square-root block + count by a charged loop, and decodes each band-array index into the original row/column and + thin-dimension coordinates. Real matrix bounds are checked before row-major multiplication. + Each shared row/column band calls the unchanged encoder once; tiles refer to its catalogue + block in the existing leaf order. The fixed 96-register standard code starts from actual + `State.init` scalar/matrix input, charges snapshot and restoration of the complete original + prefix, and preserves arbitrary outside memory. All counter/address fits follow from an + explicit original-parameter word budget; signed fits follow from the original entry bound. + Its complete instruction bound, including halt, is + `81920*(L+1)*N²/(4^m)^(1/18)` when `L=19m` and `N≥(4^m)^18`. Native power, ordered-filter, + finite-digit and row-major list lemmas live in small topic modules; the existing mathematical + layout, matrix and encoding interfaces are reused. No source machine or call/time oracle + is introduced. * **Direct 3SUM adaptation.** `Algorithms/ThreeSum/Convolution.lean` adapts the explicit Convolution-3SUM to Exact Triangle correctness reduction from Anthropic's [`formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp),