diff --git a/Complexity.lean b/Complexity.lean index d94d7824d1..5c3ecc52f1 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -3,7 +3,9 @@ module -- shake: keep-all --deprecated_module: ignore public import Complexity.Algorithms.ExactTriangle.Defs public import Complexity.Algorithms.ExactTriangle.Hashing public import Complexity.Algorithms.ExactTriangle.Lopsided +public import Complexity.Algorithms.ExactTriangle.Splitting public import Complexity.Algorithms.ExactTriangle.Weights +public import Complexity.Algorithms.ExactTriangle.Witnesses public import Complexity.Algorithms.Lean.FusionTree.Capacity public import Complexity.Algorithms.Lean.FusionTree.Compression public import Complexity.Algorithms.Lean.FusionTree.Correction @@ -87,6 +89,7 @@ public import Complexity.Algorithms.WordRAM.FusionTree.Query.Traversal public import Complexity.Algorithms.WordRAM.FusionTree.Rank public import Complexity.Algorithms.WordRAM.FusionTree.ReadHeader public import Complexity.Algorithms.WordRAM.FusionTree.Static +public import Complexity.Algorithms.WordRAM.Histogram public import Complexity.Algorithms.WordRAM.Karger.Contract public import Complexity.Algorithms.WordRAM.Karger.FastCut public import Complexity.Algorithms.WordRAM.Karger.Graph @@ -117,6 +120,9 @@ public import Complexity.Algorithms.WordRAM.Quicksort public import Complexity.Algorithms.WordRAM.Stack public import Complexity.Algorithms.WordRAM.Stack.Controller public import Complexity.Algorithms.WordRAM.ThreeSum.Cell +public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Collisions +public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Count +public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Residues public import Complexity.Algorithms.WordRAM.ThreeSum.Fill public import Complexity.Algorithms.WordRAM.ThreeSum.Memory public import Complexity.Algorithms.WordRAM.ThreeSum.Weight @@ -2100,12 +2106,14 @@ public import Complexity.Foundations.Data.Finset.Average public import Complexity.Foundations.Data.Finset.Expect public import Complexity.Foundations.Data.Finset.Markov public import Complexity.Foundations.Data.Finset.Quantile +public import Complexity.Foundations.Data.Finset.Quotient public import Complexity.Foundations.Data.Finset.Sampling public import Complexity.Foundations.Data.Finset.SubsetSum public import Complexity.Foundations.Data.Finset.SymmDiff public import Complexity.Foundations.Data.Hamming public import Complexity.Foundations.Data.Hamming.Expansion public import Complexity.Foundations.Data.Int.Order +public import Complexity.Foundations.Data.Int.Residue public import Complexity.Foundations.Data.List.Append public import Complexity.Foundations.Data.List.BalancedFold public import Complexity.Foundations.Data.List.Basic @@ -2121,6 +2129,7 @@ public import Complexity.Foundations.Data.List.Range public import Complexity.Foundations.Data.List.SelectDistinct public import Complexity.Foundations.Data.Matrix.Binary public import Complexity.Foundations.Data.Matrix.IntegerBits +public import Complexity.Foundations.Data.Matrix.PairChunks public import Complexity.Foundations.Data.Matrix.Permanent public import Complexity.Foundations.Data.Matrix.Reindex public import Complexity.Foundations.Data.Matrix.TotallyRegular @@ -2143,6 +2152,7 @@ public import Complexity.Foundations.Data.Nat.Bitwise.Projection public import Complexity.Foundations.Data.Nat.Bitwise.Sparse public import Complexity.Foundations.Data.Nat.Bitwise.SumTree public import Complexity.Foundations.Data.Nat.Card +public import Complexity.Foundations.Data.Nat.CeilDiv public import Complexity.Foundations.Data.Nat.Choose public import Complexity.Foundations.Data.Nat.Choose.Real public import Complexity.Foundations.Data.Nat.Div @@ -2157,6 +2167,7 @@ public import Complexity.Foundations.Data.Nat.Prime public import Complexity.Foundations.Data.Nat.Primrec public import Complexity.Foundations.Data.Nat.Radix public import Complexity.Foundations.Data.Nat.Radix.Measurement +public import Complexity.Foundations.Data.Nat.SqrtRatio public import Complexity.Foundations.Data.Option public import Complexity.Foundations.Data.OrderedSet public import Complexity.Foundations.Data.OrderedSet.Spec diff --git a/Complexity/Algorithms/ExactTriangle/Hashing.lean b/Complexity/Algorithms/ExactTriangle/Hashing.lean index fbb2834b6e..17ee536865 100644 --- a/Complexity/Algorithms/ExactTriangle/Hashing.lean +++ b/Complexity/Algorithms/ExactTriangle/Hashing.lean @@ -8,6 +8,7 @@ module -- upstream: candidate public import Complexity.Algorithms.ExactTriangle.Weights +public import Complexity.Foundations.Data.Int.Residue public import Complexity.Foundations.NumberTheory.PrimeWindow public import Complexity.Foundations.RingTheory.CyclicIntRing public import Mathlib.Data.Matrix.Mul @@ -74,15 +75,6 @@ def IsSelectedPrime {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D p : ℕ private theorem four_le_sqrt {D : ℕ} (hD : 16 ≤ D) : 4 ≤ Real.sqrt (D : ℝ) := Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) -private theorem natCast_toNat_emod {p : ℕ} (hp : 0 < p) (u : ℤ) : - ((u % (p : ℤ)).toNat : ℤ) = u % (p : ℤ) := - Int.toNat_of_nonneg (Int.emod_nonneg _ (by omega)) - -private theorem toNat_emod_lt {p : ℕ} (hp : 0 < p) (u : ℤ) : (u % (p : ℤ)).toNat < p := by - have := Int.emod_nonneg u (by omega : (p : ℤ) ≠ 0) - have := Int.emod_lt_of_pos u (by omega : (0 : ℤ) < p) - omega - /-! ### The primes in the range -/ @@ -134,8 +126,8 @@ private theorem eq_add_toNat_emod_iff {p r : ℕ} (hp : p ≠ 0) (hr : r < p) (u rw [← Int.natCast_inj, Int.ModEq, Int.emod_eq_of_lt (Int.natCast_nonneg r) (by exact_mod_cast hr), eq_comm] push_cast - rw [natCast_toNat_emod (Nat.pos_of_ne_zero hp), - natCast_toNat_emod (Nat.pos_of_ne_zero hp), + rw [Int.natCast_toNat_emod (Nat.pos_of_ne_zero hp), + Int.natCast_toNat_emod (Nat.pos_of_ne_zero hp), ← Int.add_emod] /-! @@ -233,8 +225,9 @@ theorem falsePositiveCount_add_zeroCount_eq_sum_coeff {p : ℕ} (hp : p ≠ 0) : refine Finset.sum_congr rfl fun a _ => ?_ rw [Fintype.sum_prod_type] refine Finset.sum_congr rfl fun b _ => ?_ - have hcast := natCast_toNat_emod (Nat.pos_of_ne_zero hp) (-wAB a b) - rw [(coeff_matP_mul_matQ wBC wAC) hp a b _ (toNat_emod_lt (Nat.pos_of_ne_zero hp) _), + have hcast := Int.natCast_toNat_emod (Nat.pos_of_ne_zero hp) (-wAB a b) + have hr : ((-wAB a b) % (p : ℤ)).toNat < p := Int.natMod_lt (m := -wAB a b) hp + rw [(coeff_matP_mul_matQ wBC wAC) hp a b _ hr, Finset.card_filter] push_cast refine Finset.sum_congr rfl fun c _ => ?_ diff --git a/Complexity/Algorithms/ExactTriangle/Splitting.lean b/Complexity/Algorithms/ExactTriangle/Splitting.lean new file mode 100644 index 0000000000..4f6f435c3a --- /dev/null +++ b/Complexity/Algorithms/ExactTriangle/Splitting.lean @@ -0,0 +1,536 @@ +/- +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.ExactTriangle.Hashing +public import Complexity.Algorithms.ExactTriangle.Lopsided +public import Complexity.Foundations.Data.Matrix.PairChunks +public import Complexity.Foundations.Data.Nat.SqrtRatio +import Mathlib.Tactic.GCongr +import Mathlib.Tactic.FieldSimp + +/-! +# Modular splitting of Exact Triangle into finite lopsided queries + +Pieces of the third part and chunks of A–B residue classes define lopsided graphs with a +middle vertex `(c, label)`. The three input matrices retain their original orientations. +All bounds below use `16 ≤ D ≤ n` and `1 ≤ g`; the sharper degree bound also uses `g ≤ sqrt D`. +No native program, fast solver or charged instruction bound is asserted. + +## Main definitions + +* `piece`, `residueClass`, `chunkOf`: vertex blocks and residue query chunks. +* `lopInstance`, `instanceIndices`, `totalChunks`: lopsided graphs and their finite indices. + +## Main results + +* `existsUnique_mem_piece`, `card_piece_le`: unique piece coverage and size. +* `totalChunks_le`, `card_instanceIndices_le`: chunk and instance bounds. +* `card_middle_lopInstance`, `hasCommonNeighbor_lopInstance_iff`: size and modular correctness. +* `ncard_nbr_lopInstance_le`: degree bound. + +## References + +The reduction is part of [AlmanWilliams2026, Section 3]. Direct adaptation of +[Instances.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Instances.lean), with only required definitions from +`PaperStatements.lean`, [Chunks.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Chunks.lean) and [Parameters.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Parameters.lean), at +`anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +Quoted docstring fragments reproduce comments in the pinned formal source; other prose describes +the local Lean API. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.ExactTriangle + +open Cslib.Matrix + +/-- The largest number of query pairs that is "at most n²/√D" (Corollary 15 and Theorem 17; +Theorem 5 has "at most N²/√D"): `⌊n²/√D⌋`. -/ +def queryCap (n D : ℕ) : ℕ := Nat.sqrt ((n ^ 2) ^ 2 / D) + +/-- The integer capacity is exactly the real-valued floor used in the reduction. -/ +theorem queryCap_eq_floor (n D : ℕ) : queryCap n D = ⌊(n : ℝ) ^ 2 / Real.sqrt D⌋₊ := by + rw [queryCap, Nat.sqrt_div_eq_floor_div_sqrt] + norm_cast + +/-- Proof of Theorem 17: "let s := ⌊√D⌋". -/ +def sOf (D : ℕ) : ℕ := Nat.sqrt D + +/-- Proof of Theorem 17: the bound "⌈s/g⌉" on the number of vertices of a piece. -/ +def pieceSize (D g : ℕ) : ℕ := sOf D ⌈/⌉ g + +/-- The integer block size equals the real ceiling, also at the zero-divisor convention. -/ +theorem pieceSize_eq_ceil (D g : ℕ) : pieceSize D g = ⌈(sOf D : ℝ) / (g : ℝ)⌉₊ := by + by_cases hg : g = 0 + · simp [pieceSize, hg] + exact (Nat.ceil_div_eq_ceilDiv _ (Nat.pos_of_ne_zero hg)).symm + +/-- Proof of Theorem 17: "split C into pieces C₁, …, C_h of at most ⌈s/g⌉ vertices each". The paper +does not say how; we take blocks of `⌈s/g⌉` consecutive vertices. `piece n D g k` is the paper's +`C_{k+1}`. -/ +def piece (n D g k : ℕ) : Finset (Fin n) := + Finset.univ.filter fun c : Fin n => c.val / pieceSize D g = k + +/-- Proof of Theorem 17: the number `h` of pieces, `⌈n / ⌈s/g⌉⌉` for the blocks of `piece` +(meaningful when `⌈s/g⌉ ≥ 1`). -/ +def numPieces (n D g : ℕ) : ℕ := n ⌈/⌉ pieceSize D g + +/-- The number of blocks is the ceiling of the real quotient by the block size. -/ +theorem numPieces_eq_ceil (n D g : ℕ) : + numPieces n D g = ⌈(n : ℝ) / (pieceSize D g : ℝ)⌉₊ := by + by_cases h : pieceSize D g = 0 + · simp [numPieces, h] + exact (Nat.ceil_div_eq_ceilDiv _ (Nat.pos_of_ne_zero h)).symm + +/-- The index of an instance built by the reduction of Theorem 17: a residue `ϱ ∈ ℤ_p`, the number +`j` of a chunk of `W_ϱ`, and the number `k` of a piece of `C`. -/ +abbrev InstanceIndex (p : ℕ) : Type := Fin p × ℕ × ℕ + + +/-- Proof of Theorem 17: "For ϱ ∈ ℤ_p let W_ϱ be the set of edges (a,b) ∈ A × B with w(a,b) ≡ ϱ (mod +p)". -/ +def residueClass {n : ℕ} (wAB : Fin n → Fin n → ℤ) (p : ℕ) (ϱ : Fin p) : Finset (Fin n × Fin n) := + Finset.univ.filter fun q : Fin n × Fin n => wAB q.1 q.2 ≡ ((ϱ : ℕ) : ℤ) [ZMOD (p : ℤ)] + +/-- Proof of Theorem 17: "and cut it into chunks of at most n²/√D query pairs". The `j`-th chunk of +`W_ϱ`. -/ +def chunkOf {n : ℕ} (wAB : Fin n → Fin n → ℤ) (D p : ℕ) (ϱ : Fin p) (j : ℕ) : + Finset (Fin n × Fin n) := + chunk (residueClass wAB p ϱ) (queryCap n D) j + +/-- Proof of Theorem 17: "For each chunk 𝒬 ⊆ W_ϱ and each piece C_k form the instance of +Lop-AE-SparseTri(n, D) with query pairs W := 𝒬 and middle part C_k × ℤ_p [...], whose middle +vertices are pairs (vertex, label) [...], with a ∼ (c, σ) ⟺ σ ≡ w(a,c) + ϱ, (c, σ) ∼ b ⟺ σ ≡ +−w(b,c) (mod p)." -/ +def lopInstance {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) : + InstanceIndex p → LopsidedGraph n + | (ϱ, j, k) => + { M := {c : Fin n // c ∈ piece n D g k} × Fin p + adjA := fun a (c, σ) => ((σ : ℕ) : ℤ) ≡ wAC a c + ((ϱ : ℕ) : ℤ) [ZMOD (p : ℤ)] + adjB := fun (c, σ) b => ((σ : ℕ) : ℤ) ≡ -wBC b c [ZMOD (p : ℤ)] + W := chunkOf wAB D p ϱ j } + +/-- Proof of Theorem 17: "For each chunk 𝒬 ⊆ W_ϱ and each piece C_k" there is an instance; these +are the indices of all the instances. This set depends only on the Exact Triangle instance and on +`D`, `g`, `p`, not on any answer of the oracle. -/ +def instanceIndices {n : ℕ} (wAB : Fin n → Fin n → ℤ) (D g p : ℕ) : + Finset (InstanceIndex p) := + Finset.univ.biUnion fun ϱ : Fin p => + ((Finset.range (numChunks (residueClass wAB p ϱ) (queryCap n D))) ×ˢ + (Finset.range (numPieces n D g))).image fun jk => (ϱ, jk) + +/-- Total number of query chunks, summed over residues. -/ +def totalChunks {n : ℕ} (wAB : Fin n → Fin n → ℤ) (D p : ℕ) : ℕ := + ∑ ϱ : Fin p, numChunks (residueClass wAB p ϱ) (queryCap n D) + +variable {n : ℕ} + +/-- `√D ≤ D ≤ n`. -/ +theorem sqrt_le_natCast {D : ℕ} (hD : 1 ≤ D) (hDn : D ≤ n) : Real.sqrt D ≤ n := + calc Real.sqrt D ≤ Real.sqrt D * Real.sqrt D := + le_mul_of_one_le_left (Real.sqrt_nonneg _) (Real.one_le_sqrt.mpr (by exact_mod_cast hD)) + _ = D := Real.mul_self_sqrt (Nat.cast_nonneg D) + _ ≤ n := by exact_mod_cast hDn + +/-- The bound `⌊n²/√D⌋` on the size of a chunk is at least 1, because `√D ≤ n ≤ n²`. -/ +theorem one_le_queryCap {D : ℕ} (hD : 1 ≤ D) (hDn : D ≤ n) : 1 ≤ queryCap n D := by + rw [queryCap_eq_floor] + have hn : (1 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn + refine Nat.le_floor ?_ + rw [Nat.cast_one, le_div_iff₀ (Real.sqrt_pos.mpr (by exact_mod_cast hD)), one_mul] + exact (sqrt_le_natCast hD hDn).trans (le_self_pow₀ hn two_ne_zero) + +private theorem four_le_sqrt {D : ℕ} (hD : 16 ≤ D) : 4 ≤ Real.sqrt (D : ℝ) := + Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) + +variable {n D g p : ℕ} + +/-! ### The pieces -/ + +/-- `s = ⌊√D⌋ ≤ √D`. -/ +theorem sOf_le_sqrt (D : ℕ) : (sOf D : ℝ) ≤ Real.sqrt D := + Real.nat_sqrt_le_real_sqrt + +/-- `√D < s + 1`. -/ +theorem sqrt_lt_sOf_add_one (D : ℕ) : Real.sqrt D < (sOf D : ℝ) + 1 := + Real.real_sqrt_lt_nat_sqrt_succ + +/-- `s > 0` because `D ≥ 16`. -/ +private theorem sOf_pos (hD : 16 ≤ D) : (0 : ℝ) < sOf D := by + linarith [four_le_sqrt hD, sqrt_lt_sOf_add_one D] + +/-- The bound `⌈s/g⌉` on the size of a piece is at least 1. -/ +theorem pieceSize_pos (hD : 16 ≤ D) (hg1 : 1 ≤ g) : 0 < pieceSize D g := by + rw [pieceSize_eq_ceil] + exact Nat.ceil_pos.mpr (div_pos (sOf_pos hD) (by exact_mod_cast hg1)) + +/-- `⌈s/g⌉ ≤ s` because `g ≥ 1`. -/ +theorem pieceSize_le_sOf (hg1 : 1 ≤ g) : pieceSize D g ≤ sOf D := by + rw [pieceSize_eq_ceil] + exact Nat.ceil_le.mpr (div_le_self (Nat.cast_nonneg _) (by exact_mod_cast hg1)) + +/-- A vertex lies in the piece number `k` exactly if its number divided by `⌈s/g⌉` is `k`. -/ +theorem mem_piece {k : ℕ} (c : Fin n) : c ∈ piece n D g k ↔ c.val / pieceSize D g = k := by + simp [piece] + +/-- Proof of Theorem 17: "split C into pieces C₁, …, C_h of at most ⌈s/g⌉ vertices each": every +vertex of `C` lies in exactly one of the `h` pieces. -/ +theorem existsUnique_mem_piece (hD : 16 ≤ D) (hg1 : 1 ≤ g) (c : Fin n) : + ∃! k, k < numPieces n D g ∧ c ∈ piece n D g k := by + refine ⟨c.val / pieceSize D g, ⟨?_, (mem_piece c).mpr rfl⟩, + fun k hk => ((mem_piece c).mp hk.2).symm⟩ + rw [numPieces_eq_ceil] + refine Nat.lt_ceil.mpr (lt_of_le_of_lt Nat.cast_div_le ?_) + exact div_lt_div_of_pos_right (by exact_mod_cast c.isLt) + (by exact_mod_cast pieceSize_pos hD hg1) + +/-- Proof of Theorem 17, "pieces [...] of at most ⌈s/g⌉ vertices each": `|C_k| ≤ ⌈s/g⌉`, because the +numbers of the vertices of the piece number `k` lie in `[k ⌈s/g⌉, (k + 1) ⌈s/g⌉)`. -/ +theorem card_piece_le (hD : 16 ≤ D) (hg1 : 1 ≤ g) (n k : ℕ) : + (piece n D g k).card ≤ pieceSize D g := + Finset.univ.card_filter_div_eq_le (fun _ _ _ _ h => Fin.ext h) (pieceSize_pos hD hg1) k + +/-- Proof of Theorem 17: "so that h ≤ ⌈ng/s⌉". -/ +theorem numPieces_le (hD : 16 ≤ D) (hg1 : 1 ≤ g) : + numPieces n D g ≤ ⌈(n : ℝ) * (g : ℝ) / (sOf D : ℝ)⌉₊ := by + have hsize : (0 : ℝ) < pieceSize D g := by exact_mod_cast pieceSize_pos hD hg1 + have hg : (0 : ℝ) < g := by exact_mod_cast hg1 + -- `s ≤ ⌈s/g⌉ g` + have hs' : (sOf D : ℝ) / (g : ℝ) ≤ pieceSize D g := by + rw [pieceSize_eq_ceil] + exact Nat.le_ceil _ + have hs := (div_le_iff₀ hg).mp hs' + rw [numPieces_eq_ceil] + refine Nat.ceil_mono ?_ + rw [div_le_div_iff₀ hsize (sOf_pos hD)] + calc (n : ℝ) * (sOf D : ℝ) ≤ (n : ℝ) * ((pieceSize D g : ℝ) * (g : ℝ)) := by gcongr + _ = (n : ℝ) * (g : ℝ) * (pieceSize D g : ℝ) := by ring + +/-! ### Residues -/ + +/-- The residue of the integer `x` modulo `p`, as an element of `Fin p`. -/ +def resFin (hp : p ≠ 0) (x : ℤ) : Fin p := + ⟨x.natMod p, Int.natMod_lt hp⟩ + +/-- The residue of `x` is congruent to `x`. -/ +theorem resFin_modEq (hp : p ≠ 0) (x : ℤ) : ((resFin hp x : ℕ) : ℤ) ≡ x [ZMOD (p : ℤ)] := by + change ((x % (p : ℤ)).toNat : ℤ) ≡ x [ZMOD (p : ℤ)] + rw [Int.natCast_toNat_emod (Nat.pos_of_ne_zero hp)] + exact Int.mod_modEq x p + +/-- An element of `Fin p` is congruent to `x` modulo `p` exactly if it is the residue of `x`. -/ +private theorem modEq_iff_eq_resFin (hp : p ≠ 0) (σ : Fin p) (x : ℤ) : + ((σ : ℕ) : ℤ) ≡ x [ZMOD (p : ℤ)] ↔ σ = resFin hp x := by + refine ⟨fun h => ?_, fun h => h ▸ resFin_modEq hp x⟩ + -- Two congruent numbers in `{0, …, p − 1}` are equal. + have hmod : (σ : ℕ) ≡ (resFin hp x : ℕ) [MOD p] := + (Int.natCast_modEq_iff).mp (h.trans (resFin_modEq hp x).symm) + exact Fin.ext (Nat.ModEq.eq_of_lt_of_lt hmod σ.isLt (resFin hp x).isLt) + +/-! ### The arithmetic behind the two counts -/ + +/-- The number `⌈|S|/cap⌉` of chunks of `S` satisfies `⌈|S|/cap⌉ cap ≤ |S| + cap − 1`. -/ +private theorem numChunks_mul_le (S : Finset (Fin n × Fin n)) (cap : ℕ) (hcap : 1 ≤ cap) : + numChunks S cap * cap + 1 ≤ S.card + cap := by + rw [numChunks] + exact Nat.ceilDiv_mul_lt hcap + +/-- `s + 1` exceeds `√D` by at least `1/(2√D + 1)`: since `(s + 1)²` and `D` are integers, +`1 ≤ (s + 1)² − D = (s + 1 − √D)(s + 1 + √D)`. -/ +private theorem one_le_gap_mul (D : ℕ) : + 1 ≤ ((sOf D : ℝ) + 1 - Real.sqrt D) * (2 * Real.sqrt D + 1) := by + have hsq : Real.sqrt D ^ 2 = D := Real.sq_sqrt (Nat.cast_nonneg D) + have hle := sOf_le_sqrt D + have hlt := sqrt_lt_sOf_add_one D + have hinteger : (D : ℝ) + 1 ≤ ((sOf D : ℝ) + 1) ^ 2 := by + have hreal : (D : ℝ) < ((sOf D : ℝ) + 1) ^ 2 := + hsq ▸ pow_lt_pow_left₀ hlt (Real.sqrt_nonneg _) two_ne_zero + have hnat : D < (sOf D + 1) ^ 2 := by exact_mod_cast hreal + exact_mod_cast hnat + calc (1 : ℝ) ≤ ((sOf D : ℝ) + 1) ^ 2 - Real.sqrt D ^ 2 := by linarith [hinteger, hsq] + _ = ((sOf D : ℝ) + 1 - Real.sqrt D) * ((sOf D : ℝ) + 1 + Real.sqrt D) := by ring + _ ≤ ((sOf D : ℝ) + 1 - Real.sqrt D) * (2 * Real.sqrt D + 1) := + mul_le_mul_of_nonneg_left (by linarith) (by linarith) + +/-- The inequality behind "at most p + √D [...] chunks": `n² < (s + 1) ⌊n²/√D⌋ + p`. + +Write `r = √D`, `x = n²/r` and `δ = s + 1 − r`. Then +`(s + 1)⌊x⌋ > (s + 1)(x − 1) = n² + xδ − (s + 1)`, so it is enough that `xδ ≥ r/2 + 1 ≥ s + 1 − p`. +This holds because `x ≥ r³` (as `D ≤ n`), `δ ≥ 1/(2r + 1)` (as `D` is an integer) and `r ≥ 4`. -/ +private theorem sq_lt_mul_queryCap_add (hD : 16 ≤ D) (hDn : D ≤ n) + (hp : Real.sqrt D / 2 ≤ (p : ℝ)) : n ^ 2 < (sOf D + 1) * queryCap n D + p := by + have hgap := one_le_gap_mul D + have hs := sOf_le_sqrt D + set r := Real.sqrt D with hr + set δ := (sOf D : ℝ) + 1 - r with hδ + set x := (n : ℝ) ^ 2 / r with hx + have hr4 : 4 ≤ r := four_le_sqrt hD + have hxr : x * r = (n : ℝ) ^ 2 := div_mul_cancel₀ _ (by linarith) + have hx3 : r ^ 3 ≤ x := by + rw [hx, le_div_iff₀ (by linarith)] + calc r ^ 3 * r = (r ^ 2) ^ 2 := by ring + _ = (D : ℝ) ^ 2 := by rw [hr, Real.sq_sqrt (Nat.cast_nonneg D)] + _ ≤ (n : ℝ) ^ 2 := by gcongr + have hxδ : r / 2 + 1 ≤ x * δ := by + refine le_of_mul_le_mul_right ?_ (show 0 < 2 * r + 1 by linarith) + calc (r / 2 + 1) * (2 * r + 1) ≤ r ^ 3 := by + -- `r³ ≥ 4r²` and `r² ≥ 4r` + linarith [mul_le_mul_of_nonneg_right hr4 (sq_nonneg r), + mul_le_mul_of_nonneg_right hr4 (show 0 ≤ r by linarith)] + _ ≤ x * 1 := by rw [mul_one]; exact hx3 + _ ≤ x * (δ * (2 * r + 1)) := + mul_le_mul_of_nonneg_left hgap (le_trans (by positivity) hx3) + _ = x * δ * (2 * r + 1) := by ring + have hfloor : x < (queryCap n D : ℝ) + 1 := by + rw [queryCap_eq_floor] + exact Nat.lt_floor_add_one _ + have hmain : (n : ℝ) ^ 2 < ((sOf D : ℝ) + 1) * (queryCap n D : ℝ) + p := + calc (n : ℝ) ^ 2 = x * r := hxr.symm + _ ≤ x * r + x * δ - ((sOf D : ℝ) + 1) + r / 2 := by + -- `s + 1 ≤ r + 1 ≤ xδ + r/2` + linarith [hxδ, hs] + _ = ((sOf D : ℝ) + 1) * (x - 1) + r / 2 := by rw [hδ]; ring + _ < ((sOf D : ℝ) + 1) * (queryCap n D : ℝ) + p := + add_lt_add_of_lt_of_le (mul_lt_mul_of_pos_left (by linarith [hfloor]) (by positivity)) hp + exact_mod_cast hmain + +/-- Proof of Theorem 17: the two facts used for the count of the instances, "s ≥ √D − 1 ≥ 3√D/4 +because D ≥ 16, and √D ≤ √n ≤ 2n/3 because D ≤ n and we may assume n ≥ 3". -/ +theorem le_sOf_and_sqrt_le (hD : 16 ≤ D) (hDn : D ≤ n) : + 3 * Real.sqrt D / 4 ≤ (sOf D : ℝ) ∧ Real.sqrt D ≤ 2 * (n : ℝ) / 3 := by + have hsqrtD := four_le_sqrt hD + have hsqrtn := four_le_sqrt (hD.trans hDn) + constructor + · calc 3 * Real.sqrt D / 4 ≤ Real.sqrt D - 1 := by linarith + _ ≤ (sOf D : ℝ) := by linarith [sqrt_lt_sOf_add_one D] + · calc Real.sqrt D ≤ Real.sqrt n := Real.sqrt_le_sqrt (by exact_mod_cast hDn) + _ ≤ 2 * (n : ℝ) / 3 := by + -- `n = √n √n ≥ 4√n` + linarith [mul_le_mul_of_nonneg_right hsqrtn (Real.sqrt_nonneg n), + Real.mul_self_sqrt (Nat.cast_nonneg n : (0 : ℝ) ≤ n)] + +section Matrices + +variable (wAB wBC wAC : Fin n → Fin n → ℤ) + +/-! ### The sets `W_ϱ` -/ + +/-- A pair lies in `W_ϱ` exactly if `ϱ` is the residue of its weight. -/ +theorem mem_residueClass (hp : p ≠ 0) (ϱ : Fin p) (q : Fin n × Fin n) : + q ∈ residueClass wAB p ϱ ↔ ϱ = resFin hp (wAB q.1 q.2) := by + rw [← modEq_iff_eq_resFin] + simp only [residueClass, Finset.mem_filter, Finset.mem_univ, true_and] + exact ⟨Int.ModEq.symm, Int.ModEq.symm⟩ + +/-- Proof of Theorem 17, the sets `W_ϱ` for `ϱ ∈ ℤ_p`: every pair `(a,b) ∈ A × B` lies in exactly +one of them. -/ +theorem existsUnique_mem_residueClass (hp : p ≠ 0) (q : Fin n × Fin n) : + ∃! ϱ : Fin p, q ∈ residueClass wAB p ϱ := + ⟨resFin hp (wAB q.1 q.2), (mem_residueClass wAB hp _ q).mpr rfl, + fun ϱ hϱ => (mem_residueClass wAB hp ϱ q).mp hϱ⟩ + +/-- The sets `W_ϱ` have `n²` pairs in all. -/ +private theorem sum_card_residueClass (hp : p ≠ 0) : + ∑ ϱ : Fin p, (residueClass wAB p ϱ).card = n ^ 2 := by + calc ∑ ϱ : Fin p, (residueClass wAB p ϱ).card + = ∑ ϱ : Fin p, ((Finset.univ : Finset (Fin n × Fin n)).filter + fun q => resFin hp (wAB q.1 q.2) = ϱ).card := by + refine Finset.sum_congr rfl fun ϱ _ => congrArg Finset.card ?_ + ext q + rw [mem_residueClass wAB hp] + simp [eq_comm] + _ = (Finset.univ : Finset (Fin n × Fin n)).card := + (Finset.card_eq_sum_card_fiberwise fun _ _ => Finset.mem_coe.mpr (Finset.mem_univ _)).symm + _ = n ^ 2 := by simp [sq] + +/-! ### The chunks -/ + +/-- Proof of Theorem 17: "cut it into chunks of at most n²/√D query pairs". (The chunks of `W_ϱ` +cover it: `biUnion_chunk`.) -/ +theorem card_chunkOf_le (hD : 1 ≤ D) (hDn : D ≤ n) (ϱ : Fin p) (j : ℕ) : + ((chunkOf wAB D p ϱ j).card : ℝ) ≤ (n : ℝ) ^ 2 / Real.sqrt D := + calc ((chunkOf wAB D p ϱ j).card : ℝ) ≤ (queryCap n D : ℝ) := by + exact_mod_cast card_chunk_le _ (one_le_queryCap hD hDn) j + _ ≤ (n : ℝ) ^ 2 / Real.sqrt D := by + rw [queryCap_eq_floor] + exact Nat.floor_le (by positivity) + +/-- There are at most `p + s` chunks in all. Each `W_ϱ` has at most `(|W_ϱ| + cap − 1)/cap` chunks, +where `cap = ⌊n²/√D⌋`, so `p + s + 1` chunks or more would need `n² ≥ (s + 1) cap + p` pairs. -/ +private theorem totalChunks_le_add_sOf (hD : 16 ≤ D) (hDn : D ≤ n) (hp0 : p ≠ 0) + (hp : Real.sqrt D / 2 ≤ (p : ℝ)) : totalChunks wAB D p ≤ p + sOf D := by + have hsum : totalChunks wAB D p * queryCap n D + p ≤ n ^ 2 + p * queryCap n D := by + have h := Finset.sum_le_sum fun ϱ (_ : ϱ ∈ (Finset.univ : Finset (Fin p))) => + numChunks_mul_le (residueClass wAB p ϱ) (queryCap n D) (one_le_queryCap (by omega) hDn) + simp only [Finset.sum_add_distrib, ← Finset.sum_mul, Finset.sum_const, Finset.card_univ, + Fintype.card_fin, smul_eq_mul, mul_one] at h + rwa [sum_card_residueClass wAB hp0] at h + have hpairs := sq_lt_mul_queryCap_add hD hDn hp + by_contra hmany + have hmul : (p + sOf D + 1) * queryCap n D ≤ totalChunks wAB D p * queryCap n D := + Nat.mul_le_mul_right _ (by omega) + rw [show (p + sOf D + 1) * queryCap n D = p * queryCap n D + (sOf D + 1) * queryCap n D by + ring] at hmul + -- `hsum` and `hmul` give `(s + 1) cap + p ≤ n²`, against `hpairs`. + omega + +/-- Proof of Theorem 17: "There are at most p + √D ≤ 2√D chunks in all." + +NOTE. A chunk holds a whole number of pairs, at most `⌊n²/√D⌋`, so the count is at most +`p + (n² − p)/⌊n²/√D⌋`, and the second term can exceed `√D` slightly. The printed bound is still +true, because the count is an integer, `D` is an integer and `D ≤ n` (`sq_lt_mul_queryCap_add`). -/ +theorem totalChunks_le (hD : 16 ≤ D) (hDn : D ≤ n) (hp : p ∈ primesInRange D) : + (totalChunks wAB D p : ℝ) ≤ (p : ℝ) + Real.sqrt D ∧ + (p : ℝ) + Real.sqrt D ≤ 2 * Real.sqrt D := by + obtain ⟨hprime, hge, hlt⟩ := mem_primesInRange.mp hp + have hnat : (totalChunks wAB D p : ℝ) ≤ (p : ℝ) + (sOf D : ℝ) := by + exact_mod_cast totalChunks_le_add_sOf wAB hD hDn hprime.ne_zero hge + exact ⟨by linarith [sOf_le_sqrt D], by linarith⟩ + +/-! ### The instances -/ + +/-- Proof of Theorem 17: "middle part C_k × ℤ_p, of size at most sp ≤ D". So each instance is an +instance of Lop-AE-SparseTri(n, D). -/ +theorem card_middle_lopInstance (hD : 16 ≤ D) (hg1 : 1 ≤ g) (hp : p ∈ primesInRange D) + (ι : InstanceIndex p) : Fintype.card (lopInstance wAB wBC wAC D g p ι).M ≤ D := by + obtain ⟨-, -, hlt⟩ := mem_primesInRange.mp hp + have hsp : (sOf D : ℝ) * (p : ℝ) ≤ D := + calc (sOf D : ℝ) * (p : ℝ) ≤ Real.sqrt D * Real.sqrt D := + mul_le_mul (sOf_le_sqrt D) hlt.le (Nat.cast_nonneg p) (Real.sqrt_nonneg _) + _ = D := Real.mul_self_sqrt (Nat.cast_nonneg D) + calc Fintype.card ({c : Fin n // c ∈ piece n D g ι.2.2} × Fin p) + = (piece n D g ι.2.2).card * p := by + rw [Fintype.card_prod, Fintype.card_coe, Fintype.card_fin] + _ ≤ sOf D * p := Nat.mul_le_mul_right p + ((card_piece_le hD hg1 n ι.2.2).trans (pieceSize_le_sOf hg1)) + _ ≤ D := by exact_mod_cast hsp + +/-- A query pair of the instance with index `(ϱ, j, k)` lies in `W_ϱ`, and `j` is the number of its +chunk. -/ +theorem mem_lopInstance_W {wAB wBC wAC : Fin n → Fin n → ℤ} + {ι : InstanceIndex p} {q : Fin n × Fin n} + (hq : q ∈ (lopInstance wAB wBC wAC D g p ι).W) : + q ∈ residueClass wAB p ι.1 ∧ rankIn (residueClass wAB p ι.1) q / queryCap n D = ι.2.1 := + Finset.mem_filter.mp hq + +/-- Proof of Theorem 17: "Within the chunk, the condition S(a,b,c) ≡ 0 (mod p) has become the +equality w(a,c) + ϱ ≡ −w(b,c) of a label of (a,c) and a label of (b,c)". -/ +theorem weightSum_modEq_zero_iff {p ϱ : ℤ} {a b : Fin n} (hϱ : wAB a b ≡ ϱ [ZMOD p]) (c : Fin n) : + weightSum wAB wBC wAC a b c ≡ 0 [ZMOD p] ↔ wAC a c + ϱ ≡ -wBC b c [ZMOD p] := by + rw [Int.modEq_iff_dvd, Int.modEq_iff_dvd, + show 0 - weightSum wAB wBC wAC a b c = -wBC b c - (wAC a c + ϱ) + (ϱ - wAB a b) by + simp only [weightSum]; ring] + exact dvd_add_left (Int.modEq_iff_dvd.mp hϱ) + +/-- Proof of Theorem 17: "a query pair (a,b) ∈ 𝒬 has a common neighbor if and only if some c ∈ C_k +has S(a,b,c) ≡ 0 (mod p)". -/ +theorem hasCommonNeighbor_lopInstance_iff (hp : p ≠ 0) + {ι : InstanceIndex p} {q : Fin n × Fin n} + (hq : q ∈ (lopInstance wAB wBC wAC D g p ι).W) : + (lopInstance wAB wBC wAC D g p ι).HasCommonNeighbor q.1 q.2 ↔ + ∃ c ∈ piece n D g ι.2.2, weightSum wAB wBC wAC q.1 q.2 c ≡ 0 [ZMOD (p : ℤ)] := by + have hϱ : wAB q.1 q.2 ≡ ((ι.1 : ℕ) : ℤ) [ZMOD (p : ℤ)] := + (Finset.mem_filter.mp (mem_lopInstance_W hq).1).2 + constructor + · -- A common neighbor `(c, σ)` carries the label `σ` of `(a,c)` and of `(b,c)`. + rintro ⟨⟨c, σ⟩, hA, hB⟩ + exact ⟨c, c.2, (weightSum_modEq_zero_iff wAB wBC wAC hϱ c).mpr (hA.symm.trans hB)⟩ + · -- If the two labels agree, the vertex `c` with that label is a common neighbor. + rintro ⟨c, hc, hS⟩ + have hB := resFin_modEq hp (-wBC q.2 c) + exact ⟨(⟨c, hc⟩, resFin hp (-wBC q.2 c)), + hB.trans ((weightSum_modEq_zero_iff wAB wBC wAC hϱ c).mp hS).symm, hB⟩ + +/-- Proof of Theorem 17, "For each chunk 𝒬 ⊆ W_ϱ and each piece C_k": `(ϱ, j, k)` is the index of an +instance exactly if `j` is the number of a chunk of `W_ϱ` and `k` the number of a piece. -/ +theorem mem_instanceIndices (ι : InstanceIndex p) : + ι ∈ instanceIndices wAB D g p ↔ + ι.2.1 < numChunks (residueClass wAB p ι.1) (queryCap n D) ∧ ι.2.2 < numPieces n D g := by + obtain ⟨ϱ, j, k⟩ := ι + simp only [instanceIndices, Finset.mem_biUnion, Finset.mem_univ, true_and, + Finset.mem_image, Finset.mem_product, Finset.mem_range, Prod.mk.injEq, Prod.exists] + exact ⟨fun ⟨_, _, _, h, h1, h2, h3⟩ => by subst h1 h2 h3; exact h, + fun h => ⟨ϱ, j, k, h, rfl, rfl, rfl⟩⟩ + +/-- There are as many instances as chunks times pieces. -/ +theorem card_instanceIndices (D g p : ℕ) : + (instanceIndices wAB D g p).card = totalChunks wAB D p * numPieces n D g := by + unfold instanceIndices totalChunks + rw [Finset.card_biUnion, Finset.sum_mul] + · refine Finset.sum_congr rfl fun ϱ _ => ?_ + rw [Finset.card_image_of_injective _ (Prod.mk_right_injective ϱ), Finset.card_product, + Finset.card_range, Finset.card_range] + · intro ϱ _ ϱ' _ hne + rw [Function.onFun, Finset.disjoint_left] + intro x hx hx' + obtain ⟨_, _, rfl⟩ := Finset.mem_image.mp hx + obtain ⟨_, _, h⟩ := Finset.mem_image.mp hx' + exact hne (congrArg Prod.fst h).symm + +/-- Proof of Theorem 17: "There are at most 2√D h ≤ 2√D(ng/s + 1) ≤ 4ng instances". -/ +theorem card_instanceIndices_le (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) + (hp : p ∈ primesInRange D) : + ((instanceIndices wAB D g p).card : ℝ) ≤ 4 * (n : ℝ) * (g : ℝ) := by + obtain ⟨hs, hsqrt⟩ := le_sOf_and_sqrt_le hD hDn + obtain ⟨hchunks, hchunks'⟩ := totalChunks_le wAB hD hDn hp + have hs0 := sOf_pos hD + have hpieces : (numPieces n D g : ℝ) ≤ (n : ℝ) * (g : ℝ) / (sOf D : ℝ) + 1 := + calc (numPieces n D g : ℝ) ≤ (⌈(n : ℝ) * (g : ℝ) / (sOf D : ℝ)⌉₊ : ℝ) := by + exact_mod_cast numPieces_le hD hg1 + _ ≤ (n : ℝ) * (g : ℝ) / (sOf D : ℝ) + 1 := (Nat.ceil_lt_add_one (by positivity)).le + -- `2√D · ng/s ≤ 8ng/3` because `s ≥ 3√D/4`, and `2√D ≤ 4n/3 ≤ 4ng/3`. + have hmain : + Real.sqrt D * ((n : ℝ) * (g : ℝ) / (sOf D : ℝ)) ≤ 4 / 3 * ((n : ℝ) * (g : ℝ)) := by + rw [← mul_div_assoc, div_le_iff₀ hs0] + calc Real.sqrt D * ((n : ℝ) * (g : ℝ)) + = 4 / 3 * ((n : ℝ) * (g : ℝ)) * (3 * Real.sqrt D / 4) := by ring + _ ≤ 4 / 3 * ((n : ℝ) * (g : ℝ)) * (sOf D : ℝ) := by gcongr + have hng : (n : ℝ) ≤ (n : ℝ) * (g : ℝ) := + le_mul_of_one_le_right (Nat.cast_nonneg n) (by exact_mod_cast hg1) + calc ((instanceIndices wAB D g p).card : ℝ) + = (totalChunks wAB D p : ℝ) * (numPieces n D g : ℝ) := by + rw [card_instanceIndices, Nat.cast_mul] + _ ≤ 2 * Real.sqrt D * (numPieces n D g : ℝ) := by gcongr; exact hchunks.trans hchunks' + _ ≤ 2 * Real.sqrt D * ((n : ℝ) * (g : ℝ) / (sOf D : ℝ) + 1) := by gcongr + _ ≤ 4 * (n : ℝ) * (g : ℝ) := by linarith [hmain, hng, hsqrt] + +end Matrices + +/-! ### The neighborhoods in an instance -/ + +/-- The graph of a function into `Fin p` has as many points as the domain. -/ +private theorem ncard_graph {α : Type} [Fintype α] (f : α → Fin p) : + {v : α × Fin p | v.2 = f v.1}.ncard = Fintype.card α := by + have h : {v : α × Fin p | v.2 = f v.1} = Set.univ.graphOn f := by + ext v + simp [Set.mem_graphOn, eq_comm] + rw [h, Set.ncard_graphOn, Set.ncard_univ, Nat.card_eq_fintype_card] + +/-- Section 3.1, "The reductions below [...] produce at most n²/√D query pairs and sets of at most +√D elements", and Remark 20, "every vertex of A has O(√D/g) neighbors in the middle part": in an +instance of the reduction, a vertex of `A` or of `B` has exactly one neighbor `(c, σ)` for every +`c ∈ C_k`, hence at most `⌈s/g⌉` neighbors. (And `⌈s/g⌉ ≤ s ≤ √D`; for `g ≤ √D` also +`⌈s/g⌉ ≤ 2√D/g`, `pieceSize_le_two_mul_sqrt_div`.) -/ +theorem ncard_nbr_lopInstance_le (wAB wBC wAC : Fin n → Fin n → ℤ) (hD : 16 ≤ D) + (hg1 : 1 ≤ g) (hp : p ≠ 0) (ι : InstanceIndex p) : + (∀ a, {v | (lopInstance wAB wBC wAC D g p ι).adjA a v}.ncard ≤ pieceSize D g) ∧ + (∀ b, {v | (lopInstance wAB wBC wAC D g p ι).adjB v b}.ncard ≤ pieceSize D g) := by + -- The neighbors of a vertex are the graph of the function that sends `c` to its label. + have hgraph : ∀ (f : {c : Fin n // c ∈ piece n D g ι.2.2} → ℤ), + {v : {c : Fin n // c ∈ piece n D g ι.2.2} × Fin p | + ((v.2 : ℕ) : ℤ) ≡ f v.1 [ZMOD (p : ℤ)]}.ncard ≤ pieceSize D g := by + intro f + refine le_trans (le_of_eq ?_) (card_piece_le hD hg1 n ι.2.2) + rw [← Fintype.card_coe, ← ncard_graph fun c => resFin hp (f c)] + congr 1 + ext v + exact modEq_iff_eq_resFin hp v.2 _ + exact ⟨fun a => hgraph fun c => wAC a c + ((ι.1 : ℕ) : ℤ), fun b => hgraph fun c => -wBC b c⟩ + +end Cslib.Algorithms.ExactTriangle diff --git a/Complexity/Algorithms/ExactTriangle/Witnesses.lean b/Complexity/Algorithms/ExactTriangle/Witnesses.lean new file mode 100644 index 0000000000..dfaeebb6aa --- /dev/null +++ b/Complexity/Algorithms/ExactTriangle/Witnesses.lean @@ -0,0 +1,402 @@ +/- +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.ExactTriangle.Splitting +public import Mathlib.Combinatorics.Enumerative.DoubleCounting +import Mathlib.Analysis.Complex.ExponentialBounds + +/-! +# Finite witness recovery from modular lopsided answers + +Correct finite answers are ordinary `Set.EqOn` premises on the actual lopsided graphs. +The recovery scans each accepted piece until it finds an exact zero. Distinct failed scans +contain distinct modular false positives. No solver, native program or time oracle is assumed. + +## Main definitions + +* `scanPiece`, `runScans`: finite witness search and number of scans. +* `AnswersCorrect`, `IsScanOrder`: correctness on queried pairs and a duplicate-free scan order. + +## Main results + +* `recovery_correct`: sound and complete recovery with at most `falsePositiveCount + 1` scans. +* `card_failedScans_le`: false-positive accounting. +* `falsePositiveCount_add_one_mul_pieceSize_le`: bound on the number of vertices scanned. + +## References + +Adapted from [Witnesses.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Witnesses.lean) and the minimal scan definitions in +`PaperStatements.lean` at `anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +Context: [AlmanWilliams2026, Section 3]. Quoted docstring fragments reproduce comments in the pinned +formal source; other prose describes the local Lean API. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.ExactTriangle + +open Cslib.Matrix + +/-- A finite search result returns a zero triangle, or certifies that no zero triangle exists. +This specifies mathematical correctness only. -/ +def IsSearchAnswer {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) + (r : Option (Fin n × Fin n × Fin n)) : Prop := + (∀ a b c, r = some (a, b, c) → weightSum wAB wBC wAC a b c = 0) ∧ + (r = none → ¬ HasZeroTriangle wAB wBC wAC) + +/-- Proof of Theorem 17: "scan the piece C_k of its instance for a c with S(a,b,c) = 0". The +vertices of the piece are tried in increasing order; the result is the first `c` found, or `none` if +the scan fails. -/ +def scanPiece {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D g k : ℕ) (a b : Fin n) : + Option (Fin n) := + ((List.finRange n).filter fun c => c ∈ piece n D g k).find? fun c => + weightSum wAB wBC wAC a b c = 0 + +/-- Proof of Theorem 17: "For every query pair that the oracle accepts". `ans ι` is the oracle's +answer to the instance with index `ι`. The result is the set of all (instance, query pair of that +instance) that the oracle accepts; each of them calls for one scan. -/ +def acceptedPairs {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (ans : InstanceIndex p → Fin n × Fin n → Bool) : Finset (InstanceIndex p × (Fin n × Fin n)) := + ((instanceIndices wAB D g p) ×ˢ Finset.univ).filter fun (ι, q) => + q ∈ (lopInstance wAB wBC wAC D g p ι).W ∧ ans ι q = true + +/-- Proof of Theorem 17: the scans, one after the other, over a list of accepted pairs: "We stop as +soon as a zero triangle is found." The result is the triangle found (or `none`) together with the +number of scans that were carried out. -/ +def runScans {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D g : ℕ) {p : ℕ} : + List (InstanceIndex p × (Fin n × Fin n)) → Option (Fin n × Fin n × Fin n) × ℕ + | [] => (none, 0) + | ((_, _, k), a, b) :: rest => + match scanPiece wAB wBC wAC D g k a b with + | some c => (some (a, b, c), 1) + | none => ((runScans wAB wBC wAC D g rest).1, (runScans wAB wBC wAC D g rest).2 + 1) + +/-- Proof of Theorem 17: "For every query pair that the oracle accepts, scan [...]". The paper does +not say in which order. `L` is an order of the scans: a list in which every accepted pair occurs +exactly once. -/ +def IsScanOrder {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (ans : InstanceIndex p → Fin n × Fin n → Bool) (L : List (InstanceIndex p × (Fin n × Fin n))) : + Prop := + L.Nodup ∧ ∀ x, x ∈ L ↔ x ∈ acceptedPairs wAB wBC wAC D g p ans + + +variable {n D g p : ℕ} + +section Matrices + +variable {wAB wBC wAC : Fin n → Fin n → ℤ} + +/-! ### Scans and accepted pairs -/ + +/-- Proof of Theorem 17, "scan the piece C_k of its instance for a c with S(a,b,c) = 0": a `c` that +a scan returns lies in the piece and has `S(a,b,c) = 0`. -/ +theorem scanPiece_eq_some {k : ℕ} {a b c : Fin n} (h : scanPiece wAB wBC wAC D g k a b = some c) : + c ∈ piece n D g k ∧ weightSum wAB wBC wAC a b c = 0 := by + have hmem := List.mem_of_find?_eq_some h + rw [List.mem_filter] at hmem + exact ⟨by simpa using hmem.2, by simpa using List.find?_some h⟩ + +/-- A scan fails exactly if the piece has no `c` with `S(a,b,c) = 0`. -/ +theorem scanPiece_eq_none_iff {k : ℕ} {a b : Fin n} : + scanPiece wAB wBC wAC D g k a b = none ↔ + ∀ c ∈ piece n D g k, weightSum wAB wBC wAC a b c ≠ 0 := by + unfold scanPiece + rw [List.find?_eq_none] + simp [List.mem_filter] + +variable (wAB wBC wAC) (D g) + +/-- The scan that belongs to `x = ((ϱ, j, k), (a, b))`, a query pair `(a, b)` of the instance with +index `(ϱ, j, k)`: the piece `C_k` is scanned for the pair `(a, b)`. -/ +abbrev scanOf (wAB wBC wAC : Fin n → Fin n → ℤ) (D g : ℕ) + (x : InstanceIndex p × (Fin n × Fin n)) : Option (Fin n) := + scanPiece wAB wBC wAC D g x.1.2.2 x.2.1 x.2.2 + +/-- What the scans over a list of pairs return: a triangle that is returned is a zero triangle; if +none is returned, then every scan failed; and all the scans that are carried out, except possibly +the last one, fail. -/ +private theorem runScans_spec (wAB wBC wAC : Fin n → Fin n → ℤ) (D g : ℕ) + (l : List (InstanceIndex p × (Fin n × Fin n))) : + (∀ a b c, (runScans wAB wBC wAC D g l).1 = some (a, b, c) → weightSum wAB wBC wAC a b c = 0) ∧ + ((runScans wAB wBC wAC D g l).1 = none → ∀ x ∈ l, scanOf wAB wBC wAC D g x = none) ∧ + (runScans wAB wBC wAC D g l).2 ≤ + (l.filter fun x => decide (scanOf wAB wBC wAC D g x = none)).length + 1 := by + induction l with + | nil => simp [runScans] + | cons x rest ih => + obtain ⟨ihzero, ihnone, ihcount⟩ := ih + cases hscan : scanOf wAB wBC wAC D g x with + | some c => + -- The first scan finds `c`, and the scans stop. + have hrun : runScans wAB wBC wAC D g (x :: rest) = (some (x.2.1, x.2.2, c), 1) := by + simp only [runScans, hscan] + rw [hrun] + refine ⟨fun a b c' h => ?_, fun h => by simp at h, by simp⟩ + simp only [Option.some.injEq, Prod.mk.injEq] at h + obtain ⟨rfl, rfl, rfl⟩ := h + exact (scanPiece_eq_some hscan).2 + | none => + -- The first scan fails, and the scans go on with the rest of the list. + have hrun : runScans wAB wBC wAC D g (x :: rest) + = ((runScans wAB wBC wAC D g rest).1, (runScans wAB wBC wAC D g rest).2 + 1) := by + simp only [runScans, hscan] + rw [hrun, List.filter_cons_of_pos (by simpa using hscan), List.length_cons] + refine ⟨ihzero, fun h y hy => ?_, by omega⟩ + rcases List.mem_cons.mp hy with rfl | hy + · exact hscan + · exact ihnone h y hy + +variable (p) (ans : InstanceIndex p → Fin n × Fin n → Bool) + +/-- The hypothesis on the oracle: `ans ι` is a correct answer to the instance with index `ι`, for +every instance of the reduction. -/ +def AnswersCorrect (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (ans : InstanceIndex p → Fin n × Fin n → Bool) : Prop := + ∀ ι ∈ instanceIndices wAB D g p, Set.EqOn (fun q => ans ι q = true) + (fun q => (lopInstance wAB wBC wAC D g p ι).HasCommonNeighbor q.1 q.2) + ((lopInstance wAB wBC wAC D g p ι).W : Set _) + +/-- Proof of Theorem 17: a "failed scan, of a piece C_k for a pair (a,b)": the oracle accepted the +pair in the instance, and the piece contains no `c` with `S(a,b,c) = 0`. -/ +def failedScans (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (ans : InstanceIndex p → Fin n × Fin n → Bool) : Finset (InstanceIndex p × (Fin n × Fin n)) := + (acceptedPairs wAB wBC wAC D g p ans).filter fun x => scanOf wAB wBC wAC D g x = none + +variable {wAB wBC wAC D g p ans} + +/-- The accepted pairs: an index of an instance, a query pair of that instance, and the answer +"yes". -/ +theorem mem_acceptedPairs {ι : InstanceIndex p} {q : Fin n × Fin n} : + (ι, q) ∈ acceptedPairs wAB wBC wAC D g p ans ↔ + ι ∈ instanceIndices wAB D g p ∧ q ∈ (lopInstance wAB wBC wAC D g p ι).W ∧ ans ι q = true := by + simp [acceptedPairs] + +/-! ### The sentences of the paragraph "Witnesses" -/ + +/-- Proof of Theorem 17: "A zero triangle is always found, since its c lies in some piece and makes +the oracle accept its pair." -/ +theorem exists_mem_acceptedPairs (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) (hp : p ≠ 0) + (hans : AnswersCorrect wAB wBC wAC D g p ans) {a b c : Fin n} + (h0 : weightSum wAB wBC wAC a b c = 0) : + ∃ ϱ j k, ((ϱ, j, k), (a, b)) ∈ acceptedPairs wAB wBC wAC D g p ans ∧ c ∈ piece n D g k := by + -- The piece of `c`, the residue class of `(a, b)` and the chunk of `(a, b)` in it. + obtain ⟨k, ⟨hk, hck⟩, -⟩ := existsUnique_mem_piece hD hg1 c + obtain ⟨ϱ, hab, -⟩ := existsUnique_mem_residueClass wAB hp (a, b) + obtain ⟨j, hj, hW⟩ := Finset.mem_biUnion.mp + ((biUnion_chunk (residueClass wAB p ϱ) (one_le_queryCap (by omega) hDn)).symm ▸ hab) + have hι : ((ϱ, j, k) : InstanceIndex p) ∈ instanceIndices wAB D g p := + (mem_instanceIndices wAB _).mpr ⟨Finset.mem_range.mp hj, hk⟩ + have hquery : (a, b) ∈ (lopInstance wAB wBC wAC D g p (ϱ, j, k)).W := hW + refine ⟨ϱ, j, k, mem_acceptedPairs.mpr ⟨hι, hquery, ?_⟩, hck⟩ + -- The vertex `c` makes the oracle accept the pair. + apply (Iff.of_eq (hans _ hι hquery)).mpr + rw [hasCommonNeighbor_lopInstance_iff wAB wBC wAC hp hquery] + exact ⟨c, hck, by rw [show weightSum wAB wBC wAC a b c = 0 from h0]⟩ + +/-- Proof of Theorem 17: "A failed scan, of a piece C_k for a pair (a,b), contains a c ∈ C_k with +S(a,b,c) ≡ 0 (mod p) but S(a,b,c) ≠ 0, that is, a false positive (a,b,c) of p". -/ +theorem exists_isFalsePositive_of_mem_failedScans (hp : p ≠ 0) + (hans : AnswersCorrect wAB wBC wAC D g p ans) + {ϱ : Fin p} {j k : ℕ} {a b : Fin n} + (hx : ((ϱ, j, k), (a, b)) ∈ failedScans wAB wBC wAC D g p ans) : + ∃ c ∈ piece n D g k, IsFalsePositive wAB wBC wAC p (a, b, c) := by + obtain ⟨hacc, hnone⟩ := Finset.mem_filter.mp hx + obtain ⟨hι, hquery, htrue⟩ := mem_acceptedPairs.mp hacc + obtain ⟨c, hc, hS⟩ := (hasCommonNeighbor_lopInstance_iff wAB wBC wAC hp hquery).mp + ((Iff.of_eq (hans _ hι hquery)).mp htrue) + exact ⟨c, hc, scanPiece_eq_none_iff.mp hnone c hc, Int.modEq_zero_iff_dvd.mp hS⟩ + +/-- Proof of Theorem 17: "distinct scans contain distinct false positives": a triple `(a,b,c)` +belongs to the scan of at most one accepted pair. -/ +theorem eq_of_mem_acceptedPairs {ϱ ϱ' : Fin p} {j j' k k' : ℕ} {q : Fin n × Fin n} + (hx : ((ϱ, j, k), q) ∈ acceptedPairs wAB wBC wAC D g p ans) + (hy : ((ϱ', j', k'), q) ∈ acceptedPairs wAB wBC wAC D g p ans) {c : Fin n} + (hc : c ∈ piece n D g k) + (hc' : c ∈ piece n D g k') : (ϱ, j, k) = (ϱ', j', k') := by + have hp : p ≠ 0 := by + rintro rfl + exact ϱ.elim0 + obtain ⟨hclass, hchunk⟩ := mem_lopInstance_W (mem_acceptedPairs.mp hx).2.1 + obtain ⟨hclass', hchunk'⟩ := mem_lopInstance_W (mem_acceptedPairs.mp hy).2.1 + -- The pair determines its residue class, then its chunk; the vertex `c` determines the piece. + obtain rfl : ϱ = ϱ' := (existsUnique_mem_residueClass wAB hp q).unique hclass hclass' + obtain rfl : j = j' := hchunk.symm.trans hchunk' + obtain rfl : k = k' := ((mem_piece c).mp hc).symm.trans ((mem_piece c).mp hc') + rfl + +/-- Proof of Theorem 17: "distinct scans contain distinct false positives, so there are at most F(p) +failed scans", for any finite family `s` of scans and any notion "the scan `x` contains the triple +`t`". This applies to the accepted-pair scans in `card_failedScans_le`. -/ +theorem card_le_falsePositiveCount_of_distinct_scans + (wAB wBC wAC : Fin n → Fin n → ℤ) {σ : Type*} (s : Finset σ) + (contains : σ → Fin n × Fin n × Fin n → Prop) + (hfalse : ∀ x ∈ s, ∃ t, contains x t ∧ IsFalsePositive wAB wBC wAC p t) + (hdistinct : ∀ t, ∀ x ∈ s, ∀ y ∈ s, contains x t → contains y t → x = y) : + s.card ≤ falsePositiveCount wAB wBC wAC p := by + classical + refine Finset.card_le_card_of_forall_subsingleton contains (fun x hx => ?_) + fun t _ x hx y hy => hdistinct t x hx.1 y hy.1 hx.2 hy.2 + obtain ⟨t, ht, hfp⟩ := hfalse x hx + exact ⟨t, Finset.mem_filter.mpr ⟨Finset.mem_univ _, hfp⟩, ht⟩ + +/-- Proof of Theorem 17: "so there are at most F(p) failed scans". -/ +theorem card_failedScans_le (hp : p ≠ 0) + (hans : AnswersCorrect wAB wBC wAC D g p ans) : + (failedScans wAB wBC wAC D g p ans).card ≤ falsePositiveCount wAB wBC wAC p := by + -- The scan of the piece `C_k` for the pair `(a, b)` contains the triples `(a, b, c)`, `c ∈ C_k`. + refine card_le_falsePositiveCount_of_distinct_scans wAB wBC wAC _ + (fun x t => x.2 = (t.1, t.2.1) ∧ t.2.2 ∈ piece n D g x.1.2.2) ?_ ?_ + · rintro ⟨⟨ϱ, j, k⟩, a, b⟩ hx + obtain ⟨c, hc, hfalse⟩ := exists_isFalsePositive_of_mem_failedScans hp hans hx + exact ⟨(a, b, c), ⟨rfl, hc⟩, hfalse⟩ + · rintro ⟨a, b, c⟩ ⟨⟨ϱ, j, k⟩, q⟩ hx ⟨⟨ϱ', j', k'⟩, q'⟩ hy ⟨rfl, hc⟩ ⟨rfl, hc'⟩ + rw [eq_of_mem_acceptedPairs (Finset.mem_filter.mp hx).1 (Finset.mem_filter.mp hy).1 hc hc'] + +end Matrices + +/-- Proof of Theorem 17: the length of a scan, the `O(√D/g)` in "the scans cost O((F(p) + 1)√D/g)": +a piece has at most `⌈s/g⌉ ≤ 2√D/g` vertices. -/ +theorem pieceSize_le_two_mul_sqrt_div (hg1 : 1 ≤ g) (hg : (g : ℝ) ≤ Real.sqrt D) : + (pieceSize D g : ℝ) ≤ 2 * Real.sqrt D / (g : ℝ) := by + have hg0 : (0 : ℝ) < g := by exact_mod_cast hg1 + calc (pieceSize D g : ℝ) ≤ (sOf D : ℝ) / (g : ℝ) + 1 := by + rw [pieceSize_eq_ceil] + exact (Nat.ceil_lt_add_one (by positivity)).le + _ ≤ Real.sqrt D / (g : ℝ) + Real.sqrt D / (g : ℝ) := + add_le_add (div_le_div_of_nonneg_right (sOf_le_sqrt D) hg0.le) + ((one_le_div hg0).mpr hg) + _ = 2 * Real.sqrt D / (g : ℝ) := by ring + +/-- For `16 ≤ D ≤ n`, `1 ≤ g ≤ sqrt D`, `κ ≥ 1` and edge weights bounded by `n^κ`, bound +`(falsePositiveCount + 1) * pieceSize` at any selected prime by +`(2 * Hashing.falsePositiveConst + 2) * κ * n^3 * log n / g`. +This arithmetic estimate counts candidate vertices rather than native instructions. -/ +theorem falsePositiveCount_add_one_mul_pieceSize_le {κ : ℝ} (wAB wBC wAC : Fin n → Fin n → ℤ) + (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) (hg : (g : ℝ) ≤ Real.sqrt D) (hκ : 1 ≤ κ) + (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) (hp : IsSelectedPrime wAB wBC wAC D p) : + ((falsePositiveCount wAB wBC wAC p : ℝ) + 1) * (pieceSize D g : ℝ) + ≤ (2 * Hashing.falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by + have hsqrt : 4 ≤ Real.sqrt (D : ℝ) := + Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) + have hn : (16 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn + have hlog : 1 ≤ Real.log n := + (Real.le_log_iff_exp_le (by linarith)).2 (Real.exp_one_lt_three.le.trans (by linarith)) + set B := κ * (n : ℝ) ^ 3 * Real.log n with hB + -- The two terms: `F(p) √D = O(B)` by the bound on `F(p)`, and `√D ≤ n ≤ n³ ≤ B`. + have hfalse : (falsePositiveCount wAB wBC wAC p : ℝ) * Real.sqrt D ≤ + Hashing.falsePositiveConst * B := by + have hF := Hashing.falsePositiveCount_le_const_mul hD hDn hκ wAB wBC wAC hT hp + rwa [← mul_div_assoc, le_div_iff₀ (by linarith)] at hF + have hone : Real.sqrt D ≤ B := + calc Real.sqrt D ≤ n := sqrt_le_natCast (by omega) hDn + _ ≤ (n : ℝ) ^ 3 := le_self_pow₀ (by linarith) (by norm_num) + _ ≤ κ * (n : ℝ) ^ 3 := le_mul_of_one_le_left (by positivity) hκ + _ ≤ B := le_mul_of_one_le_right (by positivity) hlog + calc ((falsePositiveCount wAB wBC wAC p : ℝ) + 1) * (pieceSize D g : ℝ) + ≤ ((falsePositiveCount wAB wBC wAC p : ℝ) + 1) * (2 * Real.sqrt D / (g : ℝ)) := by + gcongr + exact pieceSize_le_two_mul_sqrt_div hg1 hg + _ = (2 * ((falsePositiveCount wAB wBC wAC p : ℝ) * Real.sqrt D) + + 2 * Real.sqrt D) / (g : ℝ) := by ring + _ ≤ (2 * (Hashing.falsePositiveConst * B) + 2 * B) / (g : ℝ) := by gcongr + _ = (2 * Hashing.falsePositiveConst + 2) * (B / (g : ℝ)) := by ring + +/-! ### The output of the reduction -/ + +/-- **Theorem 17, correctness.** Given any correct answers of an oracle for +Lop-AE-SparseTri to the instances of the reduction, and whatever the order `L` of the scans, the +reduction's output is correct: it is a zero triangle if one exists, and `none` otherwise. It +carries out at most `F(p) + 1` scans ("We stop as soon as a zero triangle is found"). -/ +theorem recovery_correct (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) + (wAB wBC wAC : Fin n → Fin n → ℤ) (hp : IsSelectedPrime wAB wBC wAC D p) + {ans : InstanceIndex p → Fin n × Fin n → Bool} (hans : AnswersCorrect wAB wBC wAC D g p ans) + {L : List (InstanceIndex p × (Fin n × Fin n))} (hL : IsScanOrder wAB wBC wAC D g p ans L) : + IsSearchAnswer wAB wBC wAC (runScans wAB wBC wAC D g L).1 ∧ + (runScans wAB wBC wAC D g L).2 ≤ falsePositiveCount wAB wBC wAC p + 1 := by + have hp0 : p ≠ 0 := (mem_primesInRange.mp hp.1).1.ne_zero + obtain ⟨hzero, hnone, hcount⟩ := runScans_spec wAB wBC wAC D g L + refine ⟨⟨hzero, fun hout ⟨a, b, c, h0⟩ => ?_⟩, ?_⟩ + · -- "A zero triangle is always found": its scan cannot have failed. + obtain ⟨ϱ, j, k, hx, hc⟩ := exists_mem_acceptedPairs hD hDn hg1 hp0 hans h0 + exact scanPiece_eq_none_iff.mp (hnone hout _ ((hL.2 _).mpr hx)) c hc h0 + · -- All the scans except possibly the last one fail, and there are at most `F(p)` failed scans. + have hfailed : (L.filter fun x => decide (scanOf wAB wBC wAC D g x = none)).length + ≤ (failedScans wAB wBC wAC D g p ans).card := by + rw [← List.toFinset_card_of_nodup (hL.1.filter _)] + refine Finset.card_le_card fun x hx => ?_ + obtain ⟨hxL, hxnone⟩ := List.mem_filter.mp (List.mem_toFinset.mp hx) + exact Finset.mem_filter.mpr ⟨(hL.2 x).mp hxL, by simpa using hxnone⟩ + have hF := card_failedScans_le hp0 hans + omega + + +/-- Every finite set of accepted queries admits a duplicate-free scan order. -/ +theorem exists_isScanOrder (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (ans : InstanceIndex p → Fin n × Fin n → Bool) : + ∃ L, IsScanOrder wAB wBC wAC D g p ans L := + ⟨(acceptedPairs wAB wBC wAC D g p ans).toList, + Finset.nodup_toList _, fun _ => Finset.mem_toList⟩ + +/-- Correct finite answer tables exist. This is a finite specification, with no solver or +complexity claim. Values outside the queried edges remain unconstrained. -/ +theorem exists_answersCorrect (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) : + ∃ ans, AnswersCorrect wAB wBC wAC D g p ans := by + classical + refine ⟨fun ι q => decide ((lopInstance wAB wBC wAC D g p ι).HasCommonNeighbor q.1 q.2), ?_⟩ + intro ι _ q hq + simp + +/-- Correct finite common-neighbor counts supply the Boolean answers used by recovery. +This connects directly to the lopsided matrix/counting interface. -/ +theorem answersCorrect_of_counting (wAB wBC wAC : Fin n → Fin n → ℤ) (D g p : ℕ) + (counts : InstanceIndex p → Fin n × Fin n → ℕ) + (hcounts : ∀ ι ∈ instanceIndices wAB D g p, + Set.EqOn (counts ι) (fun q => (lopInstance wAB wBC wAC D g p ι).commonNeighborCount q.1 q.2) + ((lopInstance wAB wBC wAC D g p ι).W : Set _)) : + AnswersCorrect wAB wBC wAC D g p (fun ι q => decide (counts ι q ≠ 0)) := by + intro ι hι q hq + apply propext + simp only [decide_eq_true_eq] + exact Iff.of_eq ((lopInstance wAB wBC wAC D g p ι).detection_eqOn_of_counting + (counts ι) (hcounts ι hι) hq) + +/-- The recovery scan examines at most the uniform false-positive bound times the piece size. +This counts candidate vertices only; it is not a charged native running-time theorem. -/ +theorem recovery_vertexBound {κ : ℝ} (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) + (hg : (g : ℝ) ≤ Real.sqrt D) (hκ : 1 ≤ κ) + (wAB wBC wAC : Fin n → Fin n → ℤ) + (hw : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) + (hp : IsSelectedPrime wAB wBC wAC D p) + {ans : InstanceIndex p → Fin n × Fin n → Bool} (hans : AnswersCorrect wAB wBC wAC D g p ans) + {L : List (InstanceIndex p × (Fin n × Fin n))} (hL : IsScanOrder wAB wBC wAC D g p ans L) : + IsSearchAnswer wAB wBC wAC (runScans wAB wBC wAC D g L).1 ∧ + ((runScans wAB wBC wAC D g L).2 : ℝ) * pieceSize D g ≤ + (2 * Hashing.falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by + obtain ⟨hcorrect, hcount⟩ := recovery_correct hD hDn hg1 wAB wBC wAC hp hans hL + refine ⟨hcorrect, le_trans ?_ (falsePositiveCount_add_one_mul_pieceSize_le + wAB wBC wAC hD hDn hg1 hg hκ hw hp)⟩ + gcongr + exact_mod_cast hcount + +/-- Finite reduction equivalence: scanning all accepted queries returns `none` exactly when the +original three matrices have no zero triangle. Any duplicate-free scan order works. -/ +theorem recovery_eq_none_iff (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) + (wAB wBC wAC : Fin n → Fin n → ℤ) (hp : IsSelectedPrime wAB wBC wAC D p) + {ans : InstanceIndex p → Fin n × Fin n → Bool} (hans : AnswersCorrect wAB wBC wAC D g p ans) + {L : List (InstanceIndex p × (Fin n × Fin n))} (hL : IsScanOrder wAB wBC wAC D g p ans L) : + (runScans wAB wBC wAC D g L).1 = none ↔ ¬ HasZeroTriangle wAB wBC wAC := by + obtain ⟨hsound, hcomplete⟩ := (recovery_correct hD hDn hg1 wAB wBC wAC hp hans hL).1 + refine ⟨hcomplete, fun hn => ?_⟩ + cases h : (runScans wAB wBC wAC D g L).1 with + | none => rfl + | some t => + obtain ⟨a, b, c⟩ := t + exact False.elim (hn ⟨a, b, c, hsound a b c h⟩) + +end Cslib.Algorithms.ExactTriangle diff --git a/Complexity/Algorithms/WordRAM/Histogram.lean b/Complexity/Algorithms/WordRAM/Histogram.lean new file mode 100644 index 0000000000..da7669c784 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Histogram.lean @@ -0,0 +1,402 @@ +/- +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.Copy +public import Mathlib.Algebra.BigOperators.Group.Finset.Basic +public import Mathlib.Data.Finset.Prod +import Mathlib.Tactic.Linarith +import Mathlib.Tactic.NormNum + +/-! +# Charged word-RAM histograms + +A fixed program clears a histogram and inserts the keys of a disjoint input array. +An insertion into a bucket of size `c` creates `2c` ordered collisions. Every address +calculation, load, increment, store and loop test is charged. + +## Main definitions + +* `Histogram.count`: bucket sizes in an input prefix. +* `Histogram.pairs`: ordered collisions accumulated by prefix insertion. +* `Histogram.program`: clear the buckets, then insert the input keys. + +## Main results + +* `Histogram.exec`: histogram construction and collision counting with a linear instruction bound. + +## References + +* [T. M. Chan and Q. He, *Reducing 3SUM to Convolution-3SUM*][ChanHe2020] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Histogram + +open Cslib.WordRAM +open Finset + +/-- Number of occurrences of a key among the first `k` input positions. -/ +def count (keys : ℕ → ℕ) (k b : ℕ) : ℕ := #{i ∈ range k | keys i = b} + +@[simp] theorem count_zero (keys : ℕ → ℕ) (b : ℕ) : count keys 0 b = 0 := by + simp [count] + +/-- Inserting one key changes just its bucket. -/ +theorem count_succ (keys : ℕ → ℕ) (k b : ℕ) : + count keys (k + 1) b = count keys k b + if keys k = b then 1 else 0 := by + by_cases h : keys k = b <;> simp [count, range_add_one, filter_insert, h] + +/-- A bucket cannot contain more positions than the prefix. -/ +theorem count_le (keys : ℕ → ℕ) (k b : ℕ) : count keys k b ≤ k := + (card_filter_le _ _).trans_eq (card_range k) + +/-- Ordered collisions accumulated while inserting successive keys. -/ +def pairs (keys : ℕ → ℕ) : ℕ → ℕ + | 0 => 0 + | k + 1 => pairs keys k + 2 * count keys k (keys k) + +/-- The prefix collision counter is at most the square of the prefix length. -/ +theorem pairs_le_sq (keys : ℕ → ℕ) (k : ℕ) : pairs keys k ≤ k * k := by + induction k with + | zero => simp [pairs] + | succ k ih => + have := count_le keys k (keys k) + dsimp only [pairs] + nlinarith + +/-- The insertion counter counts exactly the ordered pairs of different positions +with the same key, including when equal keys occur at different input positions. -/ +theorem pairs_eq_card (keys : ℕ → ℕ) (k : ℕ) : + pairs keys k = #{p ∈ (range k).offDiag | keys p.1 = keys p.2} := by + induction k with + | zero => simp [pairs] + | succ k ih => + let S := range k + let A := {p ∈ S.offDiag | keys p.1 = keys p.2} + let B := {p ∈ {k} ×ˢ S | keys p.1 = keys p.2} + let C := {p ∈ S ×ˢ {k} | keys p.1 = keys p.2} + have hrow : B = {k} ×ˢ {i ∈ S | keys i = keys k} := by + ext p + simp only [B, mem_filter, mem_product, mem_singleton] + aesop + have hcol : C = {i ∈ S | keys i = keys k} ×ˢ {k} := by + ext p + simp only [C, mem_filter, mem_product, mem_singleton] + aesop + have hAB : Disjoint A B := by + rw [disjoint_left] + intro p hp hq + have hp' := (mem_offDiag.mp (mem_filter.mp hp).1).1 + have hq' := mem_singleton.mp (mem_product.mp (mem_filter.mp hq).1).1 + have : p.1 < k := mem_range.mp hp' + omega + have hAC : Disjoint (A ∪ B) C := by + rw [disjoint_left] + intro p hp hq + have hq' := mem_singleton.mp (mem_product.mp (mem_filter.mp hq).1).2 + rcases mem_union.mp hp with hp | hp + · have hp' := (mem_offDiag.mp (mem_filter.mp hp).1).2.1 + have : p.2 < k := mem_range.mp hp' + omega + · have hp' := (mem_product.mp (mem_filter.mp hp).1).2 + have : p.2 < k := mem_range.mp hp' + omega + have hsplit : {p ∈ (range (k + 1)).offDiag | keys p.1 = keys p.2} = A ∪ B ∪ C := by + rw [range_add_one, offDiag_insert (notMem_range_self : k ∉ range k), + filter_union, filter_union] + rw [hsplit, card_union_of_disjoint hAC, card_union_of_disjoint hAB, hrow, hcol, + card_product, card_product, card_singleton, one_mul, mul_one] + change pairs keys k + 2 * count keys k (keys k) = + #{p ∈ (range k).offDiag | keys p.1 = keys p.2} + count keys k (keys k) + + count keys k (keys k) + rw [ih] + omega + +/-- The original memory, with the histogram block replaced by prefix bucket counts. -/ +def memory (mem : ℕ → ℕ) (keys : ℕ → ℕ) (dst bins k : ℕ) : ℕ → ℕ := + fun a => if dst ≤ a ∧ a < dst + bins then count keys k (a - dst) else mem a + +/-- Reading a histogram bucket. -/ +theorem memory_read (mem : ℕ → ℕ) (keys : ℕ → ℕ) (dst bins k b : ℕ) + (hb : b < bins) : memory mem keys dst bins k (dst + b) = count keys k b := by + simp [memory, hb] + +/-- Histogram writes preserve all other cells. -/ +theorem memory_outside (mem : ℕ → ℕ) (keys : ℕ → ℕ) (dst bins k a : ℕ) + (ha : ¬(dst ≤ a ∧ a < dst + bins)) : memory mem keys dst bins k a = mem a := by + simp [memory, ha] + +/-- One insertion is a single histogram-cell update. -/ +theorem memory_succ (mem : ℕ → ℕ) (keys : ℕ → ℕ) (dst bins k : ℕ) + (hk : keys k < bins) : + memory mem keys dst bins (k + 1) = Function.update (memory mem keys dst bins k) + (dst + keys k) (count keys k (keys k) + 1) := by + funext a + by_cases ha : a = dst + keys k + · subst a + simp [memory_read, hk, count_succ] + · rw [Function.update_of_ne ha] + simp only [memory] + split_ifs with h + · have hne : keys k ≠ a - dst := by omega + simp [count_succ, hne] + · rfl + +/-- Input base. -/ +abbrev src : Fin 11 := 0 +/-- Histogram base. -/ +abbrev dst : Fin 11 := 1 +/-- Number of input keys. -/ +abbrev len : Fin 11 := 2 +/-- Number of buckets. -/ +abbrev bins : Fin 11 := 3 +/-- Prefix length. -/ +abbrev idx : Fin 11 := 4 +/-- Temporary address. -/ +abbrev addr : Fin 11 := 5 +/-- Current key. -/ +abbrev key : Fin 11 := 6 +/-- Current bucket count. -/ +abbrev bucketSize : Fin 11 := 7 +/-- Number of ordered collisions. -/ +abbrev total : Fin 11 := 8 +/-- Constant one. -/ +abbrev one : Fin 11 := 9 +/-- Constant zero. -/ +abbrev zero : Fin 11 := 10 + +/-- Clear one histogram cell and advance the cursor. -/ +def clearBody : Code .standard 11 := ram[ + addr ← dst + idx; *addr ← zero; idx ← idx + one ] + +/-- Insert one key and add twice its old bucket count. -/ +def body : Code .standard 11 := ram[ + addr ← src + idx; key ← *addr; addr ← dst + key; + bucketSize ← *addr; total ← total + bucketSize; total ← total + bucketSize; + bucketSize ← bucketSize + one; *addr ← bucketSize; idx ← idx + one ] + +/-- Clear the histogram, retaining the input and all four parameters. -/ +def clear : Code .standard 11 := + .seq (ram[ idx ← #0; one ← #1; zero ← #0 ]) (.loop (.lt idx bins) clearBody) + +/-- Insert all keys into a cleared histogram. -/ +def insert : Code .standard 11 := + .seq (ram[ idx ← #0; total ← #0 ]) (.loop (.lt idx len) body) + +/-- Fixed code, with all bases and lengths supplied in registers. -/ +def program : Code .standard 11 := .seq clear insert + +private theorem exec_clear {w : ℕ} (x : Data 11) (h1 : Fits w 1) + (hd : Fits w (x.regs dst + x.regs bins)) : + ∃ y, ExecWithin w clear x y (5 * x.regs bins + 5) ∧ + y.mem = Cslib.WordRAM.clear x.mem (x.regs dst) (x.regs bins) ∧ + (∀ i, i.val < 4 → y.regs i = x.regs i) ∧ y.regs one = 1 := by + let I (y : Data 11) : Prop := + (∀ i, i.val < 4 → y.regs i = x.regs i) ∧ y.regs one = 1 ∧ y.regs zero = 0 ∧ + y.regs idx ≤ x.regs bins ∧ y.mem = Cslib.WordRAM.clear x.mem (x.regs dst) (y.regs idx) + have hbody (y : Data 11) (hy : I y) (hc : (Cond.lt idx bins).Holds y) : + ∃ z, ExecWithin w clearBody y z 3 ∧ I z ∧ + x.regs bins - z.regs idx < x.regs bins - y.regs idx := by + obtain ⟨hr, ho, hz, hi, hm⟩ := hy + have hb := hr bins (by decide) + have hdst := hr dst (by decide) + have hlt : y.regs idx < x.regs bins := by simpa only [Cond.holds_lt, hb] using hc + have ha : Fits w (y.regs dst + y.regs idx) := hd.of_le (by rw [hdst]; omega) + have hn : Fits w (y.regs idx + 1) := hd.of_le (by omega) + refine ⟨?z, ?run, ?inv, ?decrease⟩ + case run => ram_eval [clearBody, ha, hn, ho] + case inv => + refine ⟨?_, ?_, ?_, ?_, ?_⟩ + · intro i hi + have hia : i ≠ addr := by intro h; subst i; norm_num [addr] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + simpa only [Data.set_regs, Data.write_regs, Function.update_of_ne hia, + Function.update_of_ne hii] using hr i hi + · simpa +decide [Function.update_apply] using ho + · simpa +decide [Function.update_apply] using hz + · simpa [ho] using Nat.succ_le_of_lt hlt + · simp only [Data.set_mem, Data.write_mem, Data.set_regs, Function.update_apply] + rw [hz, hm, hdst, ← Cslib.WordRAM.clear_succ] + rfl + case decrease => simp only [Data.set_regs, Function.update_self]; omega + let x₀ := ((x.set idx 0).set one 1).set zero 0 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_⟩ + · intro i hi + have hiz : i ≠ zero := by intro h; subst i; norm_num [zero] at hi + have hio : i ≠ one := by intro h; subst i; norm_num [one] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + simp only [x₀, Data.set_regs, Function.update_of_ne hiz, Function.update_of_ne hio, + Function.update_of_ne hii] + · simp +decide [x₀, one, zero, idx] + · simp [x₀] + · simp +decide [x₀, one, zero, idx] + · simp +decide [x₀, one, zero, idx, Cslib.WordRAM.clear_zero] + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => x.regs bins - y.regs idx) hbody hx₀ + obtain ⟨hr, ho, hz, hi, hm⟩ := hy + have hb := hr bins (by decide) + have hend : y.regs idx = x.regs bins := by + simp only [Cond.holds_lt, hb] at hc + omega + refine ⟨y, ?_, by simpa only [hend] using hm, hr, ho⟩ + have hinit : ExecWithin (B := .standard) w (ram[ idx ← #0; one ← #1; zero ← #0 ]) x x₀ 3 := by + ram_eval [x₀, h1] + convert hinit.seq he using 1 + · rfl + · simp +decide only [x₀, Data.set_regs, Function.update_apply, idx, one, zero, ↓reduceIte, + Nat.sub_zero] + omega + +private theorem exec_insert {w : ℕ} (x : Data 11) (keys : ℕ → ℕ) + (hs : Fits w (x.regs src + x.regs len)) + (hd : Fits w (x.regs dst + x.regs bins)) + (ht : Fits w (x.regs len * x.regs len + x.regs len + 1)) + (hdisj : x.regs src + x.regs len ≤ x.regs dst) + (hkeys : ∀ i < x.regs len, keys i < x.regs bins) + (hinput : ∀ i < x.regs len, x.mem (x.regs src + i) = keys i) + (hzero : ∀ b < x.regs bins, x.mem (x.regs dst + b) = 0) + (ho : x.regs one = 1) : + ∃ y, ExecWithin w insert x y (11 * x.regs len + 4) ∧ + y.mem = memory x.mem keys (x.regs dst) (x.regs bins) (x.regs len) ∧ + y.regs total = pairs keys (x.regs len) ∧ + (∀ i, i.val < 4 → y.regs i = x.regs i) ∧ y.regs one = 1 := by + let I (y : Data 11) : Prop := + (∀ i, i.val < 4 → y.regs i = x.regs i) ∧ y.regs one = 1 ∧ + y.regs idx ≤ x.regs len ∧ y.regs total = pairs keys (y.regs idx) ∧ + y.mem = memory x.mem keys (x.regs dst) (x.regs bins) (y.regs idx) + have hbody (y : Data 11) (hy : I y) (hc : (Cond.lt idx len).Holds y) : + ∃ z, ExecWithin w body y z 9 ∧ I z ∧ + x.regs len - z.regs idx < x.regs len - y.regs idx := by + obtain ⟨hr, ho, hi, hp, hm⟩ := hy + have hl := hr len (by decide) + have hsrc := hr src (by decide) + have hdst := hr dst (by decide) + have hlt : y.regs idx < x.regs len := by simpa only [Cond.holds_lt, hl] using hc + have hkey := hkeys _ hlt + have hread : y.mem (y.regs src + y.regs idx) = keys (y.regs idx) := by + rw [hm, hsrc, memory_outside] + · exact hinput _ hlt + · omega + have hbucket : y.mem (y.regs dst + keys (y.regs idx)) = + count keys (y.regs idx) (keys (y.regs idx)) := by + rw [hm, hdst, memory_read _ _ _ _ _ _ hkey] + have ha : Fits w (y.regs src + y.regs idx) := hs.of_le (by rw [hsrc]; omega) + have hb : Fits w (y.regs dst + keys (y.regs idx)) := hd.of_le (by rw [hdst]; omega) + have hn : Fits w (y.regs idx + 1) := ht.of_le (by omega) + have hc := count_le keys (y.regs idx) (keys (y.regs idx)) + have hpnext := pairs_le_sq keys (y.regs idx + 1) + have hsquare : (y.regs idx + 1) * (y.regs idx + 1) ≤ x.regs len * x.regs len := + Nat.mul_self_le_mul_self (by omega) + have ht₂ : Fits w (pairs keys (y.regs idx) + count keys (y.regs idx) (keys (y.regs idx)) + + count keys (y.regs idx) (keys (y.regs idx))) := ht.of_le (by + dsimp only [pairs] at hpnext + omega) + have ht₁ : Fits w (pairs keys (y.regs idx) + count keys (y.regs idx) (keys (y.regs idx))) := + ht₂.of_le (by omega) + have hcount : Fits w (count keys (y.regs idx) (keys (y.regs idx)) + 1) := + ht.of_le (by omega) + refine ⟨?z, ?run, ?inv, ?decrease⟩ + case run => ram_eval [body, ha, hb, hn, ho, hread, hbucket, hp, ht₁, ht₂, hcount] + case inv => + refine ⟨?_, ?_, ?_, ?_, ?_⟩ + · intro i hi + have hia : i ≠ addr := by intro h; subst i; norm_num [addr] at hi + have hik : i ≠ key := by intro h; subst i; norm_num [key] at hi + have hic : i ≠ bucketSize := by intro h; subst i; norm_num [bucketSize] at hi + have hit : i ≠ total := by intro h; subst i; norm_num [total] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + simpa only [Data.set_regs, Data.write_regs, Function.update_of_ne hia, + Function.update_of_ne hik, Function.update_of_ne hic, Function.update_of_ne hit, + Function.update_of_ne hii] using hr i hi + · simpa +decide [Function.update_apply] using ho + · simpa [ho] using Nat.succ_le_of_lt hlt + · change pairs keys (y.regs idx) + count keys (y.regs idx) (keys (y.regs idx)) + + count keys (y.regs idx) (keys (y.regs idx)) = pairs keys (y.regs idx + 1) + simp only [pairs, two_mul, Nat.add_assoc] + · change Function.update y.mem (y.regs dst + keys (y.regs idx)) + (count keys (y.regs idx) (keys (y.regs idx)) + 1) = + memory x.mem keys (x.regs dst) (x.regs bins) (y.regs idx + 1) + rw [hm, hdst, ← memory_succ _ _ _ _ _ hkey] + case decrease => simp only [Data.set_regs, Function.update_self]; omega + let x₀ := (x.set idx 0).set total 0 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_⟩ + · intro i hi + have hit : i ≠ total := by intro h; subst i; norm_num [total] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + simp only [x₀, Data.set_regs, Function.update_of_ne hit, Function.update_of_ne hii] + · simpa +decide [x₀, one, total, idx] using ho + · simp +decide [x₀, total, idx] + · simp +decide [x₀, total, idx, pairs] + · change x.mem = memory x.mem keys (x.regs dst) (x.regs bins) 0 + funext a + by_cases h : x.regs dst ≤ a ∧ a < x.regs dst + x.regs bins + · have ha : x.regs dst + (a - x.regs dst) = a := by omega + simpa [memory, h, ha] using hzero (a - x.regs dst) (by omega) + · simp only [memory, h, ↓reduceIte] + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => x.regs len - y.regs idx) hbody hx₀ + obtain ⟨hr, ho, hi, hp, hm⟩ := hy + have hl := hr len (by decide) + have hend : y.regs idx = x.regs len := by + simp only [Cond.holds_lt, hl] at hc + omega + refine ⟨y, ?_, by simpa only [hend] using hm, by simpa only [hend] using hp, hr, ho⟩ + have hinit : ExecWithin (B := .standard) w (ram[ idx ← #0; total ← #0 ]) x x₀ 2 := by + ram_eval [x₀] + convert hinit.seq he using 1 + · rfl + · simp +decide only [x₀, Data.set_regs, Function.update_apply, idx, total, ↓reduceIte, + Nat.sub_zero] + omega + +/-- Build a histogram and count ordered collisions within `5M + 11n + 9` native instructions. +The input keys lie in `[0,M)`. The input precedes the output block. Both one-past-end +addresses and the quadratic counter bound must fit. These are explicit local contracts; +constructing the keys and deriving a logarithmic width are separate client obligations. -/ +theorem exec {w : ℕ} (x : Data 11) (keys : ℕ → ℕ) (h1 : Fits w 1) + (hs : Fits w (x.regs src + x.regs len)) + (hd : Fits w (x.regs dst + x.regs bins)) + (ht : Fits w (x.regs len * x.regs len + x.regs len + 1)) + (hdisj : x.regs src + x.regs len ≤ x.regs dst) + (hkeys : ∀ i < x.regs len, keys i < x.regs bins) + (hinput : ∀ i < x.regs len, x.mem (x.regs src + i) = keys i) : + ∃ y, ExecWithin w program x y (5 * x.regs bins + 11 * x.regs len + 9) ∧ + y.mem = memory x.mem keys (x.regs dst) (x.regs bins) (x.regs len) ∧ + y.regs total = pairs keys (x.regs len) ∧ + (∀ i, i.val < 4 → y.regs i = x.regs i) := by + obtain ⟨y, he, hm, hr, ho⟩ := exec_clear x h1 hd + have hsrc := hr src (by decide) + have hdst := hr dst (by decide) + have hlen := hr len (by decide) + have hbins := hr bins (by decide) + obtain ⟨z, hz, hzm, hzt, hzr, _⟩ := exec_insert y keys + (by simpa only [hsrc, hlen] using hs) (by simpa only [hdst, hbins] using hd) + (by simpa only [hlen] using ht) (by simpa only [hsrc, hlen, hdst] using hdisj) + (by simpa only [hlen, hbins] using hkeys) (by + intro i hi + rw [hsrc, hm] + have hlt : i < x.regs len := by simpa only [hlen] using hi + simp only [Cslib.WordRAM.clear, show ¬(x.regs dst ≤ x.regs src + i ∧ + x.regs src + i < x.regs dst + x.regs bins) from by omega, ↓reduceIte] + exact hinput i hlt) (by + intro b hb + rw [hm, hdst] + have hb' : b < x.regs bins := by simpa only [hbins] using hb + simp [Cslib.WordRAM.clear, hb']) ho + refine ⟨z, ?_, ?_, by simpa only [hlen] using hzt, fun i hi => (hzr i hi).trans (hr i hi)⟩ + · exact (he.seq hz).mono (by rw [hlen]; omega) + · rw [hzm, hm, hdst, hbins, hlen] + funext a + simp only [memory, Cslib.WordRAM.clear] + split_ifs <;> rfl + +end Cslib.Algorithms.WordRAM.Histogram diff --git a/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Collisions.lean b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Collisions.lean new file mode 100644 index 0000000000..6b18f5850d --- /dev/null +++ b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Collisions.lean @@ -0,0 +1,77 @@ +/- +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.Histogram +public import Complexity.Algorithms.ThreeSum.ChanHe.Collisions + +/-! +# Histogram collision counts for Chan–He + +The native histogram counts ordered pairs of different input positions. For an injective +array of integers and its Euclidean residue keys, this is exactly the existing Chan–He +collision count. Injectivity is required because the mathematical reduction operates on sets. +Equal input values must first be handled by the distinct-position frontend. + +## Main results + +* `coll_eq_histogramPairs`: the existing `coll` is the counter computed by the histogram. + +## References + +* [T. M. Chan and Q. He, *Reducing 3SUM to Convolution-3SUM*][ChanHe2020] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.ThreeSum.ChanHe + +open Finset +open Cslib.Algorithms.ThreeSum.ChanHe + +/-- For distinct array values, histogramming Euclidean residues computes exactly the +ordered distinct-element collision count of the existing Chan–He reduction. -/ +theorem coll_eq_histogramPairs (values : ℕ → ℤ) (n M : ℕ) (hM : 0 < M) + (hinj : Set.InjOn values (↑(range n) : Set ℕ)) : + coll ((range n).image values) M = + Histogram.pairs (fun i => (values i % (M : ℤ)).toNat) n := by + rw [Histogram.pairs_eq_card] + symm + refine card_bij (fun p _ => (values p.1, values p.2)) ?_ ?_ ?_ + · intro p hp + obtain ⟨hp, heq⟩ := mem_filter.mp hp + obtain ⟨hp₁, hp₂, hne⟩ := mem_offDiag.mp hp + refine mem_filter.mpr ⟨mem_offDiag.mpr ⟨mem_image.mpr ⟨p.1, hp₁, rfl⟩, + mem_image.mpr ⟨p.2, hp₂, rfl⟩, fun h => hne (hinj hp₁ hp₂ h)⟩, ?_⟩ + apply (natCast_dvd_sub_iff_emod_eq M _ _).mpr + change values p.1 % (M : ℤ) = values p.2 % (M : ℤ) + have hpos : (0 : ℤ) < M := by exact_mod_cast hM + have h₁ := Int.emod_nonneg (values p.1) (ne_of_gt hpos) + have h₂ := Int.emod_nonneg (values p.2) (ne_of_gt hpos) + omega + · intro p hp q hq heq + have hp' := mem_offDiag.mp (mem_filter.mp hp).1 + have hq' := mem_offDiag.mp (mem_filter.mp hq).1 + exact Prod.ext (hinj hp'.1 hq'.1 (congrArg Prod.fst heq)) + (hinj hp'.2.1 hq'.2.1 (congrArg Prod.snd heq)) + · intro p hp + obtain ⟨hp, hdvd⟩ := mem_filter.mp hp + obtain ⟨hp₁, hp₂, hne⟩ := mem_offDiag.mp hp + obtain ⟨i, hi, hival⟩ := mem_image.mp hp₁ + obtain ⟨j, hj, hjval⟩ := mem_image.mp hp₂ + refine ⟨(i, j), mem_filter.mpr ⟨mem_offDiag.mpr ⟨hi, hj, ?_⟩, ?_⟩, ?_⟩ + · intro hij + change i = j at hij + apply hne + rw [← hival, ← hjval, hij] + · have heq := (natCast_dvd_sub_iff_emod_eq M _ _).mp hdvd + change (values i % (M : ℤ)).toNat = (values j % (M : ℤ)).toNat + rw [hival, hjval, heq] + · exact Prod.ext hival hjval + +end Cslib.Algorithms.WordRAM.ThreeSum.ChanHe diff --git a/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Count.lean b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Count.lean new file mode 100644 index 0000000000..5a8a697325 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Count.lean @@ -0,0 +1,220 @@ +/- +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.ThreeSum.ChanHe.Collisions +public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Residues +import Mathlib.Tactic.FinCases +import Mathlib.Tactic.NormNum + +/-! +# Native Chan–He collision counting + +This fixed native program first decodes signed Euclidean residue keys, then builds a +histogram and counts ordered collisions. Register setup between the two loops is charged. +For injective input arrays its counter is exactly the existing Chan–He `coll`. + +## Main definitions + +* `Count.program`: the composed native residue/histogram program. +* `Count.memory`: the complete output memory, including the residue scratch array. + +## Main results + +* `Count.exec`: execution in at most `34n + 5M + 18` native instructions. +* `Count.exec_coll`: the counter equals the existing mathematical collision count. + +## References + +* [T. M. Chan and Q. He, *Reducing 3SUM to Convolution-3SUM*][ChanHe2020] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.ThreeSum.ChanHe.Count + +open Cslib.WordRAM Cslib.WordRAM.Signed + +/-- Residue routine registers occupy the first sixteen registers. -/ +def residuesRegisters (i : Fin 16) : Fin 28 := ⟨i.val, by omega⟩ + +/-- Histogram registers occupy registers sixteen through twenty-six. -/ +def histogramRegisters (i : Fin 11) : Fin 28 := ⟨16 + i.val, by omega⟩ + +private theorem residues_injective : Function.Injective residuesRegisters := by + intro i j h + apply Fin.ext + exact congrArg (fun i : Fin 28 => i.val) h + +private theorem histogram_injective : Function.Injective histogramRegisters := by + intro i j h + apply Fin.ext + have := congrArg Fin.val h + dsimp [histogramRegisters] at this + omega + +/-- Encoded integer input base. -/ +abbrev src : Fin 28 := 9 +/-- Residue scratch base. -/ +abbrev scratch : Fin 28 := 10 +/-- Number of input values. -/ +abbrev len : Fin 28 := 11 +/-- Encoded positive modulus. -/ +abbrev modulus : Fin 28 := 12 +/-- Histogram destination base. -/ +abbrev dst : Fin 28 := 27 +/-- Output collision count, an ordinary unsigned word. -/ +abbrev total : Fin 28 := 24 + +/-- Histogram input base register. -/ +abbrev hsrc : Fin 28 := histogramRegisters Histogram.src +/-- Histogram output base register. -/ +abbrev hdst : Fin 28 := histogramRegisters Histogram.dst +/-- Histogram length register. -/ +abbrev hlen : Fin 28 := histogramRegisters Histogram.len +/-- Histogram modulus register. -/ +abbrev hbins : Fin 28 := histogramRegisters Histogram.bins +/-- Histogram constant register, also used to decode the encoded modulus. -/ +abbrev hone : Fin 28 := histogramRegisters Histogram.one + +/-- Set up the histogram parameters, explicitly decoding the positive modulus. -/ +def prepare : Code .standard 28 := ram[ + hsrc ← scratch; hdst ← dst; hlen ← len; hone ← #2; hbins ← modulus / hone ] + +/-- A single finite native program, independent of input length, modulus and word width. -/ +def program : Code .standard 28 := + .seq (Residues.program.rename residuesRegisters) + (.seq prepare (Histogram.program.rename histogramRegisters)) + +/-- The complete counting program preserves its input base, scratch base, length, +encoded modulus and histogram base registers. -/ +theorem preserves_parameters (i : Fin 28) + (hi : (9 ≤ i.val ∧ i.val ≤ 12) ∨ i = dst) : program.Preserves i := by + fin_cases i <;> norm_num [dst] at hi <;> decide + +/-- The residue scratch array and the complete histogram, framed by the original memory. -/ +def memory (mem : ℕ → ℕ) (values : ℕ → ℤ) (n M scratch dst : ℕ) : ℕ → ℕ := + Histogram.memory (Residues.memory mem values M scratch n) + (fun i => (values i % (M : ℤ)).toNat) dst M n + +/-- Complete native collision counting, including residue computation, clearing the histogram +and register setup. Input, scratch and histogram blocks occur in that order and are disjoint. +The final halt is excluded. Fit hypotheses are local; a caller must derive them from its +input and width contract. -/ +theorem exec {w M : ℕ} (x : Data 28) (values : ℕ → ℤ) (hM : 0 < M) + (hmod : GenericRAM.Signed.Fits w (M : ℤ)) + (hmodreg : x.regs modulus = encode (M : ℤ)) + (hfit : ∀ i < x.regs len, GenericRAM.Signed.Fits w (values i)) + (hs : Fits (w + 2) (x.regs src + x.regs len)) + (hk : Fits (w + 2) (x.regs scratch + x.regs len)) + (hd : Fits (w + 2) (x.regs dst + M)) + (ht : Fits (w + 2) (x.regs len * x.regs len + x.regs len + 1)) + (hsk : x.regs src + x.regs len ≤ x.regs scratch) + (hkd : x.regs scratch + x.regs len ≤ x.regs dst) + (hinput : ∀ i < x.regs len, x.mem (x.regs src + i) = encode (values i)) : + ∃ y, ExecWithin (w + 2) program x y (34 * x.regs len + 5 * M + 18) ∧ + y.mem = memory x.mem values (x.regs len) M (x.regs scratch) (x.regs dst) ∧ + y.regs total = Histogram.pairs (fun i => (values i % (M : ℤ)).toNat) (x.regs len) ∧ + ∀ i, ((9 ≤ i.val ∧ i.val ≤ 12) ∨ i = dst) → y.regs i = x.regs i := by + obtain ⟨a, ha, ham, har⟩ := Residues.exec (x.comap residuesRegisters) values hM hmod + hmodreg hfit hs hk hsk hinput + let y := a.extend residuesRegisters x + have hy : ExecWithin (w + 2) (Residues.program.rename residuesRegisters) x y + (23 * x.regs len + 4) := ha.rename residues_injective + have hyr (i : Fin 16) : y.regs (residuesRegisters i) = a.regs i := + residues_injective.extend_apply _ _ i + have hys : y.regs scratch = x.regs scratch := + (hyr Residues.dst).trans (har Residues.dst (by decide)) + have hyl : y.regs len = x.regs len := (hyr Residues.len).trans (har Residues.len (by decide)) + have hym : y.regs modulus = encode (M : ℤ) := + (hyr Residues.modulus).trans ((har Residues.modulus (by decide)).trans hmodreg) + have hyd : y.regs dst = x.regs dst := by + change Function.extend residuesRegisters a.regs x.regs dst = x.regs dst + apply Function.extend_apply' + rintro ⟨i, hi⟩ + have := congrArg Fin.val hi + dsimp [residuesRegisters, dst] at this + omega + let q := ((((y.set (histogramRegisters Histogram.src) (x.regs scratch)).set + (histogramRegisters Histogram.dst) (x.regs dst)).set + (histogramRegisters Histogram.len) (x.regs len)).set + (histogramRegisters Histogram.one) 2).set (histogramRegisters Histogram.bins) M + have htwo : Fits (w + 2) 2 := (show Fits 2 2 from by decide).mono (by omega) + have hraw : Fits (w + 2) M := + ((encode_fits hmod).mono (by omega)).of_le (by + simp only [encode, Int.natAbs_natCast] + omega) + have hprep : ExecWithin (w + 2) prepare y q 5 := by + ram_eval [prepare, q, hys, hyl, hym, hyd, htwo, hraw, histogramRegisters, + Histogram.src, Histogram.dst, Histogram.len, Histogram.one, Histogram.bins, + hsrc, hdst, hlen, hbins, hone, scratch, len, modulus, dst, Fin.mk.injEq, + encode_div_two, Int.natAbs_natCast] + have hqs : (q.comap histogramRegisters).regs Histogram.src = x.regs scratch := by + simp +decide [q, histogramRegisters, Histogram.src, Histogram.dst, Histogram.len, + Histogram.one, Histogram.bins] + have hqd : (q.comap histogramRegisters).regs Histogram.dst = x.regs dst := by + simp +decide [q, histogramRegisters, Histogram.src, Histogram.dst, Histogram.len, + Histogram.one, Histogram.bins] + have hql : (q.comap histogramRegisters).regs Histogram.len = x.regs len := by + simp +decide [q, histogramRegisters, Histogram.src, Histogram.dst, Histogram.len, + Histogram.one, Histogram.bins] + have hqb : (q.comap histogramRegisters).regs Histogram.bins = M := by + simp [q] + have hqm : q.mem = Residues.memory x.mem values M (x.regs scratch) (x.regs len) := ham + have h1 : Fits (w + 2) 1 := (show Fits 2 1 from by decide).mono (by omega) + obtain ⟨b, hb, hbm, hbt, _⟩ := Histogram.exec (q.comap histogramRegisters) + (fun i => (values i % (M : ℤ)).toNat) h1 + (by simpa only [hqs, hql] using hk) (by simpa only [hqd, hqb] using hd) + (by simpa only [hql] using ht) (by simpa only [hqs, hql, hqd] using hkd) + (by + intro i hi + rw [hqb] + have hnonneg := Int.emod_nonneg (values i) (show (M : ℤ) ≠ 0 by omega) + have hlt := Int.emod_lt_of_pos (values i) (show (0 : ℤ) < M by omega) + omega) + (by + intro i hi + have hi' : i < x.regs len := by simpa only [hql] using hi + rw [Data.comap_mem, hqm, hqs] + simp [Residues.memory, hi']) + have he := (hy.seq (hprep.seq (hb.rename histogram_injective))).mono + (show (23 * x.regs len + 4) + (5 + (5 * (q.comap histogramRegisters).regs Histogram.bins + + 11 * (q.comap histogramRegisters).regs Histogram.len + 9)) ≤ + 34 * x.regs len + 5 * M + 18 by rw [hql, hqb]; omega) + refine ⟨b.extend histogramRegisters q, he, ?_, ?_, fun i hi => + ExecWithin.regs_eq he (preserves_parameters i hi)⟩ + · change b.mem = _ + rw [hbm, Data.comap_mem, hqm, hqd, hqb, hql] + rfl + · change Function.extend histogramRegisters b.regs q.regs + (histogramRegisters Histogram.total) = _ + rw [histogram_injective.extend_apply, hbt, hql] + +/-- For an injective integer array, the composed program returns exactly the collision +count used by the existing Chan–He prime searches. -/ +theorem exec_coll {w M : ℕ} (x : Data 28) (values : ℕ → ℤ) (hM : 0 < M) + (hmod : GenericRAM.Signed.Fits w (M : ℤ)) + (hmodreg : x.regs modulus = encode (M : ℤ)) + (hfit : ∀ i < x.regs len, GenericRAM.Signed.Fits w (values i)) + (hs : Fits (w + 2) (x.regs src + x.regs len)) + (hk : Fits (w + 2) (x.regs scratch + x.regs len)) + (hd : Fits (w + 2) (x.regs dst + M)) + (ht : Fits (w + 2) (x.regs len * x.regs len + x.regs len + 1)) + (hsk : x.regs src + x.regs len ≤ x.regs scratch) + (hkd : x.regs scratch + x.regs len ≤ x.regs dst) + (hinput : ∀ i < x.regs len, x.mem (x.regs src + i) = encode (values i)) + (hinj : Set.InjOn values (↑(Finset.range (x.regs len)) : Set ℕ)) : + ∃ y, ExecWithin (w + 2) program x y (34 * x.regs len + 5 * M + 18) ∧ + y.mem = memory x.mem values (x.regs len) M (x.regs scratch) (x.regs dst) ∧ + y.regs total = Cslib.Algorithms.ThreeSum.ChanHe.coll + ((Finset.range (x.regs len)).image values) M ∧ + ∀ i, ((9 ≤ i.val ∧ i.val ≤ 12) ∨ i = dst) → y.regs i = x.regs i := by + obtain ⟨y, hy, hm, ht, hr⟩ := exec x values hM hmod hmodreg hfit hs hk hd ht hsk hkd hinput + exact ⟨y, hy, hm, ht.trans (coll_eq_histogramPairs values _ M hM hinj).symm, hr⟩ + +end Cslib.Algorithms.WordRAM.ThreeSum.ChanHe.Count diff --git a/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Residues.lean b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Residues.lean new file mode 100644 index 0000000000..4cdbd3eef3 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/ThreeSum/ChanHe/Residues.lean @@ -0,0 +1,211 @@ +/- +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.Signed.Arithmetic +public import Complexity.Computability.WordRAM.Copy +import Mathlib.Tactic.NormNum + +/-! +# Native Euclidean residue arrays + +The input holds signed encodings. The output holds ordinary unsigned residue keys, ready +for histogram construction. In particular a nonnegative encoded integer is explicitly decoded +by division by two before it is used as a histogram index. + +## Main definitions + +* `Residues.program`: a fixed loop producing Euclidean residue keys. + +## Main results + +* `Residues.exec`: charged execution preserving the input and all cells outside the output. + +## References + +* [T. M. Chan and Q. He, *Reducing 3SUM to Convolution-3SUM*][ChanHe2020] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.ThreeSum.ChanHe.Residues + +open Cslib.WordRAM Cslib.WordRAM.Signed + +/-- Encoded integer input base. -/ +abbrev src : Fin 16 := 9 +/-- Unsigned residue output base. -/ +abbrev dst : Fin 16 := 10 +/-- Input length. -/ +abbrev len : Fin 16 := 11 +/-- Encoded modulus. -/ +abbrev modulus : Fin 16 := 12 +/-- Loop counter. -/ +abbrev idx : Fin 16 := 13 +/-- Address scratch. -/ +abbrev addr : Fin 16 := 14 +/-- Constant one. -/ +abbrev one : Fin 16 := 15 + +/-- The signed remainder block followed by explicit unsigned decoding and a store. -/ +def body : Code .standard 16 := + .seq (ram[ addr ← src + idx; left ← *addr; right ← modulus ]) + (.seq Signed.modulo (ram[ two ← #2; result ← result / two; + addr ← dst + idx; *addr ← result; idx ← idx + one ])) + +/-- All parameters are supplied at runtime; registers and code are independent +of length and width. -/ +def program : Code .standard 16 := + .seq (ram[ idx ← #0; one ← #1 ]) (.loop (.lt idx len) body) + +/-- Memory after writing a prefix of unsigned residue keys. -/ +def memory (mem : ℕ → ℕ) (values : ℕ → ℤ) (M dst k : ℕ) : ℕ → ℕ := + fun a => if dst ≤ a ∧ a < dst + k then (values (a - dst) % (M : ℤ)).toNat else mem a + +@[simp] theorem memory_zero (mem : ℕ → ℕ) (values : ℕ → ℤ) (M dst : ℕ) : + memory mem values M dst 0 = mem := by + funext a + simp [memory] + +private theorem memory_succ (mem : ℕ → ℕ) (values : ℕ → ℤ) (M dst k : ℕ) : + memory mem values M dst (k + 1) = Function.update (memory mem values M dst k) + (dst + k) (values k % (M : ℤ)).toNat := by + funext a + by_cases ha : a = dst + k + · subst a + simp [memory] + · rw [Function.update_of_ne ha] + simp only [memory] + have h : (dst ≤ a ∧ a < dst + (k + 1)) ↔ (dst ≤ a ∧ a < dst + k) := by omega + simp only [h] + +/-- Convert `n` encoded signed integers to unsigned Euclidean residue keys in at most +`23n + 4` native instructions. The modulus is positive and encoded in its parameter register; +addresses are ordinary unsigned words. The final halt is excluded. -/ +theorem exec {w M : ℕ} (x : Data 16) (values : ℕ → ℤ) (hM : 0 < M) + (hmod : GenericRAM.Signed.Fits w (M : ℤ)) + (hmodreg : x.regs modulus = encode (M : ℤ)) + (hfit : ∀ i < x.regs len, GenericRAM.Signed.Fits w (values i)) + (hs : Fits (w + 2) (x.regs src + x.regs len)) + (hd : Fits (w + 2) (x.regs dst + x.regs len)) + (hdisj : x.regs src + x.regs len ≤ x.regs dst) + (hinput : ∀ i < x.regs len, x.mem (x.regs src + i) = encode (values i)) : + ∃ y, ExecWithin (w + 2) program x y (23 * x.regs len + 4) ∧ + y.mem = memory x.mem values M (x.regs dst) (x.regs len) ∧ + (∀ i, 9 ≤ i.val ∧ i.val ≤ 12 → y.regs i = x.regs i) := by + let I (y : Data 16) : Prop := + (∀ i, 9 ≤ i.val ∧ i.val ≤ 12 → y.regs i = x.regs i) ∧ y.regs one = 1 ∧ + y.regs idx ≤ x.regs len ∧ + y.mem = memory x.mem values M (x.regs dst) (y.regs idx) + have hbody (y : Data 16) (hy : I y) (hc : (Cond.lt idx len).Holds y) : + ∃ z, ExecWithin (w + 2) body y z 21 ∧ I z ∧ + x.regs len - z.regs idx < x.regs len - y.regs idx := by + obtain ⟨hr, ho, hi, hm⟩ := hy + have hl := hr len (by decide) + have hsrc := hr src (by decide) + have hdst := hr dst (by decide) + have hmod' := hr modulus (by decide) + have hlt : y.regs idx < x.regs len := by simpa only [Cond.holds_lt, hl] using hc + have hread : y.mem (y.regs src + y.regs idx) = encode (values (y.regs idx)) := by + rw [hm, hsrc] + simp only [memory, show ¬(x.regs dst ≤ x.regs src + y.regs idx ∧ + x.regs src + y.regs idx < x.regs dst + y.regs idx) from by omega, ↓reduceIte] + exact hinput _ hlt + have ha : Fits (w + 2) (y.regs src + y.regs idx) := hs.of_le (by rw [hsrc]; omega) + let p := ((y.set addr (y.regs src + y.regs idx)).set left + (encode (values (y.regs idx)))).set right (encode (M : ℤ)) + have hp : ExecWithin (w + 2) + (ram[ addr ← src + idx; left ← *addr; right ← modulus ]) y p 3 := by + ram_eval [p, ha, hread, hmod', hmodreg, addr, src, idx, left, right, modulus, Fin.mk.injEq] + obtain ⟨z, hz, hzval, hzmem, hzreg⟩ := Signed.exec_modulo (hfit _ hlt) hmod p + (by simp +decide [p, left, right]) (by simp [p]) + have hzr (i : Fin 16) (hi : 9 ≤ i.val) : z.regs i = p.regs i := hzreg i (Or.inr hi) + have hzi : z.regs idx = y.regs idx := by + simpa +decide [p, idx, addr, left, right] using hzr idx (by decide) + have hzo : z.regs one = 1 := by + simpa +decide [p, one, addr, left, right] using (hzr one (by decide)).trans + (show p.regs one = 1 by simpa +decide [p, one, addr, left, right] using ho) + have hzd : z.regs dst = x.regs dst := by + simpa +decide [p, dst, addr, left, right] using (hzr dst (by decide)).trans + (show p.regs dst = x.regs dst by simpa +decide [p, dst, addr, left, right] using hdst) + have hnonneg := Int.emod_nonneg (values (y.regs idx)) (show (M : ℤ) ≠ 0 by omega) + have hdecode : z.regs result / 2 = (values (y.regs idx) % (M : ℤ)).toNat := by + rw [hzval, GenericRAM.Signed.Arithmetic.eval, encode_div_two] + omega + have hb : Fits (w + 2) (z.regs dst + z.regs idx) := hd.of_le (by rw [hzd, hzi]; omega) + have hn : Fits (w + 2) (z.regs idx + 1) := hs.of_le (by rw [hzi]; omega) + have htwo : Fits (w + 2) 2 := (show Fits 2 2 from by decide).mono (by omega) + have hrem : Fits (w + 2) (z.regs result / 2) := by + rw [hdecode] + have hlt' := Int.emod_lt_of_pos (values (y.regs idx)) (show (0 : ℤ) < M by omega) + exact ((encode_fits hmod).mono (by omega)).of_le (by + simp only [encode, Int.natAbs_natCast, show ¬(M : ℤ) < 0 from by omega, ↓reduceIte] + omega) + refine ⟨?q, ?run, ?inv, ?decrease⟩ + case run => + have htail : ExecWithin (w + 2) + (ram[ two ← #2; result ← result / two; addr ← dst + idx; + *addr ← result; idx ← idx + one ]) z ?q 5 := by + ram_eval [htwo, hrem, hb, hn, hzo, two, result, addr, dst, idx, one, Fin.mk.injEq] + exact (hp.seq (hz.seq htail)) + case inv => + refine ⟨?_, ?_, ?_, ?_⟩ + · intro i hi + have hip : p.regs i = y.regs i := by + have hia : i ≠ addr := by intro h; subst i; norm_num [addr] at hi + have hil : i ≠ left := by intro h; subst i; norm_num [left] at hi + have hir : i ≠ right := by intro h; subst i; norm_num [right] at hi + simp only [p, Data.set_regs, Function.update_of_ne hia, Function.update_of_ne hil, + Function.update_of_ne hir] + have hia : i ≠ addr := by intro h; subst i; norm_num [addr] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + have hit : i ≠ two := by intro h; subst i; norm_num [two] at hi + have hir : i ≠ result := by intro h; subst i; norm_num [result] at hi + simpa only [Data.set_regs, Data.write_regs, Function.update_of_ne hia, + Function.update_of_ne hii, Function.update_of_ne hit, Function.update_of_ne hir] + using (hzr i hi.1).trans (hip.trans (hr i hi)) + · simpa +decide [Function.update_apply] using hzo + · simpa +decide [Function.update_apply, hzi, hzo] using Nat.succ_le_of_lt hlt + · simp only [Data.set_mem, Data.write_mem, Data.set_regs, Function.update_apply] + rw [hdecode, hzmem] + change Function.update y.mem (z.regs dst + z.regs idx) + (values (y.regs idx) % (M : ℤ)).toNat = _ + rw [hm, hzd, hzi, ← memory_succ] + rfl + case decrease => + change x.regs len - (z.regs idx + 1) < x.regs len - y.regs idx + rw [hzi] + omega + let x₀ := (x.set idx 0).set one 1 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_⟩ + · intro i hi + have hio : i ≠ one := by intro h; subst i; norm_num [one] at hi + have hii : i ≠ idx := by intro h; subst i; norm_num [idx] at hi + simp only [x₀, Data.set_regs, Function.update_of_ne hio, Function.update_of_ne hii] + · simp [x₀] + · simp +decide [x₀, idx, one] + · simp +decide [x₀, idx, one] + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => x.regs len - y.regs idx) hbody hx₀ + obtain ⟨hr, ho, hi, hm⟩ := hy + have hl := hr len (by decide) + have hend : y.regs idx = x.regs len := by + simp only [Cond.holds_lt, hl] at hc + omega + refine ⟨y, ?_, by simpa only [hend] using hm, hr⟩ + have h1 : Fits (w + 2) 1 := (show Fits 2 1 from by decide).mono (by omega) + have hinit : ExecWithin (B := .standard) (w + 2) (ram[ idx ← #0; one ← #1 ]) x x₀ 2 := by + ram_eval [x₀, h1] + convert hinit.seq he using 1 + · rfl + · simp +decide only [x₀, Data.set_regs, Function.update_apply, idx, one, ↓reduceIte, + Nat.sub_zero] + omega + +end Cslib.Algorithms.WordRAM.ThreeSum.ChanHe.Residues diff --git a/Complexity/Foundations/Data/Finset/Quotient.lean b/Complexity/Foundations/Data/Finset/Quotient.lean new file mode 100644 index 0000000000..875f7e52d7 --- /dev/null +++ b/Complexity/Foundations/Data/Finset/Quotient.lean @@ -0,0 +1,44 @@ +/- +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.Init +public import Mathlib.Order.Interval.Finset.Nat + +/-! +# Cardinality of a quotient fiber + +## Main results + +* `Finset.card_filter_div_eq_le`: an injectively indexed quotient fiber has at most the divisor + many members. + +## References + +Adapted from [Counting.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Util/Counting.lean) at +`anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +namespace Finset + +/-- At most `c` elements `x` of `s` have `⌊f x / c⌋ = k`, if `f` is injective on `s`: their values +lie in `[kc, kc + c)`. -/ +theorem card_filter_div_eq_le {α : Type*} (s : Finset α) {f : α → ℕ} (hf : Set.InjOn f s) + {c : ℕ} (hc : 0 < c) (k : ℕ) : (s.filter fun x => f x / c = k).card ≤ c := by + have hmaps : ∀ x ∈ s.filter fun x => f x / c = k, f x ∈ Finset.Ico (k * c) (k * c + c) := by + intro x hx + obtain ⟨-, rfl⟩ := Finset.mem_filter.mp hx + rw [Finset.mem_Ico, ← Nat.succ_mul] + exact ⟨Nat.div_mul_le_self _ _, (Nat.div_lt_iff_lt_mul hc).mp (Nat.lt_succ_self _)⟩ + calc (s.filter fun x => f x / c = k).card ≤ (Finset.Ico (k * c) (k * c + c)).card := + Finset.card_le_card_of_injOn f hmaps (hf.mono fun x hx => (Finset.mem_filter.mp hx).1) + _ = c := by rw [Nat.card_Ico, Nat.add_sub_cancel_left] + +end Finset diff --git a/Complexity/Foundations/Data/Int/Residue.lean b/Complexity/Foundations/Data/Int/Residue.lean new file mode 100644 index 0000000000..eee6215f92 --- /dev/null +++ b/Complexity/Foundations/Data/Int/Residue.lean @@ -0,0 +1,39 @@ +/- +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.Init +public import Mathlib.Algebra.Order.Ring.Int + +/-! +# Integer residues as natural numbers + +## Main results + +* `Int.natCast_toNat_emod`: the nonnegative integer residue survives a round trip through `Nat`. + +Bounds on natural residues use Mathlib's existing `Int.natMod_lt`. + +## References + +Adapted from [Index.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Util/Index.lean) at +`anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +The existing Exact Triangle hashing client now uses these shared facts. +-/ + +@[expose] public section + +namespace Int + +/-- Casting the natural Euclidean residue back to `Int` preserves its value for a positive +modulus. -/ +theorem natCast_toNat_emod {p : ℕ} (hp : 0 < p) (u : ℤ) : + ((u % (p : ℤ)).toNat : ℤ) = u % (p : ℤ) := + Int.toNat_of_nonneg (Int.emod_nonneg _ (by omega)) + +end Int diff --git a/Complexity/Foundations/Data/Matrix/PairChunks.lean b/Complexity/Foundations/Data/Matrix/PairChunks.lean new file mode 100644 index 0000000000..3f3c74fee3 --- /dev/null +++ b/Complexity/Foundations/Data/Matrix/PairChunks.lean @@ -0,0 +1,114 @@ +/- +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.Foundations.Data.Nat.CeilDiv +public import Complexity.Foundations.Data.Finset.Quotient +public import Mathlib.Logic.Equiv.Fin.Basic +public import Mathlib.Data.Fintype.Prod +import Mathlib.Order.Interval.Finset.Fin + +/-! +# Chunks of matrix index pairs + +## Main definitions + +* `pairIndex`, `rankIn`: row-major index and rank in a finite set. +* `chunk`, `numChunks`: consecutive chunks and their number. + +## Main results + +* `card_chunk_le`: each chunk has at most the capacity many pairs. +* `biUnion_chunk`: all chunks cover the input finite set. + +## References + +Adapted from [Chunks.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Chunks.lean) and the minimal definitions in +`PaperStatements.lean` at `anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +There are no graph or weight assumptions. +-/ + +@[expose] public section + +namespace Cslib.Matrix + +/-- The position of the pair `(a, b)` in the row-major order of `A × B`. -/ +def pairIndex {n : ℕ} (q : Fin n × Fin n) : ℕ := q.1.val * n + q.2.val + +/-- The number of pairs of `S` that come before `q` in row-major order. -/ +def rankIn {n : ℕ} (S : Finset (Fin n × Fin n)) (q : Fin n × Fin n) : ℕ := + (S.filter fun r => pairIndex r < pairIndex q).card + +/-- The `j`-th chunk (`j = 0, 1, …`) of `S` when `S` is cut, in row-major order, into chunks of +`cap` pairs each (the last one may be smaller). -/ +def chunk {n : ℕ} (S : Finset (Fin n × Fin n)) (cap j : ℕ) : Finset (Fin n × Fin n) := + S.filter fun q => rankIn S q / cap = j + +/-- The number `⌈|S| / cap⌉` of nonempty chunks of `S` (for `cap ≥ 1`). -/ +def numChunks {n : ℕ} (S : Finset (Fin n × Fin n)) (cap : ℕ) : ℕ := + S.card ⌈/⌉ cap + +variable {n : ℕ} {S : Finset (Fin n × Fin n)} {cap : ℕ} + +/-! ### Ranks -/ + +/-- The position in row-major order determines the pair. -/ +theorem pairIndex_injective : Function.Injective (pairIndex (n := n)) := fun q r h => + finProdFinEquiv.injective (Fin.ext (by simpa [pairIndex, add_comm, mul_comm] using h)) + +/-- Among the pairs of `S`, a pair that comes earlier in row-major order has a smaller rank. -/ +theorem rankIn_lt_rankIn {q r : Fin n × Fin n} (hq : q ∈ S) + (h : pairIndex q < pairIndex r) : rankIn S q < rankIn S r := by + refine Finset.card_lt_card ⟨fun x hx => ?_, fun hsub => ?_⟩ + · rw [Finset.mem_filter] at hx ⊢ + exact ⟨hx.1, hx.2.trans h⟩ + · exact lt_irrefl _ (Finset.mem_filter.mp (hsub (Finset.mem_filter.mpr ⟨hq, h⟩))).2 + +/-- Different pairs of `S` have different ranks. -/ +theorem rankIn_injOn (S : Finset (Fin n × Fin n)) : Set.InjOn (rankIn S) S := by + intro q hq r hr h + rcases lt_trichotomy (pairIndex q) (pairIndex r) with hlt | heq | hgt + · exact absurd h (rankIn_lt_rankIn hq hlt).ne + · exact pairIndex_injective heq + · exact absurd h (rankIn_lt_rankIn hr hgt).ne' + +/-- The rank of a pair of `S` is less than the number of pairs of `S`. -/ +theorem rankIn_lt_card {q : Fin n × Fin n} (hq : q ∈ S) : rankIn S q < S.card := + Finset.card_lt_card ⟨Finset.filter_subset _ _, + fun hsub => lt_irrefl _ (Finset.mem_filter.mp (hsub hq)).2⟩ + +/-- In the full matrix, row-major rank is the row-major index itself. This avoids enumerating +all preceding pairs when applying the chunk API to a full matrix. -/ +theorem rankIn_univ (q : Fin n × Fin n) : + rankIn (Finset.univ : Finset (Fin n × Fin n)) q = pairIndex q := by + calc rankIn (Finset.univ : Finset (Fin n × Fin n)) q + = (Finset.Iio (finProdFinEquiv q)).card := by + apply Finset.card_equiv finProdFinEquiv + intro r + simp [pairIndex, finProdFinEquiv, add_comm, mul_comm, Fin.lt_def] + _ = pairIndex q := by simp [pairIndex, finProdFinEquiv, add_comm, mul_comm] + +/-! ### Chunks -/ + +/-- A chunk has at most `cap` pairs: their ranks are different numbers in `[j cap, (j + 1) cap)`. -/ +theorem card_chunk_le (S : Finset (Fin n × Fin n)) (hcap : 1 ≤ cap) (j : ℕ) : + (chunk S cap j).card ≤ cap := + S.card_filter_div_eq_le (rankIn_injOn S) hcap j + +/-- The first `⌈|S|/cap⌉` chunks cover `S`, because `|S| ≤ ⌈|S|/cap⌉ cap`. -/ +theorem biUnion_chunk (S : Finset (Fin n × Fin n)) (hcap : 1 ≤ cap) : + (Finset.range (numChunks S cap)).biUnion (fun j => chunk S cap j) = S := by + have hcard : S.card ≤ numChunks S cap * cap := by + rw [numChunks] + exact Nat.le_ceilDiv_mul hcap + ext q + simp only [Finset.mem_biUnion, Finset.mem_range, chunk, Finset.mem_filter] + refine ⟨fun ⟨_, _, hq, _⟩ => hq, fun hq => ⟨rankIn S q / cap, ?_, hq, rfl⟩⟩ + exact (Nat.div_lt_iff_lt_mul hcap).mpr ((rankIn_lt_card hq).trans_le hcard) + +end Cslib.Matrix diff --git a/Complexity/Foundations/Data/Nat/CeilDiv.lean b/Complexity/Foundations/Data/Nat/CeilDiv.lean new file mode 100644 index 0000000000..c88317bc1b --- /dev/null +++ b/Complexity/Foundations/Data/Nat/CeilDiv.lean @@ -0,0 +1,55 @@ +/- +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.Init +public import Mathlib.Algebra.Order.Floor.Div +public import Mathlib.Algebra.Order.Archimedean.Real.Basic + +/-! +# Natural ceilings of quotients + +## Main results + +* `Nat.ceil_div_eq_ceilDiv`: real quotient ceilings equal natural ceiling division. +* `Nat.le_ceilDiv_mul`, `Nat.ceilDiv_mul_lt`: covering and overshoot bounds. + +## References + +Adapted from [Ceil.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Util/Ceil.lean) at +`anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +namespace Nat + +/-! ## The ceiling of a quotient of natural numbers -/ + +/-- `a ⌈/⌉ b ≤ k` says that `k` pieces of size `b` cover `a`. Mathlib's `ceilDiv_le_iff_le_mul` has +`b * k` on the right. -/ +theorem ceilDiv_le_iff {a b k : ℕ} (hb : 0 < b) : a ⌈/⌉ b ≤ k ↔ a ≤ k * b := by + rw [ceilDiv_le_iff_le_mul hb, Nat.mul_comm] + +/-- `a ⌈/⌉ b` pieces of size `b` cover `a`. -/ +theorem le_ceilDiv_mul {a b : ℕ} (hb : 0 < b) : a ≤ a ⌈/⌉ b * b := + (ceilDiv_le_iff hb).1 le_rfl + +/-- `a ⌈/⌉ b` pieces of size `b` overshoot `a` by less than one piece. -/ +theorem ceilDiv_mul_lt {a b : ℕ} (hb : 0 < b) : a ⌈/⌉ b * b < a + b := by + have := Nat.div_mul_le_self (a + b - 1) b + rw [Nat.ceilDiv_eq_add_pred_div] + omega + +/-- The ceiling of the real quotient of two natural numbers, in natural numbers. -/ +theorem ceil_div_eq_ceilDiv (a : ℕ) {b : ℕ} (hb : 0 < b) : ⌈(a : ℝ) / (b : ℝ)⌉₊ = a ⌈/⌉ b := by + refine eq_of_forall_ge_iff fun k => ?_ + rw [Nat.ceil_le, div_le_iff₀ (Nat.cast_pos.2 hb), ceilDiv_le_iff hb] + exact_mod_cast Iff.rfl + +end Nat diff --git a/Complexity/Foundations/Data/Nat/SqrtRatio.lean b/Complexity/Foundations/Data/Nat/SqrtRatio.lean new file mode 100644 index 0000000000..602c737373 --- /dev/null +++ b/Complexity/Foundations/Data/Nat/SqrtRatio.lean @@ -0,0 +1,47 @@ +/- +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.Analysis.Real.Sqrt +import Mathlib.Tactic.Positivity +import Mathlib.Tactic.NormNum + +/-! +# Computing a floor involving a square root + +## Main results + +* `Nat.sqrt_div_eq_floor_div_sqrt`: the floor of `a / sqrt b` is the natural square root + of `a² / b`, including the zero-denominator convention. + +This elementary arithmetic bridge is original adaptation code. It lets finite algorithms use +integer arithmetic while retaining their real-valued capacity specifications. +-/ + +public section + +namespace Nat + +/-- A computable integer formula for the floor of a quotient by a square root. -/ +theorem sqrt_div_eq_floor_div_sqrt (a b : ℕ) : + Nat.sqrt (a ^ 2 / b) = ⌊(a : ℝ) / Real.sqrt b⌋₊ := by + by_cases hb : b = 0 + · simp [hb] + have hb0 : 0 < b := Nat.pos_of_ne_zero hb + have hs : 0 < Real.sqrt (b : ℝ) := Real.sqrt_pos.mpr (by exact_mod_cast hb0) + apply eq_of_forall_le_iff + intro k + rw [Nat.le_sqrt', Nat.le_div_iff_mul_le hb0, Nat.le_floor_iff (by positivity), + le_div_iff₀ hs] + have hcast : k ^ 2 * b ≤ a ^ 2 ↔ (k : ℝ) ^ 2 * (b : ℝ) ≤ (a : ℝ) ^ 2 := by + norm_cast + rw [hcast, ← sq_le_sq₀ (show (0 : ℝ) ≤ k * Real.sqrt b by positivity) (Nat.cast_nonneg a)] + rw [mul_pow, Real.sq_sqrt (Nat.cast_nonneg b)] + +end Nat diff --git a/ComplexityTests/ExactTriangleSplitting.lean b/ComplexityTests/ExactTriangleSplitting.lean new file mode 100644 index 0000000000..1768d67f64 --- /dev/null +++ b/ComplexityTests/ExactTriangleSplitting.lean @@ -0,0 +1,158 @@ +/- +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.ExactTriangle.Witnesses +import Mathlib.Tactic.NormNum.NatSqrt + +open Cslib.Algorithms.ExactTriangle Cslib.Matrix + +-- Row-major chunks, including the short last chunk and the first empty chunk. +example : chunk (Finset.univ : Finset (Fin 2 × Fin 2)) 3 0 = + {(0, 0), (0, 1), (1, 0)} := by decide +example : chunk (Finset.univ : Finset (Fin 2 × Fin 2)) 3 1 = {(1, 1)} := by decide +example : chunk (Finset.univ : Finset (Fin 2 × Fin 2)) 3 2 = ∅ := by decide +example : chunk (∅ : Finset (Fin 2 × Fin 2)) 3 0 = ∅ := by decide +example : numChunks (Finset.univ : Finset (Fin 2 × Fin 2)) 3 = 2 := by + rw [numChunks, Nat.ceilDiv_eq_add_pred_div] + decide + +private theorem pieceSize16 : pieceSize 16 2 = 2 := by + rw [pieceSize, sOf] + norm_num [Nat.ceilDiv_eq_add_pred_div] +private theorem queryCap16 : queryCap 16 16 = 64 := by norm_num [queryCap] +example : queryCap 17 17 = 70 := by norm_num [queryCap] +example : queryCap 0 16 = 0 := by norm_num [queryCap] +example : queryCap 16 0 = 0 := by norm_num [queryCap] + +-- Boundary indices belong to adjacent pieces; out-of-range pieces are empty. +example : (1 : Fin 16) ∈ piece 16 16 2 0 ∧ (2 : Fin 16) ∈ piece 16 16 2 1 := by + simp [piece, pieceSize16] +example : piece 16 16 2 8 = ∅ := by + ext c + simp only [piece, pieceSize16, Finset.mem_filter, Finset.mem_univ, true_and, + Finset.notMem_empty, iff_false] + omega +example : piece 0 16 2 0 = ∅ := by simp [piece] +example : numPieces 16 16 2 = 8 := by + rw [numPieces, pieceSize16, Nat.ceilDiv_eq_add_pred_div] + +private def zeroWeight : Fin 16 → Fin 16 → ℤ := fun _ _ => 0 +private def falseWeight : Fin 16 → Fin 16 → ℤ := fun _ _ => 6 + +-- The prime window admits p = 3 at the smallest supported D. +private theorem prime3 : 3 ∈ primesInRange 16 := by + rw [mem_primesInRange_iff_sq] + decide + +-- All-zero matrices supply an admitted, nonempty splitting with the promised bounds. +example : (totalChunks zeroWeight 16 3 : ℝ) ≤ 3 + Real.sqrt 16 := + (totalChunks_le zeroWeight (by decide) (by decide) prime3).1 +example : ((instanceIndices zeroWeight 16 2 3).card : ℝ) ≤ 4 * 16 * 2 := + card_instanceIndices_le zeroWeight (by decide) (by decide) (by decide) prime3 +example : Fintype.card (lopInstance zeroWeight zeroWeight zeroWeight 16 2 3 (0, 0, 0)).M ≤ 16 := + card_middle_lopInstance zeroWeight zeroWeight zeroWeight (by decide) (by decide) prime3 _ + +-- Residues include negative weights, and are independent of exact zero sums. +example : (0, 0) ∈ residueClass (fun _ _ : Fin 1 => -1) 3 2 := by decide +example : IsFalsePositive falseWeight zeroWeight zeroWeight 3 (0, 0, 0) := by decide +example : ¬ HasZeroTriangle falseWeight zeroWeight zeroWeight := by + simp [HasZeroTriangle, falseWeight, zeroWeight] +example : HasZeroTriangle zeroWeight zeroWeight zeroWeight := ⟨0, 0, 0, rfl⟩ + +-- Recovery rejects a modular false positive; a true zero returns a witness on its first scan. +example : scanPiece falseWeight zeroWeight zeroWeight 16 2 0 0 0 = none := by + rw [scanPiece_eq_none_iff] + intro c _ + norm_num [weightSum, falseWeight, zeroWeight] +example : scanPiece zeroWeight zeroWeight zeroWeight 16 2 0 0 0 = some 0 := by + simp [scanPiece, piece, pieceSize16, weightSum, zeroWeight, List.finRange_succ] +example : runScans falseWeight zeroWeight zeroWeight 16 2 + [(((0 : Fin 3), 0, 0), (0, 0))] = (none, 1) := by + have hscan : scanPiece falseWeight zeroWeight zeroWeight 16 2 0 0 0 = none := by + rw [scanPiece_eq_none_iff] + intro c _ + norm_num [weightSum, falseWeight, zeroWeight] + simp [runScans, hscan] +example : runScans zeroWeight zeroWeight zeroWeight 16 2 + [(((0 : Fin 3), 0, 0), (0, 0)), ((0, 0, 1), (1, 1))] = (some (0, 0, 0), 1) := by + simp [runScans, scanPiece, piece, pieceSize16, weightSum, zeroWeight, List.finRange_succ] + +-- The query at rank 63 ends chunk 0; rank 64 starts chunk 1. +example : (3, 15) ∈ chunkOf zeroWeight 16 3 0 0 := by + have hclass : residueClass zeroWeight 3 0 = Finset.univ := by + ext q + simp [residueClass, zeroWeight] + simp [chunkOf, queryCap16, hclass, chunk, rankIn_univ, pairIndex] +example : (4, 0) ∈ chunkOf zeroWeight 16 3 0 1 := by + have hclass : residueClass zeroWeight 3 0 = Finset.univ := by + ext q + simp [residueClass, zeroWeight] + simp [chunkOf, queryCap16, hclass, chunk, rankIn_univ, pairIndex] + +private theorem falseQuery : (0, 0) ∈ + (lopInstance falseWeight zeroWeight zeroWeight 16 2 3 (0, 0, 0)).W := by + change (0, 0) ∈ chunkOf falseWeight 16 3 0 0 + have hclass : residueClass falseWeight 3 0 = Finset.univ := by + ext q + norm_num [residueClass, falseWeight, Int.ModEq] + simp [chunkOf, queryCap16, hclass, chunk, rankIn_univ, pairIndex] + +-- A modular query really is accepted despite the absence of any exact zero triangle. +example : + (lopInstance falseWeight zeroWeight zeroWeight 16 2 3 (0, 0, 0)).HasCommonNeighbor 0 0 := by + rw [hasCommonNeighbor_lopInstance_iff falseWeight zeroWeight zeroWeight (by decide) falseQuery] + refine ⟨0, ?_, ?_⟩ + · simp [piece, pieceSize16] + · decide + +-- Supply the actual finite modular answers, rather than an abstract unconstrained answer premise. +private def modularAnswers (wAB wBC wAC : Fin 16 → Fin 16 → ℤ) + (ι : InstanceIndex 3) (q : Fin 16 × Fin 16) : Bool := + decide (∃ c ∈ piece 16 16 2 ι.2.2, weightSum wAB wBC wAC q.1 q.2 c ≡ 0 [ZMOD (3 : ℤ)]) + +private theorem modularAnswers_correct (wAB wBC wAC : Fin 16 → Fin 16 → ℤ) : + AnswersCorrect wAB wBC wAC 16 2 3 (modularAnswers wAB wBC wAC) := by + intro ι _ q hq + apply propext + simp only [modularAnswers, decide_eq_true_eq] + exact (hasCommonNeighbor_lopInstance_iff wAB wBC wAC (by decide) hq).symm + +private theorem zeroSelected : IsSelectedPrime zeroWeight zeroWeight zeroWeight 16 3 := by + refine ⟨prime3, fun q _ => ?_⟩ + simp [countZeroMod, weightSum, zeroWeight] + +-- The complete recovery theorem is inhabited at its smallest admitted n and D. +example : ∃ L, IsScanOrder zeroWeight zeroWeight zeroWeight 16 2 3 + (modularAnswers zeroWeight zeroWeight zeroWeight) L ∧ + (runScans zeroWeight zeroWeight zeroWeight 16 2 L).1 ≠ none ∧ + (runScans zeroWeight zeroWeight zeroWeight 16 2 L).2 ≤ 1 := by + obtain ⟨L, hL⟩ := exists_isScanOrder zeroWeight zeroWeight zeroWeight 16 2 3 + (modularAnswers zeroWeight zeroWeight zeroWeight) + have hc := modularAnswers_correct zeroWeight zeroWeight zeroWeight + obtain ⟨hout, hcount⟩ := recovery_correct (by decide) (by decide) (by decide) + zeroWeight zeroWeight zeroWeight zeroSelected hc hL + refine ⟨L, hL, fun hn => hout.2 hn ⟨0, 0, 0, rfl⟩, ?_⟩ + simpa [falsePositiveCount, IsFalsePositive, weightSum, zeroWeight] using hcount + +private theorem falseSelected : IsSelectedPrime falseWeight zeroWeight zeroWeight 16 3 := by + refine ⟨prime3, fun q hq => ?_⟩ + obtain ⟨hprime, -, hsq⟩ := mem_primesInRange_iff_sq.mp hq + have hq : q = 2 ∨ q = 3 := by + have := hprime.two_le + have : q < 4 := by nlinarith + omega + rcases hq with rfl | rfl <;> + norm_num [countZeroMod, weightSum, falseWeight, zeroWeight, Int.ModEq] + +example : ∃ L, IsScanOrder falseWeight zeroWeight zeroWeight 16 2 3 + (modularAnswers falseWeight zeroWeight zeroWeight) L ∧ + (runScans falseWeight zeroWeight zeroWeight 16 2 L).1 = none := by + obtain ⟨L, hL⟩ := exists_isScanOrder falseWeight zeroWeight zeroWeight 16 2 3 + (modularAnswers falseWeight zeroWeight zeroWeight) + refine ⟨L, hL, (recovery_eq_none_iff (by decide) (by decide) (by decide) + falseWeight zeroWeight zeroWeight falseSelected + (modularAnswers_correct falseWeight zeroWeight zeroWeight) hL).mpr ?_⟩ + simp [HasZeroTriangle, falseWeight, zeroWeight] diff --git a/ComplexityTests/ThreeSumHistogram.lean b/ComplexityTests/ThreeSumHistogram.lean new file mode 100644 index 0000000000..d51a5340f6 --- /dev/null +++ b/ComplexityTests/ThreeSumHistogram.lean @@ -0,0 +1,151 @@ +/- +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.ThreeSum.ChanHe.Count + +/-! # Executed histogram and residue-array regressions -/ + +namespace ComplexityTests.ThreeSumHistogram + +open Cslib.WordRAM Cslib.WordRAM.Signed +open Cslib.Algorithms.WordRAM + +private def keys (i : ℕ) : ℕ := [0, 1, 0, 0].getD i 0 + +private def input : Data 11 := + (((Data.init [0, 1, 0, 0]).set Histogram.src 1).set Histogram.dst 10).set Histogram.len 4 + |>.set Histogram.bins 3 + +/-- Clearing ignores the old contents, including nonzero words in the output block. -/ +example : (Code.eval 8 4 Histogram.program (input.write 10 77)).map + (fun yn => (yn.1.mem 10, yn.1.mem 11, yn.1.mem 12, yn.1.regs Histogram.total, yn.2)) = + some (3, 1, 0, 6, 68) := by decide + +/-- Read the same result through the generic charged theorem, checking both frames. -/ +example : ∃ y, ExecWithin 8 Histogram.program input y 68 ∧ + y.mem 10 = 3 ∧ y.mem 11 = 1 ∧ y.regs Histogram.total = 6 ∧ + y.mem 0 = 4 ∧ y.mem 13 = 0 := by + obtain ⟨y, hy, hm, ht, _⟩ := Histogram.exec (w := 8) input keys + (by decide) (by decide) (by decide) (by decide) (by decide) + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + refine ⟨y, hy, ?_, ?_, ?_, ?_, ?_⟩ + · rw [hm]; decide + · rw [hm]; decide + · rw [ht]; decide + · rw [hm]; decide + · rw [hm]; decide + +/-- Empty inputs and empty histograms pay only setup and two failed loop tests. -/ +example : (Code.eval 8 0 Histogram.program + ((input.set Histogram.len 0).set Histogram.bins 0)).map + (fun yn => (yn.1.regs Histogram.total, yn.2)) = some (0, 9) := by decide + +/-- A singleton has no collision; two equal keys at different positions create two. -/ +example : Histogram.pairs (fun _ => 0) 0 = 0 ∧ + Histogram.pairs (fun _ => 0) 1 = 0 ∧ Histogram.pairs (fun _ => 0) 2 = 2 := by decide + +/-- Different keys create no collision. -/ +example : Histogram.pairs id 4 = 0 := by decide + +private def values (i : ℕ) : ℤ := [-5, -2, 1, 4].getD i 0 + +private def signedInput : Data 16 := + (((Data.init ([-5, -2, 1, 4].map encode)).set ThreeSum.ChanHe.Residues.src 1).set + ThreeSum.ChanHe.Residues.dst 10).set ThreeSum.ChanHe.Residues.len 4 + |>.set ThreeSum.ChanHe.Residues.modulus (encode 3) + +/-- Negative dividends use Euclidean residues, and the outputs are decoded unsigned keys. -/ +example : (Code.eval 8 4 ThreeSum.ChanHe.Residues.program signedInput).map + (fun yn => (List.range 4).map fun i => yn.1.mem (10 + i)) = some [1, 1, 1, 1] := by decide + +/-- The signed-array contract has satisfiable hypotheses on the executed negative input. -/ +example : ∃ y, ExecWithin 8 ThreeSum.ChanHe.Residues.program signedInput y 96 ∧ + y.mem 10 = 1 ∧ y.mem 13 = 1 ∧ y.mem 1 = encode (-5) ∧ y.mem 14 = 0 := by + obtain ⟨y, hy, hm, _⟩ := ThreeSum.ChanHe.Residues.exec (w := 6) (M := 3) + signedInput values (by decide) (by decide) rfl + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + (by decide) (by decide) (by decide) + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + refine ⟨y, hy, ?_, ?_, ?_, ?_⟩ <;> rw [hm] <;> decide + +/-- The existing mathematical collision count agrees with the histogram's counter. -/ +example : Cslib.Algorithms.ThreeSum.ChanHe.coll {-5, -2, 1, 4} 3 = 12 ∧ + Histogram.pairs (fun i => (values i % 3).toNat) 4 = 12 := by decide + +private def countInput : Data 28 := + (((((Data.init ([-5, -2, 1, 4].map encode)).set ThreeSum.ChanHe.Count.src 1).set + ThreeSum.ChanHe.Count.scratch 10).set ThreeSum.ChanHe.Count.len 4).set + ThreeSum.ChanHe.Count.modulus (encode 3)).set ThreeSum.ChanHe.Count.dst 20 + +-- Evaluating the three renamed loops unfolds over 150 native instructions and their register maps. +set_option maxRecDepth 4096 in +/-- Execute the composed signed-input program and check its unsigned count and both arrays. -/ +example : (Code.eval 8 4 ThreeSum.ChanHe.Count.program countInput).map + (fun yn => [yn.1.mem 1, yn.1.mem 10, yn.1.mem 13, yn.1.mem 20, yn.1.mem 21, + yn.1.mem 22, yn.1.regs ThreeSum.ChanHe.Count.total]) = + some [encode (-5), 1, 1, 0, 4, 0, 12] := by decide + +/-- The full composed contract is admissible for the same negative-input example. -/ +example : ∃ y, ExecWithin 8 ThreeSum.ChanHe.Count.program countInput y 169 ∧ + y.regs ThreeSum.ChanHe.Count.total = 12 ∧ y.mem 21 = 4 ∧ y.mem 1 = encode (-5) ∧ + ∀ i, ((9 ≤ i.val ∧ i.val ≤ 12) ∨ i = ThreeSum.ChanHe.Count.dst) → + y.regs i = countInput.regs i := by + obtain ⟨y, hy, hm, ht, hframe⟩ := ThreeSum.ChanHe.Count.exec (w := 6) (M := 3) + countInput values (by decide) (by decide) rfl + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + (by decide) (by decide) (by decide) (by decide) (by decide) (by decide) + (by + intro i hi + change i < 4 at hi + have : i = 0 ∨ i = 1 ∨ i = 2 ∨ i = 3 := by omega + rcases this with rfl | rfl | rfl | rfl <;> decide) + refine ⟨y, hy, ?_, ?_, ?_, hframe⟩ + · rw [ht]; decide + · rw [hm]; decide + · rw [hm]; decide + +-- The three complete program evaluations unfold the renamed loops and register maps. +set_option maxRecDepth 4096 in +/-- The composed program handles lengths zero, one and two, including an unused modulus block. -/ +example : ((Code.eval 8 3 ThreeSum.ChanHe.Count.program + (countInput.set ThreeSum.ChanHe.Count.len 0)).map + (fun yn => yn.1.regs ThreeSum.ChanHe.Count.total)) = some 0 ∧ + ((Code.eval 8 3 ThreeSum.ChanHe.Count.program + (countInput.set ThreeSum.ChanHe.Count.len 1)).map + (fun yn => yn.1.regs ThreeSum.ChanHe.Count.total)) = some 0 ∧ + ((Code.eval 8 3 ThreeSum.ChanHe.Count.program + (countInput.set ThreeSum.ChanHe.Count.len 2)).map + (fun yn => yn.1.regs ThreeSum.ChanHe.Count.total)) = some 2 := by decide + +/-- Equal integer values violate the set-enumeration premise, even at different positions. -/ +example : ¬Set.InjOn (fun _ : ℕ => (-1 : ℤ)) (↑(Finset.range 2) : Set ℕ) := by + intro h + have := h (by decide : 0 ∈ (↑(Finset.range 2) : Set ℕ)) + (by decide : 1 ∈ (↑(Finset.range 2) : Set ℕ)) rfl + omega + +end ComplexityTests.ThreeSumHistogram diff --git a/NOTICE b/NOTICE index 95b9c8d57a..e48c2ecb6c 100644 --- a/NOTICE +++ b/NOTICE @@ -1,9 +1,11 @@ The 3SUM/Exact Triangle reductions, their finite counting lemmas and native matrix builder in -Complexity/Algorithms/ExactTriangle/{Defs,Weights,Lopsided,Hashing}.lean, +Complexity/Algorithms/ExactTriangle/{Defs,Weights,Lopsided,Hashing,Splitting,Witnesses}.lean, Complexity/Algorithms/ThreeSum/, Complexity/Algorithms/WordRAM/ThreeSum/, -Complexity/Foundations/Data/Finset/Average.lean, -Complexity/Foundations/Data/Matrix/IntegerBits.lean and +Complexity/Foundations/Data/Finset/{Average,Quotient}.lean, +Complexity/Foundations/Data/Int/Residue.lean, +Complexity/Foundations/Data/Nat/CeilDiv.lean, +Complexity/Foundations/Data/Matrix/{IntegerBits,PairChunks}.lean and Complexity/Foundations/NumberTheory/{PrimeCounting,PrimeWindow}.lean and Complexity/Foundations/RingTheory/CyclicIntRing.lean are adapted from anthropics/formal-math, 3sum-apsp, revision diff --git a/ROADMAP.md b/ROADMAP.md index a57d6f3c52..2ae3bfeb3c 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -423,10 +423,11 @@ each statement against the numbered result. Organize modules by subject, not cha matrix builder and the complete finite mathematical Chan–He reduction ([#1435]), including ordinary distinct-position inputs and explicit size bounds. The finite lopsided matrix/counting interface ([#1438]) adds padding, detection and signed digit reconstruction. Modular counting - and the selected-prime false-positive bound ([#1442]) are proved mathematically; it still needs - residue/chunk/piece splitting and witness recovery, tensor/sparse-query programs, native - Chan–He preprocessing, the reduction host, the fast Exact Triangle solver and full input/time/width/ - asymptotic composition, then the improved Section 4 route to exponent 1.9992; size-based width + and the selected-prime false-positive bound ([#1442]), finite residue/chunk/piece splitting and + witness recovery ([#1447]) are proved mathematically. It still needs tensor/sparse-query programs, + native Chan–He preprocessing, the reduction host, the fast Exact Triangle solver and full + input/time/width/asymptotic composition, then the improved Section 4 route to exponent 1.9992; + size-based width selection with an input-fit guarantee and explicit overhead; expected-time rules for unbounded almost-surely terminating structured loops; random-tape equivalence and BPP/RP/ZPP simulation bridges for the randomized word RAM ([#1124]); transfer of cell-probe lower bounds; and the @@ -569,4 +570,5 @@ not planned, not as proved. [#1435]: https://github.com/cslib-community/complexitylib/issues/1435 [#1438]: https://github.com/cslib-community/complexitylib/issues/1438 [#1442]: https://github.com/cslib-community/complexitylib/issues/1442 +[#1447]: https://github.com/cslib-community/complexitylib/issues/1447 [#1449]: https://github.com/cslib-community/complexitylib/issues/1449 diff --git a/UPSTREAM.md b/UPSTREAM.md index 30a7a8947f..1023dc2cb9 100644 --- a/UPSTREAM.md +++ b/UPSTREAM.md @@ -189,11 +189,14 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he all-window prime-count bound lives in `Foundations/NumberTheory/PrimeWindow.lean`, retaining the source credit to Mathlib, Patrick Stevens and Bolton Bailey. These are mathematical statements, not native instructions or an assumed solver. See - [#1442](https://github.com/cslib-community/complexitylib/issues/1442). Residue/chunk/piece splitting, - witness recovery and tensor/sparse query programs remain in #1394, alongside native composition. - The existing signed and unsigned APIs were reused; no source machine or separate run semantics - were copied. The full reduction host, native Chan–He preprocessing, fast Exact Triangle and full - input/time/width/asymptotic composition remain in + [#1442](https://github.com/cslib-community/complexitylib/issues/1442). + `Algorithms/ExactTriangle/{Splitting,Witnesses}.lean` adapts the complete finite splitting and + witness recovery component ([#1447](https://github.com/cslib-community/complexitylib/issues/1447)), + using the existing matrices, lopsided graph and `Set.EqOn` answer API. Computable integer parameters + retain the source's exact real floor/ceiling specifications through checked bridges. The existing + signed and unsigned APIs were reused; no source machine or separate run semantics were copied. + Tensor/sparse query programs, the full reduction host, native Chan–He preprocessing, fast Exact + Triangle and full input/time/width/asymptotic composition remain in [#1394](https://github.com/cslib-community/complexitylib/issues/1394), including the improved Section 4 route needed for the published exponent 1.9992. Anthropic's copyright and attribution are retained; the native proofs and representation choices are local.