diff --git a/Complexity.lean b/Complexity.lean index e2b1bd7e08..bf3ca017d2 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -1,7 +1,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.Weights public import Complexity.Algorithms.Lean.FusionTree.Capacity public import Complexity.Algorithms.Lean.FusionTree.Compression public import Complexity.Algorithms.Lean.FusionTree.Correction @@ -48,6 +50,9 @@ public import Complexity.Algorithms.ThreeSum.ChanHe.Sizes public import Complexity.Algorithms.ThreeSum.Convolution public import Complexity.Algorithms.ThreeSum.Defs public import Complexity.Algorithms.ThreeSum.ThreeSets +public import Complexity.Algorithms.WordRAM.Deduplicate.Count +public import Complexity.Algorithms.WordRAM.Deduplicate.Model +public import Complexity.Algorithms.WordRAM.Deduplicate.Program public import Complexity.Algorithms.WordRAM.Dijkstra.Fibonacci public import Complexity.Algorithms.WordRAM.Dijkstra.FibonacciStatement public import Complexity.Algorithms.WordRAM.Dijkstra.GenericPQ @@ -112,6 +117,11 @@ public import Complexity.Algorithms.WordRAM.Queue.Dynamic public import Complexity.Algorithms.WordRAM.Queue.Model public import Complexity.Algorithms.WordRAM.Quickselect public import Complexity.Algorithms.WordRAM.Quicksort +public import Complexity.Algorithms.WordRAM.RadixSort +public import Complexity.Algorithms.WordRAM.RadixSort.Deduplicate +public import Complexity.Algorithms.WordRAM.RadixSort.Model +public import Complexity.Algorithms.WordRAM.RadixSort.Program +public import Complexity.Algorithms.WordRAM.RadixSort.Scan public import Complexity.Algorithms.WordRAM.Stack public import Complexity.Algorithms.WordRAM.Stack.Controller public import Complexity.Algorithms.WordRAM.ThreeSum.Cell @@ -1622,6 +1632,7 @@ public import Complexity.Computability.WordRAM.ExactNotation public import Complexity.Computability.WordRAM.Function public import Complexity.Computability.WordRAM.Generic public import Complexity.Computability.WordRAM.Heap +public import Complexity.Computability.WordRAM.ListArray public import Complexity.Computability.WordRAM.MapOp public import Complexity.Computability.WordRAM.Multiword public import Complexity.Computability.WordRAM.Multiword.Arithmetic @@ -2219,6 +2230,7 @@ public import Complexity.Foundations.NumberTheory.AKS.Introspective public import Complexity.Foundations.NumberTheory.AKS.SmallR public import Complexity.Foundations.NumberTheory.Harmonic public import Complexity.Foundations.NumberTheory.PrimeCounting +public import Complexity.Foundations.NumberTheory.PrimeWindow public import Complexity.Foundations.Order.BooleanSubalgebra public import Complexity.Foundations.Order.Changes public import Complexity.Foundations.Order.UpperLower @@ -2229,6 +2241,7 @@ public import Complexity.Foundations.Probability.SecondMoment public import Complexity.Foundations.Quantum.Measurement public import Complexity.Foundations.Relation.InductiveCounting public import Complexity.Foundations.Relation.RelatesWithinSteps +public import Complexity.Foundations.RingTheory.CyclicIntRing public import Complexity.Foundations.Semantics.Stateful.Amortized public import Complexity.Foundations.Semantics.Stateful.Basic public import Complexity.Foundations.Semantics.Stateful.Composition diff --git a/Complexity/Algorithms/ExactTriangle/Hashing.lean b/Complexity/Algorithms/ExactTriangle/Hashing.lean new file mode 100644 index 0000000000..fbb2834b6e --- /dev/null +++ b/Complexity/Algorithms/ExactTriangle/Hashing.lean @@ -0,0 +1,508 @@ +/- +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.Weights +public import Complexity.Foundations.NumberTheory.PrimeWindow +public import Complexity.Foundations.RingTheory.CyclicIntRing +public import Mathlib.Data.Matrix.Mul + +/-! +# Modular counts and prime selection for integer Exact Triangle + +Let `S(a,b,c)` denote `weightSum`, `F(p)` the false-positive count and `Z₀` the zero count. A +minimum-count prime in `[sqrt D / 2, sqrt D)` also minimises false positives. For weights +bounded by `n^κ`, with `16 ≤ D ≤ n` and `1 ≤ κ`, a universal constant bounds its false positives +by `C * κ * n^3 * log n / sqrt D`. The prime count used here is proved for every admitted `D`. + +## Main definitions + +* `countZeroMod`, `falsePositiveCount`, `zeroCount`: counts of finite triples. +* `IsSelectedPrime`: a minimum-count prime in the specified window. +* `matP`, `matQ`: cyclic polynomial matrices encoding edge residues. + +## Main results + +* `countZeroMod_eq`: modular count is the sum of false positives and true zeros. +* `coeff_matP_mul_matQ`, `falsePositiveCount_add_zeroCount_eq_sum_coeff`: coefficient + identities. +* `exists_isSelectedPrime`: the minimum exists. +* `Hashing.falsePositiveCount_le_const_mul`: the uniform false-positive bound. + +This is a mathematical reduction step. No native program or charged running-time bound is +asserted. The supplied three matrices retain the existing `HasZeroTriangle` semantics. + +## References + +The modular reduction is part of [AlmanWilliams2026, Section 3]. +Direct adaptation of [Anthropic's +Hashing.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Hashing.lean), +with the minimal count/selection definitions extracted from its `PaperStatements.lean`. +The declaration docstrings describe the Lean definitions and statements in local notation. +The source pin is +`anthropics/formal-math@16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.ExactTriangle + +open Classical in +/-- Prime numbers in the half-open real interval `[sqrt D / 2, sqrt D)`. +-/ +noncomputable def primesInRange (D : ℕ) : Finset ℕ := + (Finset.range D).filter fun p => p.Prime ∧ Real.sqrt D / 2 ≤ (p : ℝ) ∧ (p : ℝ) < Real.sqrt D + +/-- Count vertex triples whose three edge weights have sum congruent to zero modulo `p`. +-/ +def countZeroMod {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (p : ℕ) : ℕ := + (Finset.univ.filter fun t : Fin n × Fin n × Fin n => + weightSum wAB wBC wAC t.1 t.2.1 t.2.2 ≡ 0 [ZMOD (p : ℤ)]).card + +/-- A prime in `primesInRange D` whose modular-zero count is minimal among all primes in that +window. Tied minima all satisfy this predicate. +-/ +def IsSelectedPrime {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) (D p : ℕ) : Prop := + p ∈ primesInRange D ∧ ∀ q ∈ primesInRange D, + countZeroMod wAB wBC wAC p ≤ countZeroMod wAB wBC wAC q + +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 -/ + +/-- Window membership is equivalent to primality and the two real endpoint inequalities. The +finite enumeration bound in `primesInRange` excludes no member of the interval. +-/ +theorem mem_primesInRange {D p : ℕ} : + p ∈ primesInRange D ↔ p.Prime ∧ Real.sqrt D / 2 ≤ (p : ℝ) ∧ (p : ℝ) < Real.sqrt D := by + unfold primesInRange + rw [Finset.mem_filter, Finset.mem_range] + refine ⟨fun h => h.2, fun h => ⟨?_, h⟩⟩ + -- From `p < √D` we get `p ≤ p² < D`. + have hsq : p ^ 2 < D := by exact_mod_cast (Real.lt_sqrt (Nat.cast_nonneg p)).1 h.2.2 + exact (Nat.le_self_pow two_ne_zero p).trans_lt hsq + +/-- Characterize window membership using primality and integer-square inequalities, including the +closed lower and open upper endpoints. +-/ +theorem mem_primesInRange_iff_sq {D p : ℕ} : + p ∈ primesInRange D ↔ p.Prime ∧ D ≤ 4 * p ^ 2 ∧ p ^ 2 < D := by + rw [mem_primesInRange] + have hp0 : (0 : ℝ) ≤ p := Nat.cast_nonneg p + have hsqrt0 : 0 ≤ Real.sqrt (D : ℝ) := Real.sqrt_nonneg _ + have hsqrt : (Real.sqrt (D : ℝ)) ^ 2 = D := Real.sq_sqrt (Nat.cast_nonneg D) + constructor + · rintro ⟨hp, hlo, hhi⟩ + refine ⟨hp, ?_, ?_⟩ + · have h : (D : ℝ) ≤ 4 * (p : ℝ) ^ 2 := by nlinarith + exact_mod_cast h + · exact_mod_cast (Real.lt_sqrt hp0).1 hhi + · rintro ⟨hp, hlo, hhi⟩ + refine ⟨hp, ?_, (Real.lt_sqrt hp0).2 (by exact_mod_cast hhi)⟩ + have h : (D : ℝ) ≤ 4 * (p : ℝ) ^ 2 := by exact_mod_cast hlo + nlinarith + +/-- Enumerate the real prime window using a decidable filter of the integers below `D`. +-/ +theorem primesInRange_eq_filter (D : ℕ) : + primesInRange D = (Finset.range D).filter fun p => p.Prime ∧ D ≤ 4 * p ^ 2 ∧ p ^ 2 < D := by + ext p + rw [mem_primesInRange_iff_sq, Finset.mem_filter, Finset.mem_range] + exact ⟨fun h => ⟨(Nat.le_self_pow two_ne_zero p).trans_lt h.2.2, h⟩, fun h => h.2⟩ + +/-- A number `r < p` is the sum of the residues of `u` and `v` in `{0, …, p − 1}`, reduced modulo +`p`, exactly if `u + v ≡ r`. +-/ +private theorem eq_add_toNat_emod_iff {p r : ℕ} (hp : p ≠ 0) (hr : r < p) (u v : ℤ) : + r = ((u % (p : ℤ)).toNat + (v % (p : ℤ)).toNat) % p ↔ u + v ≡ (r : ℤ) [ZMOD (p : ℤ)] := by + 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), + ← Int.add_emod] + +/-! +### Counting the triples with `S(a,b,c) ≡ 0 (mod p)` -/ + +section Counting + +variable {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) + +/-- A vertex triple has a nonzero integer weight sum that is congruent to zero modulo `p`. +-/ +abbrev IsFalsePositive (p : ℕ) (t : Fin n × Fin n × Fin n) : Prop := + weightSum wAB wBC wAC t.1 t.2.1 t.2.2 ≠ 0 ∧ (p : ℤ) ∣ weightSum wAB wBC wAC t.1 t.2.1 t.2.2 + +/-- Count vertex triples satisfying `IsFalsePositive` for the supplied three weight matrices and +modulus. +-/ +def falsePositiveCount (p : ℕ) : ℕ := + (Finset.univ.filter fun t : Fin n × Fin n × Fin n => IsFalsePositive wAB wBC wAC p t).card + +/-- Count vertex triples whose three supplied edge weights sum to zero as an integer. +-/ +def zeroCount : ℕ := + (Finset.univ.filter fun t : Fin n × Fin n × Fin n => weightSum wAB wBC wAC t.1 t.2.1 t.2.2 = + 0).card + +/-- The zero-triple count is positive exactly when the supplied matrices satisfy +`HasZeroTriangle`. +-/ +theorem zeroCount_pos_iff : 0 < zeroCount wAB wBC wAC ↔ HasZeroTriangle wAB wBC wAC := by + classical + simp only [zeroCount, Finset.card_pos, HasZeroTriangle] + constructor + · rintro ⟨⟨a, b, c⟩, h⟩ + exact ⟨a, b, c, (Finset.mem_filter.mp h).2⟩ + · rintro ⟨a, b, c, h⟩ + exact ⟨(a, b, c), Finset.mem_filter.mpr ⟨Finset.mem_univ _, h⟩⟩ + +/-- Partition the modular-zero triples into false positives and exact zero-sum triples. The +exact-zero count is independent of the modulus. +-/ +theorem countZeroMod_eq (p : ℕ) : (countZeroMod wAB wBC wAC) p = (falsePositiveCount wAB wBC + wAC) p + (zeroCount wAB wBC wAC) := by + classical + simp only [countZeroMod, falsePositiveCount, zeroCount, Finset.card_filter, ← + Finset.sum_add_distrib] + -- A triple with `p ∣ S(a,b,c)` is a zero triangle or a false positive, and not both. + refine Finset.sum_congr rfl fun t _ => ?_ + by_cases h0 : weightSum wAB wBC wAC t.1 t.2.1 t.2.2 = 0 <;> + simp [IsFalsePositive, Int.modEq_zero_iff_dvd, h0] + +/-- Encode the A–C weight at `(a,c)` as a power of the cyclic polynomial generator, with exponent +equal to its nonnegative residue modulo `p`. +-/ +noncomputable def matP (p : ℕ) : Matrix (Fin n) (Fin n) (CyclicIntRing p) := + Matrix.of fun a c => CyclicIntRing.x p ^ (wAC a c % (p : ℤ)).toNat + +/-- Encode the B–C weight at `(b,c)` in matrix position `(c,b)`. This orientation makes the +product sum over the middle index `c`. +-/ +noncomputable def matQ (p : ℕ) : Matrix (Fin n) (Fin n) (CyclicIntRing p) := + Matrix.of fun c b => CyclicIntRing.x p ^ (wBC b c % (p : ℤ)).toNat + +/-- For a nonzero modulus and a residue `r < p`, the coefficient of degree `r` in the cyclic +matrix product counts middle vertices whose two incident weights sum to `r` modulo `p`. +-/ +theorem coeff_matP_mul_matQ {p : ℕ} (hp : p ≠ 0) (a b : Fin n) (r : ℕ) (hr : r < p) : + CyclicIntRing.coeff hp r ((matP wAC p * matQ wBC p) a b) = + ((Finset.univ.filter fun c : Fin n => + wAC a c + wBC b c ≡ (r : ℤ) [ZMOD (p : ℤ)]).card : ℤ) := by + classical + -- `(PQ)[a,b]` is the sum over `c` of `x^(w(a,c) mod p + w(b,c) mod p)`, and the coefficient of + -- `x^r` in each term is 1 or 0. + rw [Matrix.mul_apply] + simp only [matP, matQ, Matrix.of_apply, ← pow_add] + rw [map_sum] + simp only [CyclicIntRing.coeff_x_pow] + rw [Finset.card_filter] + push_cast + refine Finset.sum_congr rfl fun c _ => ?_ + simp only [eq_add_toNat_emod_iff hp hr] + +/-- Recover the total modular-zero count from the cyclic matrix product by summing, over `(a,b)`, +the coefficient at the residue of `-wAB a b`. +-/ +theorem falsePositiveCount_add_zeroCount_eq_sum_coeff {p : ℕ} (hp : p ≠ 0) : + (((falsePositiveCount wAB wBC wAC) p + (zeroCount wAB wBC wAC) : ℕ) : ℤ) = + ∑ a : Fin n, ∑ b : Fin n, + CyclicIntRing.coeff hp ((-wAB a b) % (p : ℤ)).toNat ((matP wAC p * matQ wBC p) a b) := by + classical + rw [← countZeroMod_eq] + unfold countZeroMod + rw [Finset.card_filter, Fintype.sum_prod_type] + push_cast + 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) _), + Finset.card_filter] + push_cast + refine Finset.sum_congr rfl fun c _ => ?_ + -- `S(a,b,c) ≡ 0` says the same as `w(a,c) + w(b,c) ≡ -w(a,b)`, and `-w(a,b)` is congruent to its + -- residue. + have hres : (-wAB a b) % (p : ℤ) ≡ -wAB a b [ZMOD (p : ℤ)] := Int.mod_modEq _ _ + have hiff : weightSum wAB wBC wAC a b c ≡ 0 [ZMOD (p : ℤ)] ↔ + wAC a c + wBC b c ≡ ((((-wAB a b) % (p : ℤ)).toNat : ℕ) : ℤ) [ZMOD (p : ℤ)] := by + rw [hcast, show wAC a c + wBC b c ≡ (-wAB a b) % (p : ℤ) [ZMOD (p : ℤ)] ↔ + wAC a c + wBC b c ≡ -wAB a b [ZMOD (p : ℤ)] from + ⟨fun h => h.trans hres, fun h => h.trans hres.symm⟩, + Int.modEq_iff_dvd, Int.modEq_iff_dvd, + show 0 - weightSum wAB wBC wAC a b c = -wAB a b - (wAC a c + wBC b c) by + simp only [weightSum]; ring] + simp only [hiff] + +end Counting + +/-- For `D ≥ 16`, the number of primes in `primesInRange D` is strictly less than `sqrt D`. +-/ +theorem card_primesInRange_lt (D : ℕ) (hD : 16 ≤ D) : + ((primesInRange D).card : ℝ) < Real.sqrt D := by + have hpos : 0 < Real.sqrt D := by linarith [four_le_sqrt hD] + -- The primes in the range are among the numbers `1, …, ⌈√D⌉ - 1`. + have hsub : primesInRange D ⊆ Finset.Ico 1 ⌈Real.sqrt D⌉₊ := by + intro p hp + obtain ⟨hprime, -, hlt⟩ := mem_primesInRange.1 hp + exact Finset.mem_Ico.2 ⟨hprime.one_lt.le, Nat.lt_ceil.2 hlt⟩ + have hcard : (primesInRange D).card + 1 ≤ ⌈Real.sqrt D⌉₊ := by + have hle := Finset.card_le_card hsub + have hone : 1 ≤ ⌈Real.sqrt D⌉₊ := Nat.one_le_ceil_iff.2 hpos + rw [Nat.card_Ico] at hle + omega + have hcast : ((primesInRange D).card : ℝ) + 1 ≤ (⌈Real.sqrt D⌉₊ : ℝ) := by exact_mod_cast hcard + linarith [Nat.ceil_lt_add_one hpos.le] + +/-! +### Selecting the prime -/ + +section Selection + +variable {n D p : ℕ} {κ : ℝ} (wAB wBC wAC : Fin n → Fin n → ℤ) + +/-- For every `D ≥ 16` and every three supplied weight matrices, the prime window contains a +minimum of the modular-zero count. +-/ +theorem exists_isSelectedPrime (D : ℕ) (hD : 16 ≤ D) : ∃ p, IsSelectedPrime wAB wBC wAC D p := by + obtain ⟨q, hq⟩ := + Nat.exists_prime_half_le_and_lt (Real.sqrt D) (four_le_sqrt hD) + exact Finset.exists_min_image (primesInRange D) (countZeroMod wAB wBC wAC) ⟨q, + mem_primesInRange.2 hq⟩ + +/-- A selected minimum of the modular-zero count also minimizes the false-positive count over the +same prime window. +-/ +theorem IsSelectedPrime.falsePositiveCount_le {wAB wBC wAC : Fin n → Fin n → ℤ} (hp : + IsSelectedPrime wAB wBC wAC D p) : + ∀ q ∈ primesInRange D, (falsePositiveCount wAB wBC wAC) p ≤ (falsePositiveCount wAB wBC wAC) + q := by + intro q hq + have h := hp.2 q hq + rw [countZeroMod_eq, countZeroMod_eq] at h + omega + +/-! +### The bound on the number of false positives of the selected prime -/ + +open Classical in +/-- For a fixed nonzero-sum triple and edge weights bounded by `n^κ`, bound the number of window +primes accepting that triple by the base-`sqrt D / 2` logarithm of `3*n^κ`. +-/ +theorem card_falsePositive_primes_le (hD : 16 ≤ D) (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) + (t : Fin n × Fin n × Fin n) (ht : weightSum wAB wBC wAC t.1 t.2.1 t.2.2 ≠ 0) : + (((primesInRange D).filter fun p => IsFalsePositive wAB wBC wAC p t).card : ℝ) + ≤ Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := by + have hsqrt := four_le_sqrt hD + set divisors := (primesInRange D).filter fun p => IsFalsePositive wAB wBC wAC p t + have hmem : ∀ p ∈ divisors, + p.Prime ∧ Real.sqrt D / 2 ≤ (p : ℝ) ∧ (p : ℤ) ∣ weightSum wAB wBC wAC t.1 t.2.1 t.2.2 := by + intro p hp + obtain ⟨hrange, -, hdvd⟩ := Finset.mem_filter.1 hp + obtain ⟨hprime, hge, -⟩ := mem_primesInRange.1 hrange + exact ⟨hprime, hge, hdvd⟩ + -- The product of these primes divides `S(a,b,c)`, so it is at most `|S(a,b,c)| ≤ 3n^κ`. + have hdvd : ∏ p ∈ divisors, p ∣ (weightSum wAB wBC wAC t.1 t.2.1 t.2.2).natAbs := + Finset.prod_primes_dvd _ (fun p hp => (hmem p hp).1.prime) + fun p hp => Int.natCast_dvd.1 (hmem p hp).2.2 + have hprod : ∏ p ∈ divisors, (p : ℝ) ≤ 3 * (n : ℝ) ^ κ := + calc ∏ p ∈ divisors, (p : ℝ) = ((∏ p ∈ divisors, p : ℕ) : ℝ) := (Nat.cast_prod _ _).symm + _ ≤ ((weightSum wAB wBC wAC t.1 t.2.1 t.2.2).natAbs : ℝ) := by + exact_mod_cast Nat.le_of_dvd (Int.natAbs_pos.2 ht) hdvd + _ = ((|weightSum wAB wBC wAC t.1 t.2.1 t.2.2| : ℤ) : ℝ) := Nat.cast_natAbs _ + _ ≤ 3 * (n : ℝ) ^ κ := (abs_weightSum_le wAB wBC wAC) hT _ _ _ + -- Each of these primes is at least `√D/2`. + have hpow : (Real.sqrt D / 2) ^ divisors.card ≤ ∏ p ∈ divisors, (p : ℝ) := by + rw [← Finset.prod_const] + exact Finset.prod_le_prod₀ (fun _ _ => by linarith) fun p hp => (hmem p hp).2.1 + have hY : 0 < 3 * (n : ℝ) ^ κ := lt_of_lt_of_le (by positivity) (hpow.trans hprod) + rw [Real.le_logb_iff_rpow_le (by linarith) hY, Real.rpow_natCast] + exact hpow.trans hprod + +/-- Sum the per-triple prime-divisibility bound over all `n^3` triples. For `D ≥ 16`, `κ ≥ 1` and +edge weights bounded by `n^κ`, the sum of false-positive counts is at most `n^3 * logb (sqrt D / +2) (3*n^κ)`. +-/ +theorem sum_falsePositiveCount_le (hD : 16 ≤ D) (hκ : 1 ≤ κ) (hT : WeightsBoundedBy wAB wBC wAC + ((n : ℝ) ^ κ)) : + ((∑ p ∈ primesInRange D, (falsePositiveCount wAB wBC wAC) p : ℕ) : ℝ) + ≤ (n : ℝ) ^ 3 * Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := by + classical + have hsqrt := four_le_sqrt hD + -- Count the pairs (prime, false positive of that prime) triple by triple. + have hswap : ∑ p ∈ primesInRange D, (falsePositiveCount wAB wBC wAC) p = ∑ t : Fin n × Fin n × + Fin n, + ((primesInRange D).filter fun p => IsFalsePositive wAB wBC wAC p t).card := by + simp only [falsePositiveCount, Finset.card_filter] + exact Finset.sum_comm + have hterm : ∀ t : Fin n × Fin n × Fin n, + (((primesInRange D).filter fun p => IsFalsePositive wAB wBC wAC p t).card : ℝ) + ≤ Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := by + intro t + by_cases ht : weightSum wAB wBC wAC t.1 t.2.1 t.2.2 = 0 + · -- A zero triangle is a false positive of no prime. + rw [Finset.filter_eq_empty_iff.2 fun p _ h => h.1 ht, Finset.card_empty, Nat.cast_zero] + have hn : (1 : ℝ) ≤ n := by exact_mod_cast Fin.pos t.1 + have hpow : (1 : ℝ) ≤ (n : ℝ) ^ κ := Real.one_le_rpow hn (by linarith) + exact Real.logb_nonneg (by linarith) (by linarith) + · -- Otherwise this is the bound of `card_falsePositive_primes_le` for this triple. + convert (card_falsePositive_primes_le wAB wBC wAC) hD hT t ht + rw [hswap, Nat.cast_sum] + calc ∑ t : Fin n × Fin n × Fin n, + (((primesInRange D).filter fun p => IsFalsePositive wAB wBC wAC p t).card : ℝ) + ≤ ∑ _t : Fin n × Fin n × Fin n, Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := + Finset.sum_le_sum fun t _ => hterm t + _ = (n : ℝ) ^ 3 * Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := by + rw [Finset.sum_const, Finset.card_univ, Fintype.card_prod, Fintype.card_prod, + Fintype.card_fin, nsmul_eq_mul] + push_cast + ring + +end Selection + +/-- There is a universal positive `c` such that every window with `D ≥ 16` contains at least `c * +sqrt D / log D` primes. The proof uses the all-window counting form of Bertrand’s postulate from +`PrimeWindow`, with no prime-number-theorem assumption. +-/ +theorem exists_le_card_primesInRange : + ∃ c : ℝ, 0 < c ∧ ∀ D : ℕ, 16 ≤ D → + c * (Real.sqrt D / Real.log D) ≤ ((primesInRange D).card : ℝ) := by + obtain ⟨c, hc, hcount⟩ := Nat.exists_forall_mul_sqrt_div_log_le_card_primes + exact ⟨c, hc, fun D hD => hcount D hD _ + fun _ hp hge hlt => mem_primesInRange.2 ⟨hp, hge, hlt⟩⟩ + +/-- Multiply the selected prime’s minimum false-positive count by the number of window primes to +obtain a bound by the sum of their false-positive counts. +-/ +theorem IsSelectedPrime.falsePositiveCount_mul_card_le_sum {n D p : ℕ} {wAB wBC wAC : Fin n → + Fin n → ℤ} + (hp : IsSelectedPrime wAB wBC wAC D p) : + (falsePositiveCount wAB wBC wAC) p * (primesInRange D).card ≤ ∑ q ∈ primesInRange D, + (falsePositiveCount wAB wBC wAC) q := by + rw [mul_comm, ← smul_eq_mul] + exact Finset.card_nsmul_le_sum _ _ _ hp.falsePositiveCount_le + +/-- For `D ≥ 16`, compare `log D` with four times `log (sqrt D / 2)`. +-/ +theorem log_le_four_mul_log_sqrt_div_two (D : ℕ) (hD : 16 ≤ D) : + Real.log D ≤ 4 * Real.log (Real.sqrt D / 2) := by + have hD16 : (16 : ℝ) ≤ D := by exact_mod_cast hD + have hsqrt : 0 < Real.sqrt D := by linarith [four_le_sqrt hD] + -- `log(√D/2) = (log D)/2 - log 2`, and `4 log 2 = log 16 ≤ log D`. + have hlog16 : 4 * Real.log 2 ≤ Real.log D := by + have h := Real.log_le_log (by norm_num) (show (2 : ℝ) ^ 4 ≤ D by linarith) + rw [Real.log_pow] at h + exact_mod_cast h + rw [Real.log_div hsqrt.ne' (by norm_num), Real.log_sqrt (by linarith)] + linarith [hlog16] + +/-- For `16 ≤ D ≤ n`, `κ ≥ 1` and edge weights bounded by `n^κ`, a positive pointwise prime-count +lower bound `c * sqrt D / log D` gives every selected prime at most `(4/c) * n^3 * log (3*n^κ) / +sqrt D` false positives. +-/ +theorem falsePositiveCount_le_of_le_card_primesInRange {n D p : ℕ} {κ c : ℝ} + (wAB wBC wAC : Fin n → Fin n → ℤ) (hc : 0 < c) + (hcard : c * (Real.sqrt D / Real.log D) ≤ ((primesInRange D).card : ℝ)) (hD : 16 ≤ D) + (hDn : D ≤ n) (hκ : 1 ≤ κ) (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) (hp : + IsSelectedPrime wAB wBC wAC D p) : + ((falsePositiveCount wAB wBC wAC) p : ℝ) ≤ 4 / c * ((n : ℝ) ^ 3 * Real.log (3 * (n : ℝ) ^ κ) + / Real.sqrt D) := by + -- The quantities of the paper's sentence: `r = √D`, `Y = 3n^κ`, and the number of primes in the + -- range. + set r := Real.sqrt D + set Y := 3 * (n : ℝ) ^ κ with hY + set numPrimes := (primesInRange D).card + have hr : 4 ≤ r := four_le_sqrt hD + have hn : (16 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn + have hlogD_pos : 0 < Real.log D := Real.log_pos (by exact_mod_cast (by omega : 1 < D)) + have hlogr_pos : 0 < Real.log (r / 2) := Real.log_pos (by linarith) + have hlogY_nonneg : 0 ≤ Real.log Y := by + have : 1 ≤ (n : ℝ) ^ κ := Real.one_le_rpow (by linarith) (by linarith) + exact Real.log_nonneg (by rw [hY]; linarith) + -- The selected prime is at most average, after multiplying by the number of primes. + -- False positives of all primes in the + -- range add up to at most `n³ log_{r/2} Y`. + have average : ((falsePositiveCount wAB wBC wAC) p : ℝ) * numPrimes ≤ (n : ℝ) ^ 3 * Real.logb + (r / 2) Y := + le_trans (by exact_mod_cast hp.falsePositiveCount_mul_card_le_sum) + ((sum_falsePositiveCount_le wAB wBC wAC) hD hκ hT) + -- The base of the logarithm changes: `log_{r/2} Y ≤ 4 log Y / log D`. + have hlogD_le : Real.log D ≤ 4 * Real.log (r / 2) := log_le_four_mul_log_sqrt_div_two D hD + have changeBase : Real.logb (r / 2) Y ≤ 4 * Real.log Y / Real.log D := by + rw [Real.logb, div_le_div_iff₀ hlogr_pos hlogD_pos] + linarith [mul_le_mul_of_nonneg_left hlogD_le hlogY_nonneg] + rw [show 4 / c * ((n : ℝ) ^ 3 * Real.log Y / r) = 4 * (n : ℝ) ^ 3 * Real.log Y / (c * r) by + field_simp, le_div_iff₀ (by positivity)] + calc ((falsePositiveCount wAB wBC wAC) p : ℝ) * (c * r) = (falsePositiveCount wAB wBC wAC) p * + (c * (r / Real.log D)) * Real.log D := by field_simp + _ ≤ (falsePositiveCount wAB wBC wAC) p * numPrimes * Real.log D := by gcongr + _ ≤ (n : ℝ) ^ 3 * Real.logb (r / 2) Y * Real.log D := by gcongr + _ ≤ (n : ℝ) ^ 3 * (4 * Real.log Y / Real.log D) * Real.log D := by gcongr + _ = 4 * (n : ℝ) ^ 3 * Real.log Y := by field_simp + +/-- There is a universal `C ≥ 1` such that, for `16 ≤ D ≤ n`, `κ ≥ 1` and edge weights bounded by +`n^κ`, every selected prime has at most `C * κ * n^3 * log n / sqrt D` false positives. +-/ +theorem exists_falsePositiveCount_le : + ∃ C : ℝ, 1 ≤ C ∧ ∀ {n D p : ℕ} {κ : ℝ}, 16 ≤ D → D ≤ n → 1 ≤ κ → ∀ wAB wBC wAC : Fin n → Fin + n → ℤ, + WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ) → IsSelectedPrime wAB wBC wAC D p → + ((falsePositiveCount wAB wBC wAC) p : ℝ) ≤ C * (κ * (n : ℝ) ^ 3 * Real.log n / Real.sqrt + D) := by + obtain ⟨c, hc, hcard⟩ := exists_le_card_primesInRange + refine ⟨max 1 (8 / c), le_max_left _ _, fun {n D p κ} hD hDn hκ wAB wBC wAC hT hp => ?_⟩ + have hn : (16 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn + have hlogn : 0 ≤ Real.log n := Real.log_nonneg (by linarith) + -- `log(3n^κ) = log 3 + κ log n ≤ 2κ log n`, because `3 ≤ n` and `κ ≥ 1`. + have hlog : Real.log (3 * (n : ℝ) ^ κ) ≤ 2 * κ * Real.log n := by + have hlog3 : Real.log 3 ≤ κ * Real.log n := + (Real.log_le_log (by norm_num) (by linarith)).trans (le_mul_of_one_le_left hlogn hκ) + rw [Real.log_mul (by norm_num) (by positivity), Real.log_rpow (by linarith)] + linarith [hlog3] + calc ((falsePositiveCount wAB wBC wAC) p : ℝ) ≤ 4 / c * ((n : ℝ) ^ 3 * Real.log (3 * (n : ℝ) ^ + κ) / Real.sqrt D) := + (falsePositiveCount_le_of_le_card_primesInRange wAB wBC wAC) hc (hcard D hD) hD hDn hκ hT hp + _ ≤ 4 / c * ((n : ℝ) ^ 3 * (2 * κ * Real.log n) / Real.sqrt D) := by gcongr + _ = 8 / c * (κ * (n : ℝ) ^ 3 * Real.log n / Real.sqrt D) := by ring + _ ≤ max 1 (8 / c) * (κ * (n : ℝ) ^ 3 * Real.log n / Real.sqrt D) := by + gcongr + exact le_max_right _ _ + +/-- A universal real constant witnessing `exists_falsePositiveCount_le`. It depends on none of the +input dimensions, weight matrices, exponent or selected prime, and is at least one. +-/ +noncomputable def Hashing.falsePositiveConst : ℝ := Classical.choose exists_falsePositiveCount_le + +/-- The universal false-positive constant `Hashing.falsePositiveConst` is at least one. +-/ +theorem Hashing.one_le_falsePositiveConst : 1 ≤ Hashing.falsePositiveConst := + (Classical.choose_spec exists_falsePositiveCount_le).1 + +/-- For `16 ≤ D ≤ n`, `κ ≥ 1` and edge weights bounded by `n^κ`, bound every selected prime’s +false-positive count by `Hashing.falsePositiveConst * κ * n^3 * log n / sqrt D`. +-/ +theorem Hashing.falsePositiveCount_le_const_mul {n D p : ℕ} {κ : ℝ} (hD : 16 ≤ D) (hDn : D ≤ n) + (hκ : 1 ≤ κ) (wAB wBC wAC : Fin n → Fin n → ℤ) (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) + (hp : IsSelectedPrime wAB wBC wAC D p) : + ((falsePositiveCount wAB wBC wAC) p : ℝ) ≤ Hashing.falsePositiveConst * (κ * (n : ℝ) ^ 3 * + Real.log n / Real.sqrt D) := + (Classical.choose_spec exists_falsePositiveCount_le).2 hD hDn hκ wAB wBC wAC hT hp + +end Cslib.Algorithms.ExactTriangle diff --git a/Complexity/Algorithms/ExactTriangle/Weights.lean b/Complexity/Algorithms/ExactTriangle/Weights.lean new file mode 100644 index 0000000000..d358d0cd80 --- /dev/null +++ b/Complexity/Algorithms/ExactTriangle/Weights.lean @@ -0,0 +1,61 @@ +/- +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.Defs +public import Mathlib.Algebra.Order.Group.Abs +public import Mathlib.Basic.Real.Basic +import Mathlib.Tactic.Linarith + +/-! +# Integer triangle sums and magnitude bounds + +## Main definitions + +* `weightSum`: the sum of the three supplied integer edge weights. +* `WeightsBoundedBy`: a shared real bound on their absolute values. + +## Main results + +* `abs_weightSum_le`: a triangle's magnitude is at most three times the edge bound. + +## References + +Adapted from the weight sum and bound definitions in [Anthropic's +PaperStatements.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/PaperStatements.lean#L893-L914) +and the triangle inequality in its +[Hashing.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Hashing.lean#L246-L255). +The existing `HasZeroTriangle` remains the decision predicate. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.ExactTriangle + +variable {n : ℕ} (wAB wBC wAC : Fin n → Fin n → ℤ) + +/-- Sum of the A–B, B–C and A–C weights at a triangle. -/ +def weightSum (a b c : Fin n) : ℤ := wAB a b + wBC b c + wAC a c + +/-- Every supplied edge weight has absolute value at most `U`. -/ +def WeightsBoundedBy (U : ℝ) : Prop := + (∀ a b, ((|wAB a b| : ℤ) : ℝ) ≤ U) ∧ + (∀ b c, ((|wBC b c| : ℤ) : ℝ) ≤ U) ∧ ∀ a c, ((|wAC a c| : ℤ) : ℝ) ≤ U + +/-- The three edge bounds imply a threefold bound on every triangle sum. -/ +theorem abs_weightSum_le {U : ℝ} (h : WeightsBoundedBy wAB wBC wAC U) (a b c : Fin n) : + ((|weightSum wAB wBC wAC a b c| : ℤ) : ℝ) ≤ 3 * U := by + obtain ⟨hAB, hBC, hAC⟩ := h + have habs : |weightSum wAB wBC wAC a b c| ≤ |wAB a b| + |wBC b c| + |wAC a c| := + abs_add_three _ _ _ + have hcast : ((|weightSum wAB wBC wAC a b c| : ℤ) : ℝ) ≤ + ((|wAB a b| : ℤ) : ℝ) + ((|wBC b c| : ℤ) : ℝ) + ((|wAC a c| : ℤ) : ℝ) := by + exact_mod_cast habs + linarith [hAB a b, hBC b c, hAC a c] + +end Cslib.Algorithms.ExactTriangle diff --git a/Complexity/Algorithms/WordRAM/Deduplicate/Count.lean b/Complexity/Algorithms/WordRAM/Deduplicate/Count.lean new file mode 100644 index 0000000000..7ecf8a3c0a --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Deduplicate/Count.lean @@ -0,0 +1,148 @@ +/- +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.Deduplicate.Model +public import Complexity.Computability.WordRAM.ListArray + +/-! +# Native counting of one adjacent run + +The counter tests the array bound before every read and stops at the first different word. +No sorting or counting oracle is used. The mathematical end index is a proof parameter and +does not occur in the code or a register. + +## Main definitions + +* `Count.body`, `Count.code`: the native run-counting loop. + +## Main results + +* `Count.exec`: every matching input position advances the cursor, with linear charged time. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Deduplicate + +open Cslib.WordRAM + +/-- Input length. -/ +abbrev len : Fin 16 := 0 +/-- Sorted input array. -/ +abbrev src : Fin 16 := 1 +/-- Unique-key output array. -/ +abbrev keys : Fin 16 := 2 +/-- Multiplicity output array. -/ +abbrev counts : Fin 16 := 3 +/-- Input cursor. -/ +abbrev idx : Fin 16 := 4 +/-- Number of completed output runs. -/ +abbrev out : Fin 16 := 5 +/-- Current run key. -/ +abbrev key : Fin 16 := 6 +/-- Current run multiplicity. -/ +abbrev count : Fin 16 := 7 +/-- Loaded input word. -/ +abbrev value : Fin 16 := 8 +/-- Address temporary. -/ +abbrev addr : Fin 16 := 9 +/-- Constant one. -/ +abbrev one : Fin 16 := 10 +/-- Run-counting loop flag. -/ +abbrev active : Fin 16 := 11 + +namespace Count + +/-- Read the next word only if the cursor is in bounds; count it or end the run. -/ +def body : Code .standard 16 := ram[ + if idx < len then [ + addr ← src + idx; value ← *addr; + if value = key then [ count ← count + one; idx ← idx + one ] + else [ active ← #0 ] + ] else [ active ← #0 ] ] + +/-- Count from the current cursor, starting with zero count and flag one. -/ +def code : Code .standard 16 := .loop (.eq active one) body + +/-- Count a maximal adjacent run in `8c+10` instructions, including the terminating test. +`stop` is a proof parameter: execution learns the boundary solely through charged comparisons +and loads. The source memory is unchanged. Empty runs are also admitted by this kernel. -/ +theorem exec {w s n start stop k : ℕ} (x : Data 16) + (hn : x.regs len = n) (hs : x.regs src = s) (hi : x.regs idx = start) + (hk : x.regs key = k) (hc : x.regs count = 0) + (ha : x.regs active = 1) (ho : x.regs one = 1) + (hle : start ≤ stop) (hstop : stop ≤ n) + (hmem : ∀ j, start ≤ j → j < stop → x.mem (s + j) = k) + (hnext : stop < n → x.mem (s + stop) ≠ k) + (hfit : Fits w (s + n)) : + ∃ y, ExecWithin w code x y (8 * (stop - start) + 10) ∧ + y.regs idx = stop ∧ y.regs count = stop - start ∧ y.mem = x.mem := by + let I (y : Data 16) : Prop := + y.regs len = n ∧ y.regs src = s ∧ y.regs key = k ∧ y.regs one = 1 ∧ + start ≤ y.regs idx ∧ y.regs idx ≤ stop ∧ y.regs count = y.regs idx - start ∧ + y.regs active ≤ 1 ∧ (y.regs active = 0 → y.regs idx = stop) ∧ y.mem = x.mem + have hbody (y : Data 16) (hy : I y) (hcond : (Cond.eq active one).Holds y) : + ∃ z, ExecWithin w body y z 6 ∧ I z ∧ + stop - z.regs idx + z.regs active < stop - y.regs idx + y.regs active := by + obtain ⟨hyn, hys, hyk, hyo, hlo, hhi, hyc, hya, hyend, hym⟩ := hy + have hyflag : y.regs active = 1 := by simpa only [Cond.holds_eq, hyo] using hcond + simp only [len, src, key, one, idx, count, active] at * + have hlenLink : y.regs len = y.regs 0 := rfl + have hidxLink : y.regs idx = y.regs 4 := rfl + have hcountLink : y.regs count = y.regs 7 := rfl + have hactiveLink : y.regs active = y.regs 11 := rfl + by_cases hlt : y.regs idx < stop + · have hidx : y.regs idx < y.regs len := by omega + have hread : y.mem (y.regs src + y.regs idx) = k := by + rw [hym, hys]; exact hmem _ hlo hlt + have had : Fits w (y.regs src + y.regs idx) := by + rw [hys]; exact hfit.of_le (by omega) + have hic : Fits w (y.regs idx + 1) := hfit.of_le (by omega) + have hcc : Fits w (y.regs count + 1) := hfit.of_le (by omega) + refine ⟨?_, ?run, ?_, ?_⟩ + case run => ram_eval [body, hidx, hread, hyk, hyo, had, hic, hcc] + · simp +decide only [I, Data.set_regs, Data.set_mem, Function.update_apply, + idx, count, value, addr, len, src, key, one, active, ↓reduceIte] + exact ⟨hyn, hys, hyk, hyo, by omega, by omega, by omega, by omega, by omega, hym⟩ + · simp +decide only [Data.set_regs, Function.update_apply, idx, count, value, addr, + ↓reduceIte] + omega + · have hend : y.regs idx = stop := by omega + by_cases hsn : stop < n + · have hidx : y.regs idx < y.regs len := by omega + have hread : y.mem (y.regs src + y.regs idx) ≠ y.regs key := by + simpa only [hym, hys, hend, hyk] using hnext hsn + have had : Fits w (y.regs src + y.regs idx) := by + rw [hys]; exact hfit.of_le (by omega) + refine ⟨?_, ?run, ?_, ?_⟩ + case run => ram_eval [body, hidx, hread, had] + · simp +decide only [I, Data.set_regs, Data.set_mem, Function.update_apply, + idx, count, value, addr, len, src, key, one, active, ↓reduceIte] + exact ⟨hyn, hys, hyk, hyo, hlo, hhi, hyc, by omega, fun _ => hend, hym⟩ + · simp +decide [Data.set_regs, hyflag, hend] + · have hidx : ¬y.regs idx < y.regs len := by omega + refine ⟨?_, ?run, ?_, ?_⟩ + case run => ram_eval [body, hidx] + · simp +decide only [I, Data.set_regs, Data.set_mem, Function.update_apply, + idx, count, len, src, key, one, active, ↓reduceIte] + exact ⟨hyn, hys, hyk, hyo, hlo, hhi, hyc, by omega, fun _ => hend, hym⟩ + · simp +decide [Data.set_regs, hyflag, hend] + have hx : I x := ⟨hn, hs, hk, ho, by omega, by omega, by omega, by omega, by omega, rfl⟩ + obtain ⟨y, he, hy, hnot⟩ := ExecWithin.loop_invariant I + (fun y => stop - y.regs idx + y.regs active) hbody hx + obtain ⟨_, _, _, hyo, _, _, hyc, hya, hyend, hym⟩ := hy + have hz : y.regs active = 0 := by + simp only [Cond.holds_eq, hyo] at hnot + omega + have hend := hyend hz + refine ⟨y, ?_, hend, by simpa only [hend] using hyc, hym⟩ + simpa only [code, hi, ha, Nat.add_mul, Nat.one_mul, Nat.mul_comm, Nat.add_assoc] using he + +end Count +end Cslib.Algorithms.WordRAM.Deduplicate diff --git a/Complexity/Algorithms/WordRAM/Deduplicate/Model.lean b/Complexity/Algorithms/WordRAM/Deduplicate/Model.lean new file mode 100644 index 0000000000..fb6feb225b --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Deduplicate/Model.lean @@ -0,0 +1,146 @@ +/- +Copyright (c) 2026 Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Data.List.Sort +public import Mathlib.Data.List.Nodup + +/-! +# Adjacent deduplication with full multiplicities + +The output is a list of key/count pairs with positive counts and strictly increasing keys. +Expanding those runs recovers the sorted input exactly. This specifies the whole finite set and +all multiplicities, rather than a one-way membership condition. These proofs are original. + +## Main definitions + +* `expand`: expand key/count pairs into their repeated keys. +* `RunEncoding`: positive, strictly ordered runs reconstructing the input. + +## Main results + +* `RunEncoding.keys_nodup`: each key occurs once. +* `RunEncoding.mem_iff`: the keys enumerate exactly the input set. +* `RunEncoding.multiplicity`: each count is the full input multiplicity. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Deduplicate + +/-- Expand a list of key/count pairs into repeated keys. -/ +def expand (rs : List (ℕ × ℕ)) : List ℕ := rs.flatMap (fun p => List.replicate p.2 p.1) + +/-- Runs with positive counts and strictly increasing keys reconstruct exactly `xs`. -/ +structure RunEncoding (xs : List ℕ) (rs : List (ℕ × ℕ)) : Prop where + /-- Each emitted run is nonempty. -/ + positive : ∀ p ∈ rs, 0 < p.2 + /-- A key never reappears after a different key. -/ + ordered : rs.Pairwise (fun p q => p.1 < q.1) + /-- All input positions are accounted for, with their exact order and multiplicity. -/ + reconstruct : expand rs = xs + +/-- Expansion distributes over concatenation. -/ +theorem expand_append (rs ss : List (ℕ × ℕ)) : expand (rs ++ ss) = expand rs ++ expand ss := by + simp [expand] + +/-- Expansion of one run. -/ +theorem expand_singleton (k c : ℕ) : expand [(k, c)] = List.replicate c k := by + simp [expand] + +/-- Every key in an expansion occurs in some nonempty run. -/ +theorem mem_expand {rs : List (ℕ × ℕ)} {k : ℕ} : + k ∈ expand rs ↔ ∃ p ∈ rs, p.1 = k ∧ 0 < p.2 := by + simp only [expand, List.mem_flatMap, List.mem_replicate] + constructor + · rintro ⟨p, hp, hc, hk⟩ + exact ⟨p, hp, hk.symm, by omega⟩ + · rintro ⟨p, hp, hk, hc⟩ + exact ⟨p, hp, by omega, hk.symm⟩ + +/-- Positive runs contain at least one element each. -/ +theorem length_le_expand {rs : List (ℕ × ℕ)} (hp : ∀ p ∈ rs, 0 < p.2) : + rs.length ≤ (expand rs).length := by + induction rs with + | nil => simp [expand] + | cons p rs ih => + have hc := hp p (by simp) + have ht := ih (fun q hq => hp q (by simp [hq])) + simp only [expand, List.flatMap_cons, List.length_append, List.length_replicate, + List.length_cons] at * + omega + +private theorem count_expand {rs : List (ℕ × ℕ)} + (ho : rs.Pairwise (fun p q => p.1 < q.1)) {p : ℕ × ℕ} (hp : p ∈ rs) : + (expand rs).count p.1 = p.2 := by + induction rs with + | nil => simp at hp + | cons q rs ih => + obtain ⟨hq, ht⟩ := List.pairwise_cons.mp ho + rcases List.mem_cons.mp hp with rfl | hp + · have hnot : p.1 ∉ expand rs := by + rintro hk + obtain ⟨p, hp, he, _⟩ := mem_expand.mp hk + have := hq p hp + omega + change (List.replicate p.2 p.1 ++ expand rs).count p.1 = p.2 + rw [List.count_append, List.count_replicate_self, List.count_eq_zero_of_not_mem hnot, + Nat.add_zero] + · have hlt := hq p hp + have hc := ih ht hp + simp only [expand, List.flatMap_cons, List.count_append] + rw [show (List.replicate q.2 q.1).count p.1 = 0 by + simp [List.count_replicate, Nat.ne_of_lt hlt]] + simpa only [Nat.zero_add, expand] using hc + +/-- The key array contains no repeated key. -/ +theorem RunEncoding.keys_nodup {xs : List ℕ} {rs : List (ℕ × ℕ)} + (h : RunEncoding xs rs) : (rs.map Prod.fst).Nodup := by + rw [List.nodup_iff_pairwise_ne, List.pairwise_map] + exact h.ordered.imp (fun hlt => Nat.ne_of_lt hlt) + +/-- The key array enumerates precisely the input set. -/ +theorem RunEncoding.mem_iff {xs : List ℕ} {rs : List (ℕ × ℕ)} + (h : RunEncoding xs rs) (k : ℕ) : k ∈ rs.map Prod.fst ↔ k ∈ xs := by + rw [← h.reconstruct, mem_expand] + rw [List.mem_map] + constructor + · rintro ⟨p, hp, hk⟩ + exact ⟨p, hp, hk, h.positive p hp⟩ + · rintro ⟨p, hp, hk, _⟩ + exact ⟨p, hp, hk⟩ + +/-- Each emitted count equals the full multiplicity in the input list. -/ +theorem RunEncoding.multiplicity {xs : List ℕ} {rs : List (ℕ × ℕ)} + (h : RunEncoding xs rs) {p : ℕ × ℕ} (hp : p ∈ rs) : xs.count p.1 = p.2 := by + rw [← h.reconstruct] + exact count_expand h.ordered hp + +/-- A nonempty sorted list starts with a positive run of its head, followed only by larger +keys. This is a mathematical decomposition used to verify the adjacent scan, not an oracle +operation of the program. -/ +theorem first_run (k : ℕ) (xs : List ℕ) (hs : (k :: xs).Pairwise (· ≤ ·)) : + ∃ c rest, 0 < c ∧ k :: xs = List.replicate c k ++ rest ∧ ∀ v ∈ rest, k < v := by + induction xs generalizing k with + | nil => exact ⟨1, [], by decide, rfl, by simp⟩ + | cons j xs ih => + obtain ⟨hk, ht⟩ := List.pairwise_cons.mp hs + by_cases heq : k = j + · subst j + obtain ⟨c, rest, hc, he, hr⟩ := ih k ht + refine ⟨c + 1, rest, by omega, ?_, hr⟩ + simpa only [List.replicate_succ, List.cons_append] using congrArg (List.cons k) he + · have hkj : k < j := lt_of_le_of_ne (hk j (by simp)) heq + refine ⟨1, j :: xs, by decide, rfl, ?_⟩ + intro v hv + rcases List.mem_cons.mp hv with rfl | hv + · exact hkj + · exact lt_of_lt_of_le hkj ((List.pairwise_cons.mp ht).1 v hv) + +end Cslib.Algorithms.WordRAM.Deduplicate diff --git a/Complexity/Algorithms/WordRAM/Deduplicate/Program.lean b/Complexity/Algorithms/WordRAM/Deduplicate/Program.lean new file mode 100644 index 0000000000..df4da01c3d --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Deduplicate/Program.lean @@ -0,0 +1,284 @@ +/- +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.Deduplicate.Count + +/-! +# Complete native adjacent deduplication and multiplicity arrays + +The input is a sorted array. The program counts one maximal adjacent run, emits its key and +full count, and repeats. The output consists of two adjacent arrays, each with capacity equal +to the input length. All run tests, loads, stores, address calculations and control instructions +are charged. The total time is linear, including repeated keys and empty inputs. + +## Main definitions + +* `flush`: store one key/count pair and advance the output cursor. +* `body`: count and emit one complete run. +* `code`: the full adjacent deduplication program. + +## Main results + +* `exec`: full unique-key and multiplicity output with a linear native instruction bound. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Deduplicate + +open Cslib.WordRAM + +/-- Store one key/count pair and advance the output cursor. -/ +def flush : Code .standard 16 := ram[ + addr ← keys + out; *addr ← key; + addr ← counts + out; *addr ← count; out ← out + one ] + +/-- Count and emit the next maximal run. -/ +def body : Code .standard 16 := + .seq (ram[ addr ← src + idx; key ← *addr; count ← #0; active ← #1 ]) + (.seq Count.code flush) + +/-- Initialize the two cursors and count all adjacent runs. -/ +def code : Code .standard 16 := + .seq (ram[ one ← #1; idx ← #0; out ← #0 ]) (.loop (.lt idx len) body) + +private theorem flush_exec {w : ℕ} (x : Data 16) (ho : x.regs one = 1) + (hk : Fits w (x.regs keys + x.regs out)) + (hc : Fits w (x.regs counts + x.regs out)) (hi : Fits w (x.regs out + 1)) : + ∃ y, ExecWithin w flush x y 5 ∧ y.regs out = x.regs out + 1 ∧ + y.mem = Function.update (Function.update x.mem (x.regs keys + x.regs out) (x.regs key)) + (x.regs counts + x.regs out) (x.regs count) := by + refine ⟨?_, ?run, ?_, ?_⟩ + case run => ram_eval [flush, ho, hk, hc, hi] + all_goals simp +decide [Data.set_regs, Data.write_mem] + +/-- Emit each distinct key once and its full multiplicity in at most `29n+5` instructions. +The input precedes the two output blocks; no input cells or cells outside the output blocks +change. Counts are ordinary unsigned words bounded by `n`. No restriction on key signs or +canonical encodings is needed: in particular canonical `Signed.encode` inputs are admitted. -/ +theorem exec {w s base : ℕ} (xs : List ℕ) (x : Data 16) + (hn : x.regs len = xs.length) (hs : x.regs src = s) + (hk : x.regs keys = base) (hc : x.regs counts = base + xs.length) + (hmem : Block.HoldsList s xs x.mem) (hsorted : xs.Pairwise (· ≤ ·)) + (hbelow : s + xs.length ≤ base) (hfit : Fits w (base + 2 * xs.length)) + (h1 : Fits w 1) : + ∃ y rs, ExecWithin w code x y (29 * xs.length + 5) ∧ + RunEncoding xs rs ∧ y.regs out = rs.length ∧ + Block.HoldsList base (rs.map Prod.fst) y.mem ∧ + Block.HoldsList (base + xs.length) (rs.map Prod.snd) y.mem ∧ + Block.HoldsList s xs y.mem ∧ + Modifies (Block.mk base (2 * xs.length)).footprint x.mem y.mem := by + let n := xs.length + let I (y : Data 16) : Prop := + y.regs len = n ∧ y.regs src = s ∧ y.regs keys = base ∧ + y.regs counts = base + n ∧ y.regs one = 1 ∧ y.regs idx ≤ n ∧ + Block.HoldsList s xs y.mem ∧ + Modifies (Block.mk base (2 * n)).footprint x.mem y.mem ∧ + ∃ rs, RunEncoding (xs.take (y.regs idx)) rs ∧ y.regs out = rs.length ∧ + Block.HoldsList base (rs.map Prod.fst) y.mem ∧ + Block.HoldsList (base + n) (rs.map Prod.snd) y.mem ∧ + ∀ p ∈ rs, ∀ v ∈ xs.drop (y.regs idx), p.1 < v + have hbody (y : Data 16) (hy : I y) (hcond : (Cond.lt idx len).Holds y) : + ∃ z T, ExecWithin w body y z T ∧ I z ∧ + T + 2 + 29 * (n - z.regs idx) ≤ 29 * (n - y.regs idx) := by + obtain ⟨hyn, hys, hyk, hyc, hyo, hiy, hym, hyf, rs, hrs, hyout, hyrk, hyrc, hrlt⟩ := hy + have hi : y.regs idx < xs.length := by simpa only [Cond.holds_lt, hyn] using hcond + have htail := hsorted.sublist (List.drop_sublist (y.regs idx) xs) + rw [List.drop_eq_getElem_cons hi] at htail + obtain ⟨c, rest, hcp, heq, hrest⟩ := first_run xs[y.regs idx] (xs.drop (y.regs idx + 1)) htail + have hdrop : xs.drop (y.regs idx) = List.replicate c xs[y.regs idx] ++ rest := by + rw [List.drop_eq_getElem_cons hi, heq] + have hlen := congrArg List.length hdrop + simp only [List.length_drop, List.length_append, List.length_replicate] at hlen + have hstop : y.regs idx + c ≤ n := by omega + have hdrop' : xs.drop (y.regs idx + c) = rest := by + rw [← List.drop_drop, hdrop] + exact List.drop_left' (by simp) + have htake : xs.take (y.regs idx + c) = + xs.take (y.regs idx) ++ List.replicate c xs[y.regs idx] := by + rw [List.take_add, hdrop, List.take_left' (by simp)] + have hrun (j : ℕ) (hlo : y.regs idx ≤ j) (hhi : j < y.regs idx + c) : + y.mem (s + j) = xs[y.regs idx] := by + rw [hym.get j (by omega)] + have hh : (xs.drop (y.regs idx))[j - y.regs idx]'(by simp only [List.length_drop]; omega) = + xs[y.regs idx] := by + simp [hdrop, show j - y.regs idx < c from by omega] + simpa only [List.getElem_drop, show y.regs idx + (j - y.regs idx) = j from by omega] using hh + have hnext (hlt : y.regs idx + c < n) : + y.mem (s + (y.regs idx + c)) ≠ xs[y.regs idx] := by + rw [hym.get _ hlt] + have hmember : xs[y.regs idx + c] ∈ rest := by + rw [← hdrop'] + have hh := List.getElem_mem (l := xs.drop (y.regs idx + c)) + (n := 0) (by simp only [List.length_drop]; omega) + simpa using hh + have := hrest _ hmember + omega + let y₀ := (((y.set addr (s + y.regs idx)).set key xs[y.regs idx]).set count 0).set active 1 + have hinit : ExecWithin (B := .standard) w + (ram[ addr ← src + idx; key ← *addr; count ← #0; active ← #1 ]) y y₀ 4 := by + have had : Fits w (y.regs src + y.regs idx) := by + rw [hys]; exact hfit.of_le (by omega) + have hread : y.mem (y.regs src + y.regs idx) = xs[y.regs idx] := by + rw [hys]; exact hym.get _ hi + have had' : Fits w (s + y.regs idx) := hfit.of_le (by omega) + have hread' : y.mem (s + y.regs idx) = xs[y.regs idx] := hym.get _ hi + ram_eval [y₀, had, had', hread, hread', hys, h1] + obtain ⟨u, hu, hui, huc, hum⟩ := Count.exec (k := xs[y.regs idx]) y₀ + (by simpa +decide [y₀] using hyn) (by simpa +decide [y₀] using hys) + (by simp +decide [y₀]) (by simp +decide [y₀]) + (by simp +decide [y₀]) (by simp +decide [y₀]) + (by simpa +decide [y₀] using hyo) (show y.regs idx ≤ y.regs idx + c by omega) + hstop (by simpa only [y₀, Data.set_mem] using hrun) + (by simpa only [y₀, Data.set_mem] using hnext) (hfit.of_le (by omega)) + have huc' : u.regs count = c := by simpa using huc + have hur (j : Fin 16) (hj : Count.code.Preserves j) : u.regs j = y₀.regs j := hu.regs_eq hj + have hun : u.regs len = n := by + simpa +decide [y₀] + using (hur len (by decide)).trans + (by simpa +decide [y₀] using hyn) + have hus : u.regs src = s := by + simpa +decide [y₀] + using (hur src (by decide)).trans + (by simpa +decide [y₀] using hys) + have huk : u.regs keys = base := by + simpa +decide [y₀] + using (hur keys (by decide)).trans + (by simpa +decide [y₀] using hyk) + have hucb : u.regs counts = base + n := by + simpa +decide [y₀] + using (hur counts (by decide)).trans + (by simpa +decide [y₀] using hyc) + have huo : u.regs one = 1 := by + simpa +decide [y₀] + using (hur one (by decide)).trans + (by simpa +decide [y₀] using hyo) + have huout : u.regs out = rs.length := by + simpa +decide [y₀] + using (hur out (by decide)).trans + (by simpa +decide [y₀] using hyout) + have hukey : u.regs key = xs[y.regs idx] := by simpa +decide [y₀] using hur key (by decide) + have hrl : rs.length ≤ y.regs idx := by + have hh := length_le_expand hrs.positive + rw [hrs.reconstruct, List.length_take_of_le hiy] at hh + exact hh + obtain ⟨z, hz, hzout, hzm⟩ := flush_exec u huo + (by rw [huk, huout]; exact hfit.of_le (by omega)) + (by rw [hucb, huout]; exact hfit.of_le (by omega)) + (by rw [huout]; exact hfit.of_le (by omega)) + have hzmem : z.mem = Function.update + (Function.update y.mem (base + rs.length) xs[y.regs idx]) (base + n + rs.length) c := by + simpa only [huk, hucb, huout, hukey, huc', hum, y₀, Data.set_mem] using hzm + have hzi : z.regs idx = y.regs idx + c := (hz.regs_eq (by decide)).trans hui + have hzf : Modifies (Block.mk base (2 * n)).footprint y.mem z.mem := by + intro a hnot + have ha : ¬(base ≤ a ∧ a < base + 2 * n) := by + simpa only [Block.mem_footprint, Block.Owns] using hnot + rw [hzmem, Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + have hkeyrep : Block.HoldsList base + ((rs ++ [(xs[y.regs idx], c)]).map Prod.fst) z.mem := by + rw [List.map_append] + simp only [List.map_cons, List.map_nil] + rw [hzmem] + have hkr : Block.HoldsList base (rs.map Prod.fst ++ [xs[y.regs idx]]) + (Function.update y.mem (base + rs.length) xs[y.regs idx]) := by + simpa only [List.length_map] using hyrk.append_one xs[y.regs idx] + apply hkr.frame + intro a hlo hhi + rw [Function.update_of_ne (by + simp only [List.length_append, List.length_map, List.length_singleton] at hhi + omega)] + have hcountrep : Block.HoldsList (base + n) + ((rs ++ [(xs[y.regs idx], c)]).map Prod.snd) z.mem := by + rw [List.map_append] + simp only [List.map_cons, List.map_nil] + rw [hzmem] + have hcr : Block.HoldsList (base + n) (rs.map Prod.snd) + (Function.update y.mem (base + rs.length) xs[y.regs idx]) := by + apply hyrc.frame + intro a hlo hhi + rw [Function.update_of_ne (by omega)] + simpa only [List.length_map] using hcr.append_one c + have hrs' : RunEncoding (xs.take (y.regs idx + c)) (rs ++ [(xs[y.regs idx], c)]) := by + refine ⟨?_, ?_, ?_⟩ + · intro p hp + rcases List.mem_append.mp hp with hp | hp + · exact hrs.positive p hp + · have hp' := List.mem_singleton.mp hp + subst p + exact hcp + · apply List.pairwise_append.mpr + refine ⟨hrs.ordered, by simp, ?_⟩ + intro p hp q hq + have hq' : q = (xs[y.regs idx], c) := List.mem_singleton.mp hq + subst q + apply hrlt p hp + rw [List.drop_eq_getElem_cons hi] + exact List.mem_cons_self + · rw [expand_append, expand_singleton, hrs.reconstruct, htake] + refine ⟨z, 8 * c + 19, ?_, ?_, ?_⟩ + · have hall := hinit.seq (hu.seq hz) + exact hall.mono (by omega) + · refine ⟨(hz.regs_eq (by decide)).trans hun, (hz.regs_eq (by decide)).trans hus, + (hz.regs_eq (by decide)).trans huk, (hz.regs_eq (by decide)).trans hucb, + (hz.regs_eq (by decide)).trans huo, by omega, ?_, ?_, + rs ++ [(xs[y.regs idx], c)], by simpa only [hzi] using hrs', ?_, hkeyrep, hcountrep, ?_⟩ + · apply hym.frame + intro a hlo hhi + apply hzf + rw [Block.mem_footprint] + change ¬(base ≤ a ∧ a < base + 2 * n) + omega + · intro a ha + exact (hzf a ha).trans (hyf a ha) + · rw [hzout, huout, List.length_append, List.length_singleton] + · intro p hp v hv + rw [hzi, hdrop'] at hv + rcases List.mem_append.mp hp with hp | hp + · apply hrlt p hp + rw [hdrop] + exact List.mem_append_right _ hv + · have hp' : p = (xs[y.regs idx], c) := List.mem_singleton.mp hp + subst p + exact hrest v hv + · rw [hzi] + omega + let x₀ := ((x.set one 1).set idx 0).set out 0 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, hmem, Modifies.refl _ _, [], ?_, ?_, ?_, ?_, ?_⟩ + · simpa +decide [x₀] using hn + · simpa +decide [x₀] using hs + · simpa +decide [x₀] using hk + · simpa +decide [x₀] using hc + · simp +decide [x₀] + · simp +decide [x₀] + · refine ⟨by simp, by simp, ?_⟩ + simp +decide [x₀, expand] + · simp [x₀] + · exact Block.holdsList_nil _ _ + · exact Block.holdsList_nil _ _ + · simp + obtain ⟨y, he, hy, hcond⟩ := ExecWithin.loop_potential I + (fun y => 29 * (n - y.regs idx)) hbody hx₀ + obtain ⟨hyn, _, _, _, _, hiy, hym, hyf, rs, hrs, hyout, hyrk, hyrc, _⟩ := hy + have hend : y.regs idx = xs.length := by + simp only [Cond.holds_lt, hyn] at hcond + omega + have hinit : ExecWithin (B := .standard) w (ram[ one ← #1; idx ← #0; out ← #0 ]) x x₀ 3 := by + ram_eval [x₀, h1] + refine ⟨y, rs, ?_, by simpa only [hend, List.take_length] using hrs, + hyout, hyrk, hyrc, hym, hyf⟩ + have hall := hinit.seq he + convert hall using 1 + · rfl + · simp +decide only [x₀, Data.set_regs, Function.update_apply, idx, out, one, + ↓reduceIte, Nat.sub_zero, n] + omega + +end Cslib.Algorithms.WordRAM.Deduplicate diff --git a/Complexity/Algorithms/WordRAM/RadixSort.lean b/Complexity/Algorithms/WordRAM/RadixSort.lean new file mode 100644 index 0000000000..a7409b44cd --- /dev/null +++ b/Complexity/Algorithms/WordRAM/RadixSort.lean @@ -0,0 +1,108 @@ +/- +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.RadixSort.Deduplicate +public import Complexity.Computability.WordRAM.Time + +/-! +# Fixed-code native sorting and multiplicity construction from machine input + +The input is the ordinary native list of unsigned words. The code reads its length and the +physical word width, and uses that width as an explicit key-bit bound. The lower-level +`deduplicateCode` accepts a smaller caller-supplied bit bound when known. +The standard-input wrapper's charged bound is `O(n * w)` for positive `n` and `w`. + +The key array starts at `1+2n`, the count array at `1+3n`; the return value is their common +logical length. The two sorting buffers occupy `[1,1+2n)`. The implementation accepts repeated +keys, zero keys and canonical signed integer encodings. Its order is unsigned encoding order. + +## Main definitions + +* `inputCode`, `program`: the complete structured and flat native programs. +* `layoutWidth`: a logarithmic width sufficient for all four array blocks. + +## Main results + +* `input_exec`: charged execution from ordinary native input, with exact key/count output. +* `layout_fits`: logarithmic-width layout fit, including empty input. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.RadixSort + +open Cslib.WordRAM + +/-- Read the native input header and width, then sort and construct both output arrays. -/ +def inputCode : Code .standard 16 := + .seq (ram[ addr ← #0; len ← *addr; src ← #1; dst ← src + len; bits ← width% ]) + deduplicateCode + +/-- The fixed native instruction list, returning the number of unique keys. -/ +def program : Program .standard 16 := inputCode.program Deduplicate.out + +/-- A logarithmic layout-width schedule; the key magnitude bound is independent. -/ +def layoutWidth (n : ℕ) : ℕ := logWidth 3 (n + 1) + +/-- Four array blocks and both initialization constants fit at the layout-width threshold. +This does not by itself bound input key magnitudes. -/ +theorem layout_fits (n : ℕ) : Fits (layoutWidth n) (1 + 4 * n) ∧ Fits (layoutWidth n) 2 := by + have hs : 4 ≤ (n + 2) ^ 2 := Nat.pow_le_pow_left (show 2 ≤ n + 2 by omega) 2 + have hm := Nat.mul_le_mul_right (n + 2) hs + have hp : 4 * (n + 2) ≤ (n + 2) ^ 3 := by simpa only [Nat.pow_succ] using hm + constructor + · apply fits_of_lt_pow + have he : n + 1 + 1 = n + 2 := by omega + rw [he] + omega + · apply fits_of_lt_pow + have he : n + 1 + 1 = n + 2 := by omega + rw [he] + omega + +/-- Execute the complete finite native program on the standard input convention. The bound +includes the final halt and every setup/address/control instruction. Input fit and the explicit +`4n+1` layout fit are independent requirements. Empty and singleton inputs need no exceptions. +The output words enumerate exactly the input set and each stored count is its full multiplicity. -/ +theorem input_exec {w : ℕ} (xs : List ℕ) (hinput : InputFits w xs) + (hfit : Fits w (1 + 4 * xs.length)) (h2 : Fits w 2) : + ∃ (y : Data 16) (rs : List (ℕ × ℕ)) (steps : ℕ), + Exec w inputCode (Data.init xs) y steps ∧ + Executes w program (State.init xs) rs.length (steps + 1) ∧ + steps + 1 ≤ deduplicateTime xs.length w + 6 ∧ + Block.HoldsList (1 + 2 * xs.length) (rs.map Prod.fst) y.mem ∧ + Block.HoldsList (1 + 3 * xs.length) (rs.map Prod.snd) y.mem ∧ + (rs.map Prod.fst).Nodup ∧ + (rs.map Prod.fst).Pairwise (· < ·) ∧ + (∀ k, k ∈ rs.map Prod.fst ↔ k ∈ xs) ∧ + (∀ p ∈ rs, 0 < p.2 ∧ xs.count p.1 = p.2) ∧ + Modifies (Block.mk 1 (4 * xs.length)).footprint (Data.init xs : Data 16).mem y.mem := by + let x := Data.init xs (r := 16) + have h1 := h2.of_le (show 1 ≤ 2 by omega) + have had : Fits w (1 + xs.length) := hfit.of_le (by omega) + let u := ((((x.set addr 0).set len xs.length).set src 1).set dst (1 + xs.length)).set bits w + have hi : ExecWithin w + (ram[ addr ← #0; len ← *addr; src ← #1; dst ← src + len; bits ← width% ]) x u 5 := by + ram_eval [u, x, Data.init_mem_zero, h1, had] + have hn : u.regs len = xs.length := by simp +decide [u] + have hs : u.regs src = 1 := by simp +decide [u] + have hd : u.regs dst = 1 + xs.length := by simp +decide [u] + have hb : u.regs bits = w := by simp [u] + have hm : Block.HoldsList 1 xs u.mem := by + intro i + change (Data.init xs : Data 16).mem (1 + i.val) = xs[i.val] + rw [Nat.add_comm, Data.init_mem_succ] + simp only [List.getD, List.getElem?_eq_getElem i.isLt, Option.getD_some] + obtain ⟨y, rs, he, hyout, hyrk, hyrc, hnodup, hordered, hset, hcount, hframe⟩ := + deduplicate_exec xs u hn hs hd hb hm hinput.2 hinput.2 hfit (fits_self w) h2 + obtain ⟨steps, hsteps, hexec⟩ := hi.seq he + refine ⟨y, rs, steps, hexec, ?_, by omega, hyrk, hyrc, hnodup, hordered, hset, hcount, hframe⟩ + simpa only [program, inputCode, hyout, State.init, x] using hexec.executes Deduplicate.out + +end Cslib.Algorithms.WordRAM.RadixSort diff --git a/Complexity/Algorithms/WordRAM/RadixSort/Deduplicate.lean b/Complexity/Algorithms/WordRAM/RadixSort/Deduplicate.lean new file mode 100644 index 0000000000..49d3ce5bb2 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/RadixSort/Deduplicate.lean @@ -0,0 +1,136 @@ +/- +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.RadixSort.Program +public import Complexity.Algorithms.WordRAM.Deduplicate.Program +public import Complexity.Computability.WordRAM.Signed.Encoding + +/-! +# Native radix sorting followed by unique-key/multiplicity array construction + +This composes the complete native sorter and adjacent-run builder. Two length-`n` sorting +buffers are followed by length-`n` key and count output buffers. All four blocks lie in one +explicit frame. The returned keys are in unsigned encoding order; decoding signed keys does +not turn that order into signed integer order. + +## Main definitions + +* `deduplicateCode`: fixed 16-register sorting and multiplicity program. +* `deduplicateTime`: complete structured instruction bound. + +## Main results + +* `deduplicate_exec`: exact set enumeration and full multiplicities, with no sorting premise. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.RadixSort + +open Cslib.WordRAM + +/-- Original base retained while radix sorting swaps its two pointers. -/ +abbrev origin : Fin 16 := 14 + +/-- Install the two output pointers after the two sorting buffers. -/ +def outputLayout : Code .standard 16 := ram[ + Deduplicate.keys ← len + len; Deduplicate.keys ← Deduplicate.keys + origin; + Deduplicate.counts ← Deduplicate.keys + len ] + +/-- Complete sorting and deduplication; `Deduplicate.out` returns the number of unique words. -/ +def deduplicateCode : Code .standard 16 := + .seq (ram[ origin ← src ]) (.seq code (.seq outputLayout Deduplicate.code)) + +/-- All instructions, including pointer construction and both complete algorithms, excluding +only the final halt. For positive `n` and key-bit bound this is `O(n * keyBits)`. -/ +def deduplicateTime (n keyBits : ℕ) : ℕ := keyBits * (22 * n + 15) + 29 * n + 14 + +/-- Complete fixed-code radix sorting and adjacent deduplication, with full multiplicities. +Keys below `2^keyBits` may repeat arbitrarily. Only the four length-`n` buffers change. The +conditions constrain arithmetic/address fit, not sorting or deduplication behavior. -/ +theorem deduplicate_exec {w base keyBits : ℕ} (xs : List ℕ) (x : Data 16) + (hn : x.regs len = xs.length) (hs : x.regs src = base) + (hd : x.regs dst = base + xs.length) (hb : x.regs bits = keyBits) + (hmem : Block.HoldsList base xs x.mem) + (hkeys : ∀ v ∈ xs, v < 2 ^ keyBits) (hwords : ∀ v ∈ xs, Fits w v) + (hfit : Fits w (base + 4 * xs.length)) (hbits : Fits w keyBits) (h2 : Fits w 2) : + ∃ (y : Data 16) (rs : List (ℕ × ℕ)), + ExecWithin w deduplicateCode x y (deduplicateTime xs.length keyBits) ∧ + y.regs Deduplicate.out = rs.length ∧ + Block.HoldsList (base + 2 * xs.length) (rs.map Prod.fst) y.mem ∧ + Block.HoldsList (base + 3 * xs.length) (rs.map Prod.snd) y.mem ∧ + (rs.map Prod.fst).Nodup ∧ + (rs.map Prod.fst).Pairwise (· < ·) ∧ + (∀ k, k ∈ rs.map Prod.fst ↔ k ∈ xs) ∧ + (∀ p ∈ rs, 0 < p.2 ∧ xs.count p.1 = p.2) ∧ + Modifies (Block.mk base (4 * xs.length)).footprint x.mem y.mem := by + let x₀ := x.set origin base + have hinit : ExecWithin (B := .standard) w (ram[ origin ← src ]) x x₀ 1 := by ram_eval [x₀, hs] + obtain ⟨u, hu, hum, hp, hsorted, hptr, huf⟩ := exec xs x₀ + (by simpa +decide [x₀] using hn) (by simpa +decide [x₀] using hs) + (by simpa +decide [x₀] using hd) (by simpa +decide [x₀] using hb) + hmem hkeys hwords (hfit.of_le (by omega)) hbits h2 + have hun : u.regs len = xs.length := by + simpa +decide [x₀] using (hu.regs_eq (i := len) (by decide)).trans + (show x₀.regs len = xs.length by simpa +decide [x₀] using hn) + have huo : u.regs origin = base := by + simpa only [x₀, Data.set_regs, Function.update_self] using hu.regs_eq (i := origin) (by decide) + have htwice : Fits w (u.regs len + u.regs len) := hfit.of_le (by rw [hun]; omega) + have hkeyaddr : Fits w (u.regs len + u.regs len + base) := hfit.of_le (by rw [hun]; omega) + have hcntaddr : Fits w (u.regs len + u.regs len + base + u.regs len) := + hfit.of_le (by rw [hun]; omega) + let v := ((u.set Deduplicate.keys (u.regs len + u.regs len)).set + Deduplicate.keys (u.regs len + u.regs len + base)).set + Deduplicate.counts (u.regs len + u.regs len + base + u.regs len) + have hl : ExecWithin w outputLayout u v 3 := by + ram_eval [v, outputLayout, huo, htwice, hkeyaddr, hcntaddr] + have hvl : v.regs Deduplicate.len = (passes keyBits xs).length := by + change u.regs len = (passes keyBits xs).length + rw [hun, hp.length_eq] + have hvs : v.regs Deduplicate.src = u.regs src := rfl + have hvk : v.regs Deduplicate.keys = base + 2 * xs.length := by + change u.regs len + u.regs len + base = base + 2 * xs.length + omega + have hvc : v.regs Deduplicate.counts = base + 3 * xs.length := by + change u.regs len + u.regs len + base + u.regs len = base + 3 * xs.length + omega + have hlen : (passes keyBits xs).length = xs.length := hp.length_eq + obtain ⟨y, rs, hy, hrs, hyout, hyrk, hyrc, _, hyf⟩ := Deduplicate.exec (passes keyBits xs) v + hvl hvs hvk (by rw [hvc, hlen]; omega) hum hsorted + (by rw [hlen]; rcases hptr with hptr | hptr <;> omega) + (by rw [hlen]; exact hfit.of_le (by omega)) (h2.of_le (by omega)) + have heq : base + 2 * xs.length + xs.length = base + 3 * xs.length := by omega + refine ⟨y, rs, ?_, hyout, hyrk, by simpa only [hlen, heq] using hyrc, + hrs.keys_nodup, List.pairwise_map.mpr hrs.ordered, ?_, ?_, ?_⟩ + · have hall := hinit.seq (hu.seq (hl.seq hy)) + apply hall.mono + simp only [deduplicateTime, timeBound, hlen] + omega + · intro k + exact (hrs.mem_iff k).trans hp.mem_iff + · intro p hp' + exact ⟨hrs.positive p hp', (hp.count_eq p.1).symm.trans (hrs.multiplicity hp')⟩ + · intro a hnot + have hnotSort : a ∉ (Block.mk base (2 * xs.length)).footprint := by + intro ha + apply hnot + rw [Block.mem_footprint] at ha ⊢ + dsimp only [Block.Owns] at ha ⊢ + omega + have hnotOut : a ∉ (Block.mk (base + 2 * xs.length) + (2 * (passes keyBits xs).length)).footprint := by + intro ha + apply hnot + rw [Block.mem_footprint] at ha ⊢ + dsimp only [Block.Owns] at ha ⊢ + rw [hlen] at ha + omega + exact (hyf a hnotOut).trans (huf a hnotSort) + +end Cslib.Algorithms.WordRAM.RadixSort diff --git a/Complexity/Algorithms/WordRAM/RadixSort/Model.lean b/Complexity/Algorithms/WordRAM/RadixSort/Model.lean new file mode 100644 index 0000000000..c92aeef96c --- /dev/null +++ b/Complexity/Algorithms/WordRAM/RadixSort/Model.lean @@ -0,0 +1,134 @@ +/- +Copyright (c) 2026 Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Data.List.Sort + +/-! +# Stable binary radix sorting: mathematical specification + +The order is the unsigned order on encoded words. In particular this does not claim that +sign-and-magnitude integer encodings are ordered as integers. The construction and proofs here +are original; the executable native program is supplied separately. + +## Main definitions + +* `digit`: one binary digit, represented as zero or one. +* `select`: the stable subsequence with a specified digit. +* `partition`: zeros followed by ones, preserving the order within each part. +* `passes`: least-significant-digit radix passes. + +## Main results + +* `partition_perm`: a pass preserves the multiset. +* `passes_correct`: sufficiently many passes give a sorted permutation. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.RadixSort + +/-- Binary digit `b` of unsigned word `v`. -/ +def digit (b v : ℕ) : ℕ := v / 2 ^ b % 2 + +/-- The stable subsequence of words whose digit `b` equals `t`. -/ +def select (b t : ℕ) (xs : List ℕ) : List ℕ := + xs.filter (fun v => decide (digit b v = t)) + +/-- Stable zero/one partition at digit `b`. -/ +def partition (b : ℕ) (xs : List ℕ) : List ℕ := select b 0 xs ++ select b 1 xs + +/-- Apply the first `bits` binary radix passes, starting with the least significant bit. -/ +def passes : ℕ → List ℕ → List ℕ + | 0, xs => xs + | b + 1, xs => partition b (passes b xs) + +/-- Every binary digit is below two. -/ +theorem digit_lt_two (b v : ℕ) : digit b v < 2 := Nat.mod_lt _ (by decide) + +/-- The native shift/remainder implementation reads precisely `digit`. -/ +theorem digit_eq_shift (b v : ℕ) : digit b v = (v >>> b) % 2 := by + rw [Nat.shiftRight_eq_div_pow, digit] + +/-- A selection never increases the length. -/ +theorem length_select_le (b t : ℕ) (xs : List ℕ) : (select b t xs).length ≤ xs.length := + List.length_filter_le _ _ + +/-- A selected word occurred in the original list and has the requested digit. -/ +theorem mem_select {b t v : ℕ} {xs : List ℕ} : + v ∈ select b t xs ↔ v ∈ xs ∧ digit b v = t := by + simp [select] + +/-- Appending an input word appends it to precisely its selected subsequence. -/ +theorem select_append (b t : ℕ) (xs ys : List ℕ) : + select b t (xs ++ ys) = select b t xs ++ select b t ys := by + simp [select] + +/-- A stable bit partition is a permutation of the input, including repeated words. -/ +theorem partition_perm (b : ℕ) (xs : List ℕ) : (partition b xs).Perm xs := by + have h (v : ℕ) : decide (digit b v = 1) = !decide (digit b v = 0) := by + have := digit_lt_two b v + by_cases hz : digit b v = 0 + · simp [hz] + · have ho : digit b v = 1 := by omega + simp [ho] + simpa only [partition, select, h] using + List.filter_append_perm (fun v => decide (digit b v = 0)) xs + +/-- Each pass keeps exactly the same number of words. -/ +theorem length_partition (b : ℕ) (xs : List ℕ) : (partition b xs).length = xs.length := + (partition_perm b xs).length_eq + +/-- All radix passes preserve the input multiset. -/ +theorem passes_perm (bits : ℕ) (xs : List ℕ) : (passes bits xs).Perm xs := by + induction bits with + | zero => exact .refl _ + | succ b ih => exact (partition_perm b _).trans ih + +/-- A pass extends sortedness modulo `2^b` to sortedness modulo `2^(b+1)`. -/ +theorem partition_pairwise {b : ℕ} {xs : List ℕ} + (hs : xs.Pairwise (fun u v => u % 2 ^ b ≤ v % 2 ^ b)) : + (partition b xs).Pairwise (fun u v => u % 2 ^ (b + 1) ≤ v % 2 ^ (b + 1)) := by + have he (v : ℕ) : v % 2 ^ (b + 1) = v % 2 ^ b + 2 ^ b * digit b v := by + rw [Nat.pow_succ, Nat.mod_mul, digit] + have hp (t : ℕ) : + (select b t xs).Pairwise (fun u v => u % 2 ^ (b + 1) ≤ v % 2 ^ (b + 1)) := by + apply (hs.filter _).imp_of_mem + intro u v hu hv huv + change u ∈ select b t xs at hu + change v ∈ select b t xs at hv + rw [mem_select] at hu hv + rw [he, he, hu.2, hv.2] + omega + apply List.pairwise_append.mpr + refine ⟨hp 0, hp 1, ?_⟩ + intro u hu v hv + rw [mem_select] at hu hv + rw [he, he, hu.2, hv.2] + have hmod := Nat.mod_lt u (Nat.two_pow_pos b) + omega + +/-- The output is sorted on the low `bits` bits. -/ +theorem passes_pairwise (bits : ℕ) (xs : List ℕ) : + (passes bits xs).Pairwise (fun u v => u % 2 ^ bits ≤ v % 2 ^ bits) := by + induction bits with + | zero => simp [passes, List.pairwise_iff_getElem, Nat.mod_one] + | succ b ih => exact partition_pairwise ih + +/-- Keys strictly below `2^bits` are sorted in unsigned order after `bits` passes. -/ +theorem passes_correct (bits : ℕ) (xs : List ℕ) + (hbound : ∀ v ∈ xs, v < 2 ^ bits) : + (passes bits xs).Perm xs ∧ (passes bits xs).Pairwise (· ≤ ·) := by + refine ⟨passes_perm bits xs, (passes_pairwise bits xs).imp_of_mem ?_⟩ + intro u v hu hv huv + have hu' := hbound u ((passes_perm bits xs).mem_iff.mp hu) + have hv' := hbound v ((passes_perm bits xs).mem_iff.mp hv) + simpa only [Nat.mod_eq_of_lt hu', Nat.mod_eq_of_lt hv'] using huv + +end Cslib.Algorithms.WordRAM.RadixSort diff --git a/Complexity/Algorithms/WordRAM/RadixSort/Program.lean b/Complexity/Algorithms/WordRAM/RadixSort/Program.lean new file mode 100644 index 0000000000..a4d8a6f851 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/RadixSort/Program.lean @@ -0,0 +1,204 @@ +/- +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.RadixSort.Scan +public import Complexity.Computability.WordRAM.Compile + +/-! +# Fixed-code native unsigned radix sorter + +Two stable scans per bit partition the source into zeros followed by ones in scratch. Swapping +the pointers avoids a copy-back scan. Code and registers are independent of length and key width. +The source pointer on exit identifies the sorted array; the other block is scratch. + +## Main definitions + +* `pass`: two digit scans and pointer swap. +* `code`: the complete radix sorter with a register-supplied bit bound. +* `timeBound`: explicit instruction bound, excluding halt. + +## Main results + +* `exec`: complete charged sorting, with permutation, sortedness and a two-array frame. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.RadixSort + +open Cslib.WordRAM + +/-- Exchange source and scratch, then advance to the next bit. -/ +def advance : Code .standard 16 := ram[ + swap ← src; src ← dst; dst ← swap; bit ← bit + one ] + +/-- Exactly two stable scans per bit. -/ +def pass : Code .standard 16 := + .seq (ram[ out ← #0; tag ← #0 ]) + (ram[ ⟪Scan.code⟫; tag ← #1; ⟪Scan.code⟫; ⟪advance⟫ ]) + +/-- Complete unsigned radix sorter. Parameters are length, two disjoint array bases, and bits. +The returned source register identifies the sorted block. -/ +def code : Code .standard 16 := + .seq (ram[ one ← #1; two ← #2; bit ← #0 ]) (.loop (.lt bit bits) pass) + +/-- All instructions of radix sorting, excluding a final halt. -/ +def timeBound (n keyBits : ℕ) : ℕ := keyBits * (22 * n + 15) + 5 + +private theorem pass_exec {w b s d : ℕ} (xs : List ℕ) (x : Data 16) + (hn : x.regs len = xs.length) (hs : x.regs src = s) (hd : x.regs dst = d) + (hb : x.regs bit = b) (ho : x.regs one = 1) (h2 : x.regs two = 2) + (hmem : Block.HoldsList s xs x.mem) (hwords : ∀ v ∈ xs, Fits w v) + (hdisj : s + xs.length ≤ d ∨ d + xs.length ≤ s) + (hsfit : Fits w (s + xs.length)) (hdfit : Fits w (d + xs.length)) + (hbit : Fits w (b + 1)) (h1 : Fits w 1) : + ∃ y, ExecWithin w pass x y (22 * xs.length + 13) ∧ + Block.HoldsList d (partition b xs) y.mem ∧ + y.regs src = d ∧ y.regs dst = s ∧ y.regs bit = b + 1 ∧ + Modifies (Block.mk d xs.length).footprint x.mem y.mem := by + let x₀ := (x.set out 0).set tag 0 + have hinit : ExecWithin (B := .standard) w (ram[ out ← #0; tag ← #0 ]) x x₀ 2 := by + ram_eval [x₀] + obtain ⟨y, hy, hym, hysrc, hyout, hyf⟩ := Scan.exec (t := 0) xs [] x₀ + (by simpa +decide [x₀] using hn) (by simpa +decide [x₀] using hs) + (by simpa +decide [x₀] using hd) (by simpa +decide [x₀] using hb) + (by simp [x₀]) (by simpa +decide [x₀] using ho) + (by simpa +decide [x₀] using h2) (by simp +decide [x₀]) + hmem (Block.holdsList_nil _ _) hwords + (by simpa using length_select_le b 0 xs) hdisj hsfit hdfit + have hyr (j : Fin 16) (hj : Scan.code.Preserves j) : y.regs j = x₀.regs j := + hy.regs_eq hj + let y₀ := y.set tag 1 + have hn₀ : y₀.regs len = xs.length := by + simpa +decide [y₀, hyr len (by decide), x₀] using hn + have hs₀ : y₀.regs src = s := by + simpa +decide [y₀, hyr src (by decide), x₀] using hs + have hd₀ : y₀.regs dst = d := by + simpa +decide [y₀, hyr dst (by decide), x₀] using hd + have hb₀ : y₀.regs bit = b := by + simpa +decide [y₀, hyr bit (by decide), x₀] using hb + have ho₀ : y₀.regs one = 1 := by + simpa +decide [y₀, hyr one (by decide), x₀] using ho + have h2₀ : y₀.regs two = 2 := by + simpa +decide [y₀, hyr two (by decide), x₀] using h2 + obtain ⟨z, hz, hzm, _, _, hzf⟩ := Scan.exec (t := 1) xs (select b 0 xs) y₀ + hn₀ hs₀ hd₀ hb₀ (by simp [y₀]) ho₀ h2₀ + (by simpa +decide [y₀] using hyout) hysrc + (by simpa only [y₀, Data.set_mem, List.nil_append] using hym) hwords + (by rw [← partition, length_partition]) hdisj hsfit hdfit + have hzs : z.regs src = s := (hz.regs_eq (by decide)).trans hs₀ + have hzd : z.regs dst = d := (hz.regs_eq (by decide)).trans hd₀ + have hzb : z.regs bit = b := (hz.regs_eq (by decide)).trans hb₀ + have hzo : z.regs one = 1 := (hz.regs_eq (by decide)).trans ho₀ + refine ⟨?u, ?run, ?_, ?_, ?_, ?_, ?_⟩ + case run => + have ha : ExecWithin w advance z ?u 4 := by + ram_eval [advance, hzs, hzd, hzb, hzo, hbit] + have hall := hinit.seq (hy.seq ((ExecWithin.const (d := tag) (x := y) h1).seq (hz.seq ha))) + exact hall.mono (by omega) + · exact hzm + · simp +decide [Data.set_regs] + · simp +decide [Data.set_regs] + · simp +decide [Data.set_regs] + · intro a hnot + exact (hzf a hnot).trans (hyf a hnot) + +/-- Sort all keys below `2^keyBits` in `keyBits * (22n+15)+5` instructions. The two arrays +have length `n`, are adjacent, and start at caller-supplied `base`. Only these `2n` cells may +change. The sorted block is returned in `src`; it is one of the two original blocks. This +includes empty lists and zero bit bounds. The key-bit bound and both one-past-end addresses +must fit; constant two must fit. No ordering contract is assumed. -/ +theorem exec {w base keyBits : ℕ} (xs : List ℕ) (x : Data 16) + (hn : x.regs len = xs.length) (hs : x.regs src = base) + (hd : x.regs dst = base + xs.length) (hb : x.regs bits = keyBits) + (hmem : Block.HoldsList base xs x.mem) + (hkeys : ∀ v ∈ xs, v < 2 ^ keyBits) + (hwords : ∀ v ∈ xs, Fits w v) + (hfit : Fits w (base + 2 * xs.length)) (hbits : Fits w keyBits) (h2 : Fits w 2) : + ∃ y, ExecWithin w code x y (timeBound xs.length keyBits) ∧ + Block.HoldsList (y.regs src) (passes keyBits xs) y.mem ∧ + (passes keyBits xs).Perm xs ∧ (passes keyBits xs).Pairwise (· ≤ ·) ∧ + (y.regs src = base ∨ y.regs src = base + xs.length) ∧ + Modifies (Block.mk base (2 * xs.length)).footprint x.mem y.mem := by + let n := xs.length + let I (y : Data 16) : Prop := + y.regs len = n ∧ y.regs bits = keyBits ∧ y.regs one = 1 ∧ y.regs two = 2 ∧ + y.regs bit ≤ keyBits ∧ + ((y.regs src = base ∧ y.regs dst = base + n) ∨ + (y.regs src = base + n ∧ y.regs dst = base)) ∧ + Block.HoldsList (y.regs src) (passes (y.regs bit) xs) y.mem ∧ + Modifies (Block.mk base (2 * n)).footprint x.mem y.mem + have hbody (y : Data 16) (hy : I y) (hc : (Cond.lt bit bits).Holds y) : + ∃ z, ExecWithin w pass y z (22 * n + 13) ∧ I z ∧ + keyBits - z.regs bit < keyBits - y.regs bit := by + obtain ⟨hyn, hyB, hyo, hy2, hyb, hyp, hym, hyf⟩ := hy + have hlt : y.regs bit < keyBits := by simpa only [Cond.holds_lt, hyB] using hc + have hl : (passes (y.regs bit) xs).length = n := (passes_perm _ _).length_eq + have hdis : y.regs src + n ≤ y.regs dst ∨ y.regs dst + n ≤ y.regs src := by + rcases hyp with ⟨ha, hb⟩ | ⟨ha, hb⟩ <;> omega + have hsf : Fits w (y.regs src + n) := hfit.of_le (by + rcases hyp with ⟨ha, hb⟩ | ⟨ha, hb⟩ <;> omega) + have hdf : Fits w (y.regs dst + n) := hfit.of_le (by + rcases hyp with ⟨ha, hb⟩ | ⟨ha, hb⟩ <;> omega) + have hwds : ∀ v ∈ passes (y.regs bit) xs, Fits w v := fun v hv => + hwords v ((passes_perm _ _).mem_iff.mp hv) + obtain ⟨z, he, hzm, hzs, hzd, hzb, hzf⟩ := pass_exec (passes (y.regs bit) xs) y + (hyn.trans hl.symm) rfl rfl rfl hyo hy2 hym hwds + (by simpa only [hl] using hdis) (by simpa only [hl] using hsf) + (by simpa only [hl] using hdf) (hbits.of_le (by omega)) (h2.of_le (by omega)) + refine ⟨z, by simpa only [hl] using he, ?_, by omega⟩ + refine ⟨(he.regs_eq (by decide)).trans hyn, (he.regs_eq (by decide)).trans hyB, + (he.regs_eq (by decide)).trans hyo, (he.regs_eq (by decide)).trans hy2, + by omega, ?_, ?_, ?_⟩ + · rw [hzs, hzd] + rcases hyp with ⟨ha, hb⟩ | ⟨ha, hb⟩ + · exact Or.inr ⟨hb, ha⟩ + · exact Or.inl ⟨hb, ha⟩ + · rw [hzs, hzb, passes] + exact hzm + · intro a hnot + have hnot' : a ∉ (Block.mk (y.regs dst) (passes (y.regs bit) xs).length).footprint := by + intro ha + apply hnot + rw [Block.mem_footprint] at ha ⊢ + change y.regs dst ≤ a ∧ + a < y.regs dst + (passes (y.regs bit) xs).length at ha + change base ≤ a ∧ a < base + 2 * n + rw [hl] at ha + rcases hyp with ⟨ha', hb'⟩ | ⟨ha', hb'⟩ <;> omega + exact (hzf a hnot').trans (hyf a hnot) + let x₀ := ((x.set one 1).set two 2).set bit 0 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, Modifies.refl _ _⟩ + · simpa +decide [x₀] using hn + · simpa +decide [x₀] using hb + · simp +decide [x₀] + · simp +decide [x₀] + · simp [x₀] + · left; simpa +decide [x₀, n] using And.intro hs hd + · simpa +decide [x₀, passes, hs] using hmem + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => keyBits - y.regs bit) hbody hx₀ + obtain ⟨hyn, hyB, hyo, hy2, hyb, hyp, hym, hyf⟩ := hy + have hend : y.regs bit = keyBits := by + simp only [Cond.holds_lt, hyB] at hc + omega + have hinit : ExecWithin (B := .standard) w (ram[ one ← #1; two ← #2; bit ← #0 ]) x x₀ 3 := by + ram_eval [x₀, h2, h2.of_le (show 1 ≤ 2 by omega)] + refine ⟨y, ?_, by simpa only [hend] using hym, + (passes_correct keyBits xs hkeys).1, (passes_correct keyBits xs hkeys).2, + hyp.imp And.left And.left, hyf⟩ + have hall := hinit.seq he + convert hall using 1 + · rfl + · simp +decide only [timeBound, n, x₀, Data.set_regs, Function.update_apply, + ↓reduceIte, Nat.sub_zero, Nat.add_assoc, Nat.reduceAdd] + omega + +end Cslib.Algorithms.WordRAM.RadixSort diff --git a/Complexity/Algorithms/WordRAM/RadixSort/Scan.lean b/Complexity/Algorithms/WordRAM/RadixSort/Scan.lean new file mode 100644 index 0000000000..8ebb617df9 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/RadixSort/Scan.lean @@ -0,0 +1,254 @@ +/- +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.RadixSort.Model +public import Complexity.Computability.WordRAM.ListArray + +/-! +# Native stable digit-selection scan + +The scan appends a stable digit-selected subsequence to a disjoint destination array. Both +orders of the two disjoint blocks are allowed, so consecutive radix passes can swap pointers. +All loads, address calculations, stores and loop control are native instructions. + +## Main definitions + +* `Scan.body`: read, test, conditionally append, and advance. +* `Scan.code`: reset the input cursor and scan the array. + +## Main results + +* `Scan.exec`: stable selection with a linear instruction bound and memory frame. +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.RadixSort + +open Cslib.WordRAM + +/-- Array length. -/ +abbrev len : Fin 16 := 0 +/-- Source pointer; the sorter returns the final sorted pointer here. -/ +abbrev src : Fin 16 := 1 +/-- Scratch pointer. -/ +abbrev dst : Fin 16 := 2 +/-- Current bit index. -/ +abbrev bit : Fin 16 := 3 +/-- Digit selected by the current scan. -/ +abbrev tag : Fin 16 := 4 +/-- Constant one. -/ +abbrev one : Fin 16 := 5 +/-- Constant two. -/ +abbrev two : Fin 16 := 6 +/-- Input cursor. -/ +abbrev idx : Fin 16 := 7 +/-- Output cursor. -/ +abbrev out : Fin 16 := 8 +/-- Address temporary. -/ +abbrev addr : Fin 16 := 9 +/-- Loaded key. -/ +abbrev value : Fin 16 := 10 +/-- Current key digit. -/ +abbrev keyDigit : Fin 16 := 11 +/-- Number of key bits, supplied by the caller. -/ +abbrev bits : Fin 16 := 12 +/-- Temporary used to swap the two pointers. -/ +abbrev swap : Fin 16 := 13 + +namespace Scan + +/-- One stable-selection iteration, with every address and counter operation charged. -/ +def body : Code .standard 16 := ram[ + addr ← src + idx; value ← *addr; + keyDigit ← value >>> bit; keyDigit ← keyDigit % two; + if keyDigit = tag then [ addr ← dst + out; *addr ← value; out ← out + one ] + else [ addr ← addr ]; + idx ← idx + one ] + +/-- Reset the input cursor, then perform one stable-selection scan. -/ +def code : Code .standard 16 := ram[ idx ← #0; while idx < len do [ ⟪body⟫ ] ] + +private theorem preserves_low (j : Fin 16) (hj : j.val ≤ 6) : body.Preserves j := by + have hne (k : Fin 16) (hk : 7 ≤ k.val) : (some k : Option (Fin 16)) ≠ some j := by + intro h + have hv := congrArg Fin.val (Option.some.inj h) + omega + change (some addr ≠ some j) ∧ (some value ≠ some j) ∧ + (some keyDigit ≠ some j) ∧ (some keyDigit ≠ some j) ∧ + (((some addr ≠ some j) ∧ (none ≠ some j) ∧ (some out ≠ some j)) ∧ + (some addr ≠ some j)) ∧ (some idx ≠ some j) + have hnone : (none : Option (Fin 16)) ≠ some j := by intro h; cases h + exact ⟨hne addr (by decide), hne value (by decide), hne keyDigit (by decide), + hne keyDigit (by decide), + ⟨⟨hne addr (by decide), hnone, hne out (by decide)⟩, + hne addr (by decide)⟩, + hne idx (by decide)⟩ + +private theorem step {w : ℕ} (x : Data 16) + (ho : x.regs one = 1) (ht : x.regs two = 2) + (hs : Fits w (x.regs src + x.regs idx)) + (hd : Fits w (x.regs dst + x.regs out)) + (hi : Fits w (x.regs idx + 1)) + (hu : digit (x.regs bit) (x.mem (x.regs src + x.regs idx)) = x.regs tag → + Fits w (x.regs out + 1)) + (hv : Fits w (x.mem (x.regs src + x.regs idx))) : + ∃ y, ExecWithin w body x y 9 ∧ + y.regs idx = x.regs idx + 1 ∧ + y.regs out = x.regs out + + (if digit (x.regs bit) (x.mem (x.regs src + x.regs idx)) = x.regs tag then 1 else 0) ∧ + y.mem = (if digit (x.regs bit) (x.mem (x.regs src + x.regs idx)) = x.regs tag then + Function.update x.mem (x.regs dst + x.regs out) (x.mem (x.regs src + x.regs idx)) + else x.mem) := by + have hshift := hv.of_le (Nat.shiftRight_le _ (x.regs bit)) + have hmod := hshift.of_le (Nat.mod_le _ 2) + by_cases h : digit (x.regs bit) (x.mem (x.regs src + x.regs idx)) = x.regs tag + · have heq : (x.mem (x.regs src + x.regs idx) >>> x.regs bit) % 2 = x.regs tag := + (digit_eq_shift _ _).symm.trans h + have hinc := hu h + have htagfit : Fits w (x.regs tag) := heq ▸ hmod + refine ⟨?_, ?run, ?_, ?_, ?_⟩ + case run => ram_eval [body, ho, ht, hs, hd, hi, hinc, hshift, hmod, heq, htagfit] + all_goals simp +decide [Data.set_regs, h] + · have hne : (x.mem (x.regs src + x.regs idx) >>> x.regs bit) % 2 ≠ x.regs tag := by + simpa only [← digit_eq_shift] using h + refine ⟨?_, ?run, ?_, ?_, ?_⟩ + case run => ram_eval [body, ho, ht, hs, hd, hi, hshift, hmod, hne] + all_goals simp +decide [Data.set_regs, h] + +/-- Append a stable digit-selected subsequence in `11n+3` instructions. The destination has +capacity `n`; source and destination blocks are disjoint, in either order. The caller's existing +destination prefix is retained. Only destination cells may change. -/ +theorem exec {w s d b t : ℕ} (xs acc : List ℕ) (x : Data 16) + (hn : x.regs len = xs.length) (hs : x.regs src = s) (hd : x.regs dst = d) + (hb : x.regs bit = b) (ht : x.regs tag = t) + (ho : x.regs one = 1) (h2 : x.regs two = 2) (hout : x.regs out = acc.length) + (hsrc : Block.HoldsList s xs x.mem) (hacc : Block.HoldsList d acc x.mem) + (hwords : ∀ v ∈ xs, Fits w v) + (hcap : (acc ++ select b t xs).length ≤ xs.length) + (hdisj : s + xs.length ≤ d ∨ d + xs.length ≤ s) + (hsfit : Fits w (s + xs.length)) (hdfit : Fits w (d + xs.length)) : + ∃ y, ExecWithin w code x y (11 * xs.length + 3) ∧ + Block.HoldsList d (acc ++ select b t xs) y.mem ∧ + Block.HoldsList s xs y.mem ∧ + y.regs out = (acc ++ select b t xs).length ∧ + Modifies (Block.mk d xs.length).footprint x.mem y.mem := by + let selected (i : ℕ) := acc ++ select b t (xs.take i) + have hselected (i : ℕ) : (selected i).length ≤ xs.length := by + have hsub : (select b t (xs.take i)).Sublist (select b t xs) := + (List.take_sublist i xs).filter _ + simp only [selected, List.length_append] at hcap ⊢ + have := hsub.length_le + omega + let I (y : Data 16) : Prop := + (∀ j : Fin 16, j.val ≤ 6 → y.regs j = x.regs j) ∧ + y.regs idx ≤ xs.length ∧ y.regs out = (selected (y.regs idx)).length ∧ + Block.HoldsList s xs y.mem ∧ Block.HoldsList d (selected (y.regs idx)) y.mem ∧ + Modifies (Block.mk d xs.length).footprint x.mem y.mem + have hbody (y : Data 16) (hy : I y) (hc : (Cond.lt idx len).Holds y) : + ∃ z, ExecWithin w body y z 9 ∧ I z ∧ + xs.length - z.regs idx < xs.length - y.regs idx := by + obtain ⟨hr, hle, hu, hm, ha, hf⟩ := hy + have hyl : y.regs len = xs.length := (hr len (by decide)).trans hn + have hys : y.regs src = s := (hr src (by decide)).trans hs + have hyd : y.regs dst = d := (hr dst (by decide)).trans hd + have hyb : y.regs bit = b := (hr bit (by decide)).trans hb + have hyt : y.regs tag = t := (hr tag (by decide)).trans ht + have hyo : y.regs one = 1 := (hr one (by decide)).trans ho + have hy2 : y.regs two = 2 := (hr two (by decide)).trans h2 + have hi : y.regs idx < xs.length := by simpa only [Cond.holds_lt, hyl] using hc + have hread : y.mem (y.regs src + y.regs idx) = xs[y.regs idx] := by + rw [hys]; exact hm.get _ hi + have hv : Fits w (y.mem (y.regs src + y.regs idx)) := by + rw [hread]; exact hwords _ (List.getElem_mem hi) + have hsa : Fits w (y.regs src + y.regs idx) := by + rw [hys]; exact hsfit.of_le (by omega) + have hda : Fits w (y.regs dst + y.regs out) := by + rw [hyd, hu]; exact hdfit.of_le (by have := hselected (y.regs idx); omega) + have hic : Fits w (y.regs idx + 1) := hsfit.of_le (by omega) + have hsuc : selected (y.regs idx + 1) = selected (y.regs idx) ++ + (if digit b xs[y.regs idx] = t then [xs[y.regs idx]] else []) := by + dsimp only [selected] + rw [List.take_succ_eq_append_getElem hi, select_append] + simp only [select, List.filter_cons, List.filter_nil, decide_eq_true_eq] + split_ifs <;> simp only [List.append_assoc, List.append_nil] + by_cases heq : digit b xs[y.regs idx] = t + · have hlt : y.regs out < xs.length := by + have hlim := hselected (y.regs idx + 1) + rw [hsuc, ite_eq_left heq, List.length_append, List.length_singleton] at hlim + omega + obtain ⟨z, he, hzi, hzu, hzm⟩ := step y hyo hy2 hsa hda hic + (fun _ => hdfit.of_le (by omega)) hv + rw [hyb, hread, hyt, ite_eq_left heq] at hzu hzm + refine ⟨z, he, ?_, by omega⟩ + refine ⟨?_, by omega, ?_, ?_, ?_, ?_⟩ + · intro j hj + exact (he.regs_eq (preserves_low j hj)).trans (hr j hj) + · rw [hzi, hsuc, ite_eq_left heq, List.length_append, List.length_singleton] + omega + · apply hm.frame + intro a hlo hhi + rw [hzm, hyd, Function.update_of_ne (by rcases hdisj with h | h <;> omega)] + · rw [hzi, hsuc, ite_eq_left heq, hzm, hyd, hu] + exact ha.append_one _ + · intro a hnot + have hne : a ≠ d + y.regs out := by + intro heqa + apply hnot + apply Block.mem_footprint.mpr + change d ≤ a ∧ a < d + xs.length + omega + rw [hzm, hyd, Function.update_of_ne hne] + exact hf a hnot + · obtain ⟨z, he, hzi, hzu, hzm⟩ := step y hyo hy2 hsa hda hic + (by simpa only [hyb, hread, hyt] using fun h => (heq h).elim) hv + rw [hyb, hread, hyt, ite_eq_right heq] at hzu hzm + have hsel : selected (y.regs idx + 1) = selected (y.regs idx) := by + rw [hsuc, ite_eq_right heq, List.append_nil] + refine ⟨z, he, ?_, by omega⟩ + refine ⟨?_, by omega, ?_, ?_, ?_, ?_⟩ + · intro j hj + exact (he.regs_eq (preserves_low j hj)).trans (hr j hj) + · rw [hzi, hsel, hzu, Nat.add_zero, hu] + · simpa only [hzm] using hm + · simpa only [hzi, hsel, hzm] using ha + · simpa only [hzm] using hf + let x₀ := x.set idx 0 + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_⟩ + · intro j hj + have hne : j ≠ idx := by + intro h + have hv := congrArg Fin.val h + change j.val = 7 at hv + omega + exact Function.update_of_ne hne _ _ + · simp [x₀] + · simpa +decide [x₀, selected, select] using hout + · exact hsrc + · simpa [x₀, selected, select] using hacc + · exact Modifies.refl _ _ + obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => xs.length - y.regs idx) hbody hx₀ + obtain ⟨hr, hle, hu, hm, ha, hf⟩ := hy + have hend : y.regs idx = xs.length := by + have hyl : y.regs len = xs.length := (hr len (by decide)).trans hn + simp only [Cond.holds_lt, hyl] at hc + omega + have hi := (ExecWithin.const (d := idx) (x := x) (fits_zero w)).seq he + refine ⟨y, ?_, ?_, hm, ?_, hf⟩ + · convert hi using 1 + · rfl + · simp only [x₀, Data.set_regs, Function.update_self, Nat.sub_zero] + omega + · simpa only [hend, selected, List.take_length] using ha + · simpa only [hend, selected, List.take_length] using hu + +end Scan +end Cslib.Algorithms.WordRAM.RadixSort diff --git a/Complexity/Computability/WordRAM/ListArray.lean b/Complexity/Computability/WordRAM/ListArray.lean new file mode 100644 index 0000000000..5d90770831 --- /dev/null +++ b/Complexity/Computability/WordRAM/ListArray.lean @@ -0,0 +1,94 @@ +/- +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.Heap + +/-! +# List-valued block assertions + +These are specializations of `Block.Realizes`, not an additional memory model. + +## Main definitions + +* `Block.HoldsList`: a contiguous block realizes a list. + +## Main results + +* `Block.HoldsList.append_one`: a store immediately after a list appends one element. +* `Block.HoldsList.frame`: writes outside the block preserve its list. +-/ + +@[expose] public section + +namespace Cslib.WordRAM.Block + +/-- The block starting at `base` realizes exactly the entries of `xs`. -/ +abbrev HoldsList (base : ℕ) (xs : List ℕ) (mem : ℕ → ℕ) : Prop := + (Block.mk base xs.length).Realizes mem (fun i => xs[i]) + +/-- Read an entry of a represented list. -/ +theorem HoldsList.get {base : ℕ} {xs : List ℕ} {mem : ℕ → ℕ} + (h : HoldsList base xs mem) (i : ℕ) (hi : i < xs.length) : + mem (base + i) = xs[i] := h ⟨i, hi⟩ + +/-- An empty list imposes no memory constraints. -/ +theorem holdsList_nil (base : ℕ) (mem : ℕ → ℕ) : HoldsList base [] mem := by + intro i + exact Fin.elim0 i + +/-- The prefix of a represented list is represented in the same memory. -/ +theorem HoldsList.take {base : ℕ} {xs : List ℕ} {mem : ℕ → ℕ} + (h : HoldsList base xs mem) (k : ℕ) : HoldsList base (xs.take k) mem := by + intro i + have hi : i.val < xs.length := lt_of_lt_of_le i.isLt (List.take_sublist k xs).length_le + change mem (base + i.val) = (xs.take k)[i.val] + rw [List.getElem_take] + exact h.get i.val hi + +/-- A store immediately following a represented list appends an element. -/ +theorem HoldsList.append_one {base : ℕ} {xs : List ℕ} {mem : ℕ → ℕ} + (h : HoldsList base xs mem) (v : ℕ) : + HoldsList base (xs ++ [v]) (Function.update mem (base + xs.length) v) := by + intro i + change Function.update mem (base + xs.length) v (base + i.val) = (xs ++ [v])[i.val] + by_cases hi : i.val < xs.length + · rw [Function.update_of_ne (by omega), List.getElem_append_left hi] + exact h.get i.val hi + · have hil := i.isLt + simp only [List.length_append, List.length_singleton] at hil + have he : i.val = xs.length := by omega + simp [he] + +/-- Writes outside the represented block preserve the list. -/ +theorem HoldsList.frame {base : ℕ} {xs : List ℕ} {mem mem' : ℕ → ℕ} + (h : HoldsList base xs mem) + (heq : ∀ a, base ≤ a → a < base + xs.length → mem' a = mem a) : + HoldsList base xs mem' := by + intro i + have hi : i.val < xs.length := i.isLt + change mem' (base + i.val) = xs[i.val] + rw [heq _ (by omega) (by omega)] + exact h.get i.val i.isLt + +/-- Append two adjacent represented blocks. -/ +theorem HoldsList.append {base : ℕ} {xs ys : List ℕ} {mem : ℕ → ℕ} + (h : HoldsList base xs mem) (h' : HoldsList (base + xs.length) ys mem) : + HoldsList base (xs ++ ys) mem := by + intro i + change mem (base + i.val) = (xs ++ ys)[i.val] + by_cases hi : i.val < xs.length + · rw [List.getElem_append_left hi] + exact h.get i.val hi + · rw [List.getElem_append_right (by omega)] + have hil := i.isLt + simp only [List.length_append] at hil + have he : base + xs.length + (i.val - xs.length) = base + i.val := by omega + simpa only [he] using h'.get (i.val - xs.length) (by omega) + +end Cslib.WordRAM.Block diff --git a/Complexity/Foundations/NumberTheory/PrimeWindow.lean b/Complexity/Foundations/NumberTheory/PrimeWindow.lean new file mode 100644 index 0000000000..9693a53028 --- /dev/null +++ b/Complexity/Foundations/NumberTheory/PrimeWindow.lean @@ -0,0 +1,344 @@ +/- +Copyright (c) 2020 Patrick Stevens, 2026 Anthropic, PBC. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Patrick Stevens, Bolton Bailey, Anthropic, PBC +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Analysis.SpecialFunctions.Log.Base +public import Mathlib.Analysis.SpecialFunctions.Pow.Real +public import Mathlib.NumberTheory.PrimeCounting +import Mathlib.Analysis.Complex.ExponentialBounds +import Mathlib.Data.Nat.Prime.Factorial +import Mathlib.NumberTheory.Bertrand +import Mathlib.NumberTheory.Primorial +import Mathlib.Analysis.Convex.SpecificFunctions.Deriv + +/-! +# Explicit prime counts in a window of ratio two + +For every real `x ≥ 4`, the half-open window `[x/2, x)` contains at least `c * x / log x` primes +for one universal positive constant `c`. The proof counts prime factors of a central binomial +coefficient, so no prime-number-theorem hypothesis is assumed. + +## Main results + +* `Nat.exists_prime_half_le_and_lt`: a prime exists in every such window. +* `Nat.exists_forall_mul_div_log_le_card_primes`: the uniform counting bound. +* `Nat.exists_forall_mul_sqrt_div_log_le_card_primes`: the bound at `x = sqrt D`, for all `D ≥ + 16`. + +## References + +Adapted from [Anthropic's +PrimesInWindow.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Util/PrimesInWindow.lean). +Its central-binomial and concavity estimates adapt `Mathlib/NumberTheory/Bertrand.lean`, +originally by Patrick Stevens and Bolton Bailey; that source credit is retained. +-/ + +@[expose] public section + +open Finset + +namespace Nat + +private theorem one_le_log_of_three_le {x : ℝ} (hx : 3 ≤ x) : 1 ≤ Real.log x := + (Real.le_log_iff_exp_le (by linarith)).2 (Real.exp_one_lt_three.le.trans hx) + +/-- The primes `p` with `n < p ≤ 2n`: the primes that Bertrand's postulate is about. -/ +private def bertrandPrimes (n : ℕ) : Finset ℕ := (Ioc n (2 * n)).filter Nat.Prime + +/-- Membership in `bertrandPrimes`, spelled out. -/ +private theorem mem_bertrandPrimes {n p : ℕ} : + p ∈ bertrandPrimes n ↔ (n < p ∧ p ≤ 2 * n) ∧ p.Prime := by + simp only [bertrandPrimes, mem_filter, mem_Ioc] + +/-! +## 1. Erdős's upper bound for the central binomial coefficient -/ + +/-- Every prime power in `C(2n, n)` is at most `2n`, so `|s|` of them multiply to at most +`(2n)^|s|`. +-/ +private theorem prod_pow_factorization_le {n : ℕ} (hn : 0 < n) (s : Finset ℕ) : + ∏ p ∈ s, p ^ n.centralBinom.factorization p ≤ (2 * n) ^ s.card := + prod_le_pow_card _ _ _ fun _ _ => Nat.pow_factorization_choose_le (by omega) + +/-- The primes in `(√(2n), n]` contribute at most `4^(2n/3)` to `C(2n, n)`: their exponent is at +most 1, it is 0 above `2n/3`, and the product of the primes up to `2n/3` is at most `4^(2n/3)`. +Adapted from Mathlib, Mathlib/NumberTheory/Bertrand.lean, Apache 2.0. +-/ +private theorem prod_middle_le {n : ℕ} (hn : 2 < n) (s : Finset ℕ) + (hs : ∀ p ∈ s, p.Prime ∧ Nat.sqrt (2 * n) < p ∧ p ≤ n) : + ∏ p ∈ s, p ^ n.centralBinom.factorization p ≤ 4 ^ (2 * n / 3) := + calc ∏ p ∈ s, p ^ n.centralBinom.factorization p + ≤ ∏ p ∈ s, (if p ≤ 2 * n / 3 then p else 1) := by + refine prod_le_prod fun p hp => ?_ + obtain ⟨hprime, hsqrt, hpn⟩ := hs p hp + split_ifs with h3 + · exact (pow_right_mono₀ hprime.one_lt.le + (Nat.factorization_choose_le_one (Nat.sqrt_lt'.mp hsqrt))).trans (pow_one p).le + · rw [Nat.factorization_centralBinom_of_two_mul_self_lt_three_mul hn hpn (by omega), + pow_zero] + _ = ∏ p ∈ s.filter (· ≤ 2 * n / 3), p := (prod_filter _ _).symm + _ ≤ primorial (2 * n / 3) := by + refine prod_le_prod_of_subset_of_one_le (fun p hp => ?_) + fun p hp _ => (mem_filter.1 hp).2.one_lt.le + obtain ⟨hps, hp3⟩ := mem_filter.1 hp + exact mem_filter.2 ⟨mem_range.2 (by omega), (hs p hps).1⟩ + _ ≤ 4 ^ (2 * n / 3) := primorial_le_four_pow _ + +/-- Erdős's bound, keeping the primes of `(n, 2n]`: `C(2n, n) ≤ (2n)^√(2n) · 4^(2n/3) · (2n)^k`, +where `k` is the number of primes in `(n, 2n]`. Mathlib proves the case `k = 0` on the way to +Bertrand's postulate (`centralBinom_le_of_no_bertrand_prime`). Adapted from Mathlib, +Mathlib/NumberTheory/Bertrand.lean, Apache 2.0. +-/ +private theorem centralBinom_le (n : ℕ) (hn : 2 < n) : + n.centralBinom ≤ + (2 * n) ^ Nat.sqrt (2 * n) * 4 ^ (2 * n / 3) * (2 * n) ^ (bertrandPrimes n).card := by + -- `C(2n, n)` is the product of its prime powers over the set `T` of the primes up to `2n`. + set T := (range (2 * n + 1)).filter Nat.Prime with hT + have hprod : ∏ p ∈ T, p ^ n.centralBinom.factorization p = n.centralBinom := + (prod_filter_of_ne fun p _ h => by + contrapose h + rw [Nat.factorization_eq_zero_of_not_prime _ h, pow_zero]).trans + n.prod_pow_factorization_centralBinom + -- Split the primes at `n`, and those up to `n` at `√(2n)`. + rw [← hprod, ← prod_filter_mul_prod_filter_not T (· ≤ n), + ← prod_filter_mul_prod_filter_not (T.filter (· ≤ n)) (· ≤ Nat.sqrt (2 * n))] + refine mul_le_mul' (mul_le_mul' ?_ ?_) ?_ + · -- There are at most `√(2n)` primes up to `√(2n)`. + have hsub : (T.filter (· ≤ n)).filter (· ≤ Nat.sqrt (2 * n)) ⊆ Icc 1 (Nat.sqrt (2 * n)) := by + intro p hp + simp only [hT, mem_filter, mem_range] at hp + obtain ⟨⟨⟨-, hprime⟩, -⟩, hsqrt⟩ := hp + exact mem_Icc.mpr ⟨hprime.one_lt.le, hsqrt⟩ + exact (prod_pow_factorization_le (by omega) _).trans (pow_right_mono₀ (by omega) + ((card_le_card hsub).trans_eq (by rw [Nat.card_Icc, Nat.add_sub_cancel]))) + · -- The primes in `(√(2n), n]`. + refine prod_middle_le hn _ fun p hp => ?_ + simp only [hT, mem_filter, mem_range, not_le] at hp + obtain ⟨⟨⟨-, hprime⟩, hpn⟩, hsqrt⟩ := hp + exact ⟨hprime, hsqrt, hpn⟩ + · -- The primes above `n` are those of `(n, 2n]`. + refine (prod_pow_factorization_le (by omega) _).trans + (pow_right_mono₀ (by omega) (card_le_card fun p hp => ?_)) + simp only [hT, mem_filter, mem_range] at hp + obtain ⟨⟨hp2n, hprime⟩, hnp⟩ := hp + exact mem_bertrandPrimes.2 ⟨⟨by omega, by omega⟩, hprime⟩ + +/-! +## 2. The elementary estimate `n^6 · (2n)^(6√(2n)) ≤ 2^n` for `n ≥ 16562` -/ + +/-- For `a ≥ 0` the function `log x + √(2x) · log (2x) - a x` is concave for `x > 1/2`. Adapted from +Mathlib, Mathlib/NumberTheory/Bertrand.lean, Apache 2.0 (there inside +`Bertrand.real_main_inequality`, for `a = log 4 / 3`). +-/ +private theorem concaveOn_log_add_sqrt_mul_log_sub (a : ℝ) (ha : 0 ≤ a) : + ConcaveOn ℝ (Set.Ioi (1 / 2 : ℝ)) + (fun x : ℝ => Real.log x + √(2 * x) * Real.log (2 * x) - a * x) := by + apply ConcaveOn.sub + · apply ConcaveOn.add + · exact strictConcaveOn_log_Ioi.concaveOn.subset + (Set.Ioi_subset_Ioi (by norm_num)) (convex_Ioi _) + · have h := strictConcaveOn_sqrt_mul_log_Ioi.concaveOn.comp_linearMap + ((2 : ℝ) • (LinearMap.id : ℝ →ₗ[ℝ] ℝ)) + have hset : (⇑((2 : ℝ) • (LinearMap.id : ℝ →ₗ[ℝ] ℝ))) ⁻¹' Set.Ioi 1 = Set.Ioi (1 / 2) := by + ext x + simp only [Set.mem_Ioi, Set.mem_preimage, LinearMap.smul_apply, LinearMap.id_coe, id_eq, + smul_eq_mul] + constructor <;> intro h <;> linarith + rw [hset] at h + exact h + · exact ConvexOn.smul ha (convexOn_id (convex_Ioi _)) + +/-- `log x + √(2x) · log (2x) ≤ (log 2 / 6) · x` for `x ≥ 16562`. The difference of the two sides +is concave, it is nonnegative at `x = 2` and nonpositive at `x = 16562 = 2 · 91²` (where `√(2x) += 182`; the only numerical input is `91^13 ≤ 2^85`), so it is nonpositive from there on. Adapted +from Mathlib, Mathlib/NumberTheory/Bertrand.lean, Apache 2.0 (the same argument proves +`Bertrand.real_main_inequality` there, with other numbers). +-/ +private theorem log_add_sqrt_mul_log_le {x : ℝ} (hx : 16562 ≤ x) : + Real.log x + √(2 * x) * Real.log (2 * x) ≤ Real.log 2 / 6 * x := by + have hlog2 : 0 < Real.log 2 := Real.log_pos (by norm_num) + have hat2 : (0 : ℝ) ≤ + Real.log 2 + √(2 * 2) * Real.log (2 * 2) - Real.log 2 / 6 * 2 := by + rw [Real.sqrt_mul_self (by norm_num), Real.log_mul (by norm_num) (by norm_num)] + linarith + have hat16562 : Real.log 16562 + √(2 * 16562) * Real.log (2 * 16562) + - Real.log 2 / 6 * 16562 ≤ 0 := by + have hlog91 : 13 * Real.log 91 ≤ 85 * Real.log 2 := by + have h := Real.log_le_log (by positivity) (show (91 : ℝ) ^ 13 ≤ 2 ^ 85 by norm_num) + rw [Real.log_pow, Real.log_pow] at h + exact_mod_cast h + have hlogx : Real.log 16562 = Real.log 2 + 2 * Real.log 91 := by + rw [show (16562 : ℝ) = 2 * 91 ^ 2 by norm_num, Real.log_mul (by norm_num) (by norm_num), + Real.log_pow] + norm_num + have hlog2x : Real.log (182 * 182) = 2 * Real.log 2 + 2 * Real.log 91 := by + rw [show (182 * 182 : ℝ) = 2 ^ 2 * 91 ^ 2 by norm_num, + Real.log_mul (by norm_num) (by norm_num), Real.log_pow, Real.log_pow] + norm_num + rw [show (2 * 16562 : ℝ) = 182 * 182 by norm_num, Real.sqrt_mul_self (by norm_num), hlogx, + hlog2x] + -- `366 log 91 ≤ 366 · (85/13) log 2 < (16562/6 - 365) log 2` + linarith only [hlog91, hlog2] + have hconcave := concaveOn_log_add_sqrt_mul_log_sub (Real.log 2 / 6) (by positivity) + have hatx := hconcave.right_le_of_le_left'' + (x := 2) (y := 16562) (z := x) (by norm_num [Set.mem_Ioi]) (by rw [Set.mem_Ioi]; linarith) + (by norm_num) hx (hat16562.trans hat2) + linarith only [hatx, hat16562] + +/-- The same estimate for natural numbers, without logarithms. -/ +private theorem pow_mul_pow_sqrt_le_two_pow {n : ℕ} (hn : 16562 ≤ n) : + n ^ 6 * (2 * n) ^ (6 * Nat.sqrt (2 * n)) ≤ 2 ^ n := by + have hx : (16562 : ℝ) ≤ n := by exact_mod_cast hn + have hreal := log_add_sqrt_mul_log_le hx + have hsqrt : (Nat.sqrt (2 * n) : ℝ) ≤ √(2 * n) := by + exact_mod_cast Real.nat_sqrt_le_real_sqrt (a := 2 * n) + have hlog : 0 ≤ Real.log (2 * n) := Real.log_nonneg (by linarith) + have hpos : (0 : ℝ) < n := by linarith + rw [← Nat.cast_le (α := ℝ)] + push_cast + rw [← Real.log_le_log_iff (by positivity) (by positivity), + Real.log_mul (by positivity) (by positivity), Real.log_pow, Real.log_pow, Real.log_pow] + push_cast + -- `6 log n + 6 ⌊√(2n)⌋ log(2n) ≤ 6 (log n + √(2n) log(2n)) ≤ n log 2` + linarith [hreal, mul_le_mul_of_nonneg_right hsqrt hlog] + +/-! +## 3. The count for `n ≥ 16562` -/ + +/-- `2^n ≤ (2n)^(2k)` for `n ≥ 16562`. If it failed, the sixth power of `4^n < n · (2n)^√(2n) · +4^(2n/3) · (2n)^k` would give `2^(12n) < (n^6 (2n)^(6√(2n))) · 2^(8n) · 2^(3n) ≤ 2^(12n)`. +-/ +private theorem two_pow_le_of_large {n : ℕ} (hn : 16562 ≤ n) : + 2 ^ n ≤ (2 * n) ^ (2 * (bertrandPrimes n).card) := by + by_contra hcon + have hfew : (2 * n) ^ (2 * (bertrandPrimes n).card) ≤ 2 ^ n := (not_le.mp hcon).le + have herdos : 4 ^ n < n * ((2 * n) ^ Nat.sqrt (2 * n) * 4 ^ (2 * n / 3) + * (2 * n) ^ (bertrandPrimes n).card) := + (Nat.four_pow_lt_mul_centralBinom n (by omega)).trans_le + (Nat.mul_le_mul_left _ (centralBinom_le n (by omega))) + have hsmall := pow_mul_pow_sqrt_le_two_pow hn + have hmiddle : 4 ^ (6 * (2 * n / 3)) ≤ 4 ^ (4 * n) := + Nat.pow_le_pow_right (by norm_num) (by omega) + refine lt_irrefl ((4 ^ n) ^ 6) ?_ + calc (4 ^ n) ^ 6 + < (n * ((2 * n) ^ Nat.sqrt (2 * n) * 4 ^ (2 * n / 3) + * (2 * n) ^ (bertrandPrimes n).card)) ^ 6 := Nat.pow_lt_pow_left herdos (by norm_num) + _ = (n ^ 6 * (2 * n) ^ (6 * Nat.sqrt (2 * n))) * 4 ^ (6 * (2 * n / 3)) + * ((2 * n) ^ (2 * (bertrandPrimes n).card)) ^ 3 := by ring + _ ≤ 2 ^ n * 4 ^ (4 * n) * (2 ^ n) ^ 3 := by gcongr + _ = (4 ^ n) ^ 6 := by + rw [mul_comm 4 n, pow_mul 4 n 4, show (4 : ℕ) ^ n = 2 ^ n * 2 ^ n by rw [← mul_pow]; rfl] + ring + +/-- The count of step 3 with natural logarithms: `n log 2 ≤ 2 k log (2n)`. -/ +private theorem mul_log_two_le_card_bertrandPrimes_mul_log {n : ℕ} (hn : 16562 ≤ n) : + (n : ℝ) * Real.log 2 ≤ 2 * (bertrandPrimes n).card * Real.log (2 * n) := by + have hpow : (2 : ℝ) ^ n ≤ (2 * (n : ℝ)) ^ (2 * (bertrandPrimes n).card) := by + exact_mod_cast two_pow_le_of_large hn + have hlog := Real.log_le_log (by positivity) hpow + rw [Real.log_pow, Real.log_pow] at hlog + exact_mod_cast hlog + +/-! +## 4. Half-open windows `[x/2, x)` with real endpoints -/ + +/-- For a real `x ≥ 4` the natural number `n = ⌈x/2⌉ - 1`, written below as `(⌈x⌉ - 1) / 2`, +satisfies `n ≥ 1`, `n ≥ x/4` and `2n < x`, and every natural number in `(n, 2n]` lies in `[x/2, +x)`. +-/ +private theorem exists_nat_window {x : ℝ} (hx : 4 ≤ x) : + ∃ n : ℕ, 1 ≤ n ∧ x / 4 ≤ (n : ℝ) ∧ 2 * (n : ℝ) < x ∧ + ∀ p : ℕ, n < p → p ≤ 2 * n → x / 2 ≤ (p : ℝ) ∧ (p : ℝ) < x := by + have hle : x ≤ (⌈x⌉₊ : ℝ) := Nat.le_ceil x + have hlt : (⌈x⌉₊ : ℝ) < x + 1 := Nat.ceil_lt_add_one (by linarith) + have hceil : 4 ≤ ⌈x⌉₊ := by exact_mod_cast hx.trans hle + obtain ⟨n, hn⟩ : ∃ n, n = (⌈x⌉₊ - 1) / 2 := ⟨_, rfl⟩ + have hlow : (⌈x⌉₊ : ℝ) ≤ 2 * n + 2 := by exact_mod_cast (by omega : ⌈x⌉₊ ≤ 2 * n + 2) + have hhigh : (2 * n + 1 : ℝ) ≤ (⌈x⌉₊ : ℝ) := by exact_mod_cast (by omega : 2 * n + 1 ≤ ⌈x⌉₊) + refine ⟨n, by omega, by linarith, by linarith, fun p hnp hp2n => ?_⟩ + have hnp' : (n : ℝ) + 1 ≤ p := by exact_mod_cast hnp + have hp2n' : (p : ℝ) ≤ 2 * n := by exact_mod_cast hp2n + exact ⟨by linarith, by linarith⟩ + +/-- For every real `x ≥ 4` there is a prime `p` with `x/2 ≤ p < x`, by Bertrand's postulate. -/ +theorem exists_prime_half_le_and_lt (x : ℝ) (hx : 4 ≤ x) : + ∃ p : ℕ, p.Prime ∧ x / 2 ≤ (p : ℝ) ∧ (p : ℝ) < x := by + obtain ⟨n, hn1, -, -, hwindow⟩ := exists_nat_window hx + obtain ⟨p, hprime, hnp, hp2n⟩ := Nat.exists_prime_lt_and_le_two_mul n (by omega) + exact ⟨p, hprime, hwindow p hnp hp2n⟩ + +/-- If `(n, 2n]` lies in the window and `n ≥ x/4` is large, the window has at least `x / (12 log x)` +primes: `(x/4) log 2 ≤ n log 2 ≤ 2 k log (2n) ≤ 2 k log x`, and `log 2 ≥ 2/3`. +-/ +private theorem div_log_le_card_of_large {x : ℝ} {n : ℕ} {S : Finset ℕ} (hx : 4 ≤ x) + (hn : 16562 ≤ n) (hn4 : x / 4 ≤ (n : ℝ)) (h2n : 2 * (n : ℝ) < x) (hsub : bertrandPrimes n ⊆ S) : + x / (12 * Real.log x) ≤ (S.card : ℝ) := by + have hcard : ((bertrandPrimes n).card : ℝ) ≤ S.card := by exact_mod_cast card_le_card hsub + have hn1 : (1 : ℝ) ≤ n := by exact_mod_cast (by omega : 1 ≤ n) + have hlog : Real.log (2 * n) ≤ Real.log x := Real.log_le_log (by linarith) h2n.le + have hlog0 : 0 ≤ Real.log (2 * n) := Real.log_nonneg (by linarith) + have hlog2 : 2 / 3 ≤ Real.log 2 := by linarith [Real.log_two_gt_d9] + rw [div_le_iff₀ (mul_pos (by norm_num) (Real.log_pos (by linarith)))] + calc x = 6 * (x / 4 * (2 / 3)) := by ring + _ ≤ 6 * (n * Real.log 2) := by gcongr + _ ≤ 6 * (2 * (bertrandPrimes n).card * Real.log (2 * n)) := + mul_le_mul_of_nonneg_left (mul_log_two_le_card_bertrandPrimes_mul_log hn) (by norm_num) + _ ≤ 6 * (2 * S.card * Real.log x) := by gcongr + _ = S.card * (12 * Real.log x) := by ring + +/-- There is a `c > 0` such that, for every real `x ≥ 4`, every finite set of natural numbers that +contains all primes `p` with `x/2 ≤ p < x` has at least `c · x / log x` elements. (This form +applies to the set of those primes however it is written down.) +-/ +theorem exists_forall_mul_div_log_le_card_primes : + ∃ c : ℝ, 0 < c ∧ ∀ x : ℝ, 4 ≤ x → ∀ S : Finset ℕ, + (∀ p : ℕ, p.Prime → x / 2 ≤ (p : ℝ) → (p : ℝ) < x → p ∈ S) → + c * (x / Real.log x) ≤ (S.card : ℝ) := by + -- `66248 = 4 · 16562`: for `x` below it one prime is enough + refine ⟨1 / 66248, by norm_num, fun x hx S hS => ?_⟩ + obtain ⟨n, hn1, hn4, h2n, hwindow⟩ := exists_nat_window hx + have hlogx : 1 ≤ Real.log x := one_le_log_of_three_le (by linarith) + rcases le_or_gt 16562 n with hlarge | hsmall + · -- The primes of `(n, 2n]` lie in the window. + have hsub : bertrandPrimes n ⊆ S := fun p hp => by + obtain ⟨⟨hnp, hp2n⟩, hprime⟩ := mem_bertrandPrimes.1 hp + exact hS p hprime (hwindow p hnp hp2n).1 (hwindow p hnp hp2n).2 + calc 1 / 66248 * (x / Real.log x) ≤ 1 / 12 * (x / Real.log x) := by gcongr; norm_num + _ = x / (12 * Real.log x) := by field_simp + _ ≤ (S.card : ℝ) := div_log_le_card_of_large hx hlarge hn4 h2n hsub + · -- Here `x ≤ 4n < 66248`, and one prime is enough. + obtain ⟨p, hprime, hlow, hhigh⟩ := exists_prime_half_le_and_lt x hx + have hone : (1 : ℝ) ≤ S.card := by + exact_mod_cast Finset.card_pos.mpr ⟨p, hS p hprime hlow hhigh⟩ + have hn : (n : ℝ) ≤ 16562 := by exact_mod_cast hsmall.le + calc 1 / 66248 * (x / Real.log x) ≤ 1 / 66248 * (x / 1) := by gcongr + _ ≤ 1 := by linarith + _ ≤ (S.card : ℝ) := hone + +/-! +## 5. The window `[√D/2, √D)` used in the modular reduction -/ + +/-- There is a `c > 0` such that, for every natural number `D ≥ 16`, every finite set of natural +numbers that contains all primes `p` with `√D/2 ≤ p < √D` has at least `c · √D / log D` +elements. +-/ +theorem exists_forall_mul_sqrt_div_log_le_card_primes : + ∃ c : ℝ, 0 < c ∧ ∀ D : ℕ, 16 ≤ D → ∀ S : Finset ℕ, + (∀ p : ℕ, p.Prime → √D / 2 ≤ (p : ℝ) → (p : ℝ) < √D → p ∈ S) → + c * (√D / Real.log D) ≤ (S.card : ℝ) := by + obtain ⟨c, hc, hcount⟩ := exists_forall_mul_div_log_le_card_primes + refine ⟨c, hc, fun D hD S hS => ?_⟩ + have hlog : 0 < Real.log D := Real.log_pos (by exact_mod_cast (by omega : 1 < D)) + calc c * (√D / Real.log D) ≤ c * (√D / (Real.log D / 2)) := by + gcongr + linarith + _ = c * (√D / Real.log (√D)) := by rw [Real.log_sqrt (Nat.cast_nonneg D)] + _ ≤ (S.card : ℝ) := hcount _ ((Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD))) S hS + +end Nat diff --git a/Complexity/Foundations/RingTheory/CyclicIntRing.lean b/Complexity/Foundations/RingTheory/CyclicIntRing.lean new file mode 100644 index 0000000000..d15f95b193 --- /dev/null +++ b/Complexity/Foundations/RingTheory/CyclicIntRing.lean @@ -0,0 +1,83 @@ +/- +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.RingTheory.AdjoinRoot + +/-! +# Integer polynomials modulo `X^p - 1` + +## Main definitions + +* `CyclicIntRing`: the integer polynomial quotient. +* `CyclicIntRing.coeff`: coefficients in the unique representative of degree less than `p`. + +## Main results + +* `CyclicIntRing.coeff_x_pow`: powers have a single coefficient at the residue of their + exponent. + +## References + +Adapted from [Anthropic's +Hashing.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Hashing.lean#L57-L101). +-/ + +@[expose] public section + +namespace Cslib + +open Polynomial in +/-- The quotient ring `ℤ[X]/(X^p - 1)`. -/ +abbrev CyclicIntRing (p : ℕ) : Type := AdjoinRoot ((X : ℤ[X]) ^ p - C 1) + +namespace CyclicIntRing + +open Polynomial + +variable {p : ℕ} + +/-- The element `x` of `ℤ[x]/(x^p − 1)`. -/ +noncomputable def x (p : ℕ) : CyclicIntRing p := AdjoinRoot.root _ + +/-- The coefficient of degree `r` in an element `z` of the quotient, as a linear map. It is the +coefficient of `x^r` in the unique polynomial of degree less than `p` that represents `z` +(Mathlib's `AdjoinRoot.modByMonicHom` gives that polynomial; it asks for the fact that `x^p − 1` +is monic, which holds as `p ≠ 0`). +-/ +noncomputable def coeff (hp : p ≠ 0) (r : ℕ) : CyclicIntRing p →ₗ[ℤ] ℤ := + lcoeff ℤ r ∘ₗ AdjoinRoot.modByMonicHom (monic_X_pow_sub_C (1 : ℤ) hp) + +/-- The coefficient, in terms of the representing polynomial. -/ +theorem coeff_apply (hp : p ≠ 0) (r : ℕ) (z : CyclicIntRing p) : + coeff hp r z = (AdjoinRoot.modByMonicHom (monic_X_pow_sub_C (1 : ℤ) hp) z).coeff r := + rfl + +/-- In `ℤ[x]/(x^p − 1)` we have `x^p = 1`. -/ +private theorem x_pow_self (p : ℕ) : x p ^ p = 1 := by + have h : AdjoinRoot.mk ((X : ℤ[X]) ^ p - C 1) ((X : ℤ[X]) ^ p - C 1) = 0 := AdjoinRoot.mk_self + rw [map_sub, map_pow, AdjoinRoot.mk_X, C_1, map_one] at h + exact sub_eq_zero.1 h + +/-- The coefficient of `x^r` in the power `x^k` of `ℤ[x]/(x^p − 1)` is 1 if `r = k mod p` and 0 +otherwise. +-/ +theorem coeff_x_pow (hp : p ≠ 0) (k r : ℕ) : + coeff hp r (x p ^ k) = if r = k % p then 1 else 0 := by + have hred : x p ^ k = x p ^ (k % p) := by + conv_lhs => rw [← Nat.div_add_mod k p, pow_add, pow_mul, x_pow_self, one_pow, one_mul] + have hdeg : ((X : ℤ[X]) ^ (k % p)).degree < ((X : ℤ[X]) ^ p - C 1).degree := by + rw [degree_X_pow, degree_X_pow_sub_C (Nat.pos_of_ne_zero hp)] + exact_mod_cast Nat.mod_lt k (Nat.pos_of_ne_zero hp) + rw [coeff_apply, hred, x, ← AdjoinRoot.mk_X, ← map_pow, AdjoinRoot.modByMonicHom_mk, + (modByMonic_eq_self_iff (monic_X_pow_sub_C (1 : ℤ) hp)).2 hdeg, coeff_X_pow] + +end CyclicIntRing + +end Cslib diff --git a/ComplexityTests/ExactTriangleHashing.lean b/ComplexityTests/ExactTriangleHashing.lean new file mode 100644 index 0000000000..41df09a874 --- /dev/null +++ b/ComplexityTests/ExactTriangleHashing.lean @@ -0,0 +1,115 @@ +/- +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.Hashing + +/-! +# Modular Exact Triangle counting and prime-selection regressions + +Concrete zeros, false positives, coefficient wraparound, minimum selection, the integer window +boundary and an admitted instance of the universal bound. These assert no native running time. +-/ + +namespace ComplexityTests.ExactTriangleHashing + +open Cslib Cslib.Algorithms.ExactTriangle + +private def ab (a _ : Fin 2) : ℤ := if a = 0 then 0 else 2 +private def z (_ _ : Fin 2) : ℤ := 0 + +private theorem window16 : primesInRange 16 = {2, 3} := by + rw [primesInRange_eq_filter] + decide + +-- Both endpoints matter: prime 2 is admitted at the lower endpoint for D=16; +-- prime 5 is excluded at the upper endpoint for D=25, and prime 3 is admitted there. +example : 2 ∈ primesInRange 16 ∧ 5 ∉ primesInRange 25 ∧ 3 ∈ primesInRange 25 := by + simp only [mem_primesInRange_iff_sq] + decide + +example : 2 ∉ primesInRange 17 := by rw [mem_primesInRange_iff_sq]; decide + +example : countZeroMod ab z z 2 = 8 ∧ countZeroMod ab z z 3 = 4 := by + unfold countZeroMod weightSum ab z + decide + +example : zeroCount ab z z = 4 ∧ falsePositiveCount ab z z 2 = 4 ∧ + falsePositiveCount ab z z 3 = 0 := by + unfold zeroCount falsePositiveCount IsFalsePositive weightSum ab z + decide + +example : countZeroMod ab z z 0 = 4 ∧ falsePositiveCount ab z z 0 = 0 := by + unfold countZeroMod falsePositiveCount IsFalsePositive weightSum ab z + decide + +-- Counting distinguishes the two candidate primes; the minimum is 3. +example : IsSelectedPrime ab z z 16 3 ∧ ¬IsSelectedPrime ab z z 16 2 := by + simp only [IsSelectedPrime, window16] + unfold countZeroMod weightSum ab z + decide + +-- A tie permits either minimum prime; the predicate does not impose a tie-break. +example : IsSelectedPrime (fun _ _ : Fin 1 => 0) (fun _ _ => 0) (fun _ _ => 0) 16 2 ∧ + IsSelectedPrime (fun _ _ : Fin 1 => 0) (fun _ _ => 0) (fun _ _ => 0) 16 3 := by + constructor <;> refine ⟨by rw [window16]; decide, ?_⟩ <;> intro q _ <;> + simp [countZeroMod, weightSum] + +example : HasZeroTriangle ab z z := ⟨0, 0, 0, rfl⟩ +example : 0 < zeroCount ab z z := (zeroCount_pos_iff ab z z).mpr ⟨0, 0, 0, rfl⟩ + +private def one (_ _ : Fin 2) : ℤ := 1 + +example : ¬HasZeroTriangle one one one := by + rintro ⟨a, b, c, h⟩ + norm_num [one] at h + +example : zeroCount one one one = 0 ∧ countZeroMod one one one 3 = 8 ∧ + falsePositiveCount one one one 3 = 8 := by + unfold zeroCount countZeroMod falsePositiveCount IsFalsePositive weightSum one + decide + +-- A negative residue and a positive residue wrap to exponent zero modulo 3. +example : CyclicIntRing.coeff (by decide : 3 ≠ 0) 0 + ((matP (fun _ _ : Fin 2 => (-1 : ℤ)) 3 * + matQ (fun _ _ : Fin 2 => (1 : ℤ)) 3) 0 0) = 2 := by + rw [coeff_matP_mul_matQ (fun _ _ : Fin 2 => (1 : ℤ)) + (fun _ _ : Fin 2 => (-1 : ℤ)) (by decide) 0 0 0 (by decide)] + decide + +-- These matrices are asymmetric: transposing either input changes the selected count. +private def ac (a c : Fin 2) : ℤ := if a = 0 then if c = 0 then -1 else 2 else 0 +private def bc (b c : Fin 2) : ℤ := if b = 0 ∧ c = 1 then 1 else 0 + +example : CyclicIntRing.coeff (by decide : 3 ≠ 0) 0 ((matP ac 3 * matQ bc 3) 0 0) = 1 := by + rw [coeff_matP_mul_matQ bc ac (by decide) 0 0 0 (by decide)] + unfold ac bc + decide + +example : CyclicIntRing.coeff (by decide : 2 ≠ 0) 1 (CyclicIntRing.x 2 ^ 5) = 1 := by + rw [CyclicIntRing.coeff_x_pow] + norm_num + +example : CyclicIntRing.coeff (by decide : 2 ≠ 0) 3 (CyclicIntRing.x 2 ^ 5) = 0 := by + rw [CyclicIntRing.coeff_x_pow] + norm_num + +example : zeroCount (fun _ _ : Fin 0 => 0) (fun _ _ => 0) (fun _ _ => 0) = 0 := by + simp [zeroCount] + +example : ¬WeightsBoundedBy (fun _ _ : Fin 1 => (0 : ℤ)) (fun _ _ => 0) (fun _ _ => 0) (-1) := by + norm_num [WeightsBoundedBy] + +-- The full bound's n,D,kappa hypotheses and input bound have a concrete witness. +example : (falsePositiveCount (fun _ _ : Fin 16 => 0) (fun _ _ => 0) (fun _ _ => 0) 2 : ℝ) ≤ + Hashing.falsePositiveConst * (1 * (16 : ℝ) ^ 3 * Real.log 16 / Real.sqrt 16) := by + apply Hashing.falsePositiveCount_le_const_mul (by decide : 16 ≤ 16) (by decide : 16 ≤ 16) + (by norm_num : (1 : ℝ) ≤ 1) + · norm_num [WeightsBoundedBy] + · refine ⟨by rw [window16]; decide, ?_⟩ + intro q _ + simp [countZeroMod, weightSum] + +end ComplexityTests.ExactTriangleHashing diff --git a/ComplexityTests/WordRAMRadixSort.lean b/ComplexityTests/WordRAMRadixSort.lean new file mode 100644 index 0000000000..28c8c96cd4 --- /dev/null +++ b/ComplexityTests/WordRAMRadixSort.lean @@ -0,0 +1,67 @@ +import Complexity.Algorithms.WordRAM.RadixSort + +open Cslib.WordRAM +open Cslib.Algorithms.WordRAM + +private def observe (w fuel : ℕ) (xs : List ℕ) : Option (List ℕ × List ℕ × ℕ) := + (Code.eval w fuel RadixSort.inputCode (Data.init xs)).map fun (y, _) => + let m := y.regs Deduplicate.out + ((List.range m).map (fun i => y.mem (1 + 2 * xs.length + i)), + (List.range m).map (fun i => y.mem (1 + 3 * xs.length + i)), m) + +-- Empty input still executes setup and exits both array scans safely. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 3 4 [] = some ([], [], 0)) + +-- The largest fitting four-bit canonical signed encoding occurs at the key boundary. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 4 5 [Signed.encode (-7)] = some ([15], [1], 1)) + +-- Zero is canonical; a repeated zero is one key with full multiplicity two. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 4 5 [0, 0] = some ([0], [2], 1)) + +-- Counts are full multiplicities, including values above three. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 5 6 (List.replicate 4 (Signed.encode (-1))) = some ([3], [4], 1)) + +-- Positive and negative input integers, repeated keys, and zero in unsigned encoding order. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 5 6 (([2, -1, 2, 0, -1] : List ℤ).map Signed.encode) = + some ([0, 3, 4], [1, 2, 2], 3)) + +-- Encoding order deliberately differs from signed integer order. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 4 5 (([-2, 1, -1] : List ℤ).map Signed.encode) = + some ([2, 3, 5], [1, 1, 1], 3)) + +-- Two physical bits cannot address this singleton's count output: layout fit matters. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 2 4 [2] = none) + +-- Key overflow faults even when the singleton's four-buffer layout fits. +/-- info: true -/ +#guard_msgs in +#eval decide (observe 4 5 [16] = none) + +-- The assembled program is executed using the ordinary native initial state. +/-- info: true -/ +#guard_msgs in +#eval decide ((run 4 RadixSort.program 500 (State.init [4, 3, 3])).map Prod.fst = some 2) + +-- The array layout threshold includes zero and small input lengths. +example : RadixSort.layoutWidth 0 = 3 ∧ RadixSort.layoutWidth 2 = 6 := by decide +kernel + +-- A bit bound of zero is meaningful exactly when all unsigned keys are zero. +example : RadixSort.passes 0 [0, 0] = [0, 0] ∧ + ¬([2, 0] : List ℕ).Pairwise (· ≤ ·) := by decide +kernel + +-- One finite instruction list is used at every input length and word width. +example : RadixSort.program.length = 81 := by decide +kernel diff --git a/NOTICE b/NOTICE index 3ab01a2df3..95b9c8d57a 100644 --- a/NOTICE +++ b/NOTICE @@ -1,10 +1,11 @@ The 3SUM/Exact Triangle reductions, their finite counting lemmas and native matrix builder in -Complexity/Algorithms/ExactTriangle/{Defs,Lopsided}.lean, +Complexity/Algorithms/ExactTriangle/{Defs,Weights,Lopsided,Hashing}.lean, Complexity/Algorithms/ThreeSum/, Complexity/Algorithms/WordRAM/ThreeSum/, Complexity/Foundations/Data/Finset/Average.lean, Complexity/Foundations/Data/Matrix/IntegerBits.lean and -Complexity/Foundations/NumberTheory/PrimeCounting.lean are adapted from +Complexity/Foundations/NumberTheory/{PrimeCounting,PrimeWindow}.lean and +Complexity/Foundations/RingTheory/CyclicIntRing.lean are adapted from anthropics/formal-math, 3sum-apsp, revision 16a5ce41304fcdaea4e0be7c132e1c361dc26da5. @@ -13,3 +14,7 @@ Sparse Lopsided Graphs" Copyright 2026 Anthropic, PBC This product is licensed under the Apache License, Version 2.0 (see LICENSE). + +The prime-window proof also adapts Mathlib/NumberTheory/Bertrand.lean, +Copyright (c) 2020 Patrick Stevens; authors Patrick Stevens, Bolton Bailey, +under the Apache License, Version 2.0. diff --git a/ROADMAP.md b/ROADMAP.md index 9dbb4ebb35..a57d6f3c52 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -422,8 +422,9 @@ each statement against the numbered result. Organize modules by subject, not cha has the Convolution-3SUM/Exact Triangle correctness reduction, a charged native single-instance 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; it still needs - modular Exact Triangle hashing/false-positive splitting, tensor/sparse-query programs, native + 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 selection with an input-fit guarantee and explicit overhead; expected-time rules for unbounded @@ -567,4 +568,5 @@ not planned, not as proved. [#1394]: https://github.com/cslib-community/complexitylib/issues/1394 [#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 [#1449]: https://github.com/cslib-community/complexitylib/issues/1449 diff --git a/UPSTREAM.md b/UPSTREAM.md index 87d9cb5288..30a7a8947f 100644 --- a/UPSTREAM.md +++ b/UPSTREAM.md @@ -158,6 +158,13 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he Given its arguments in registers, it writes one instance in at most `131*t² + 5` instructions at width `w + 2`, excluding the halt; its theorem identifies the three encoded matrices and preserves the input. Layout, signed index/filler bounds and target address bounds are explicit. + 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 + [#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. + `Algorithms/ThreeSum/ThreeSets.lean` and `ArrayPacking.lean` provide the generic distinct-position frontend and three-array packing. `Algorithms/ThreeSum/ChanHe/` now adapts the core mathematical reduction: minimum-prime searches, collision/heavy-element bounds, light-element arrays, bounded @@ -174,8 +181,16 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he `Set.EqOn`; no additional weighted-triangle predicate or source machine is introduced. This is the mathematical interface milestone [#1438](https://github.com/cslib-community/complexitylib/issues/1438), with no sparsity assumption, - preprocessing/time claim or solver. Modular hashing/false-positive splitting and tensor/sparse - query programs remain in #1394, alongside the native and asymptotic composition below. + preprocessing/time claim or solver. `Algorithms/ExactTriangle/{Weights,Hashing}.lean` now adapts + the modular counting/coefficient identities, minimum-count prime and uniform false-positive + bound from `Sec3/Theorem17/Hashing.lean`, using supplied matrices and the existing decision + predicate. Its integer-square window characterisation supports later native enumeration. + The quotient/coefficient API lives in `Foundations/RingTheory/CyclicIntRing.lean`; the uniform + 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 diff --git a/references.bib b/references.bib index 56438b2a80..10c958dee3 100644 --- a/references.bib +++ b/references.bib @@ -125,6 +125,18 @@ @article{AllenderRubinstein1988 doi = {10.1137/0217075} } +@misc{AlmanWilliams2026, + author = {Alman, Josh and Vassilevska Williams, Virginia}, + title = {Truly Subquadratic {3SUM} and Truly Subcubic {APSP} via Triangles in + Sparse Lopsided Graphs}, + year = {2026}, + eprint = {2610.06783}, + archivePrefix = {arXiv}, + primaryClass = {cs.DS}, + doi = {10.48550/arXiv.2610.06783}, + url = {https://arxiv.org/abs/2610.06783v1} +} + @article{AlonBoppana1987, author = {Alon, Noga and Boppana, Ravi B.}, title = {The Monotone Circuit Complexity of {B}oolean Functions},