diff --git a/Complexity.lean b/Complexity.lean index c2a3c1647c..867430ea54 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -31,9 +31,18 @@ public import Complexity.Algorithms.Lean.VanEmdeBoas.Insert public import Complexity.Algorithms.Lean.VanEmdeBoas.Layout public import Complexity.Algorithms.Lean.VanEmdeBoas.Successor public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Alphabets +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Batch public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Full public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Identity +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Layout +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Levels +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Magnitude +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Orders +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Parameters public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Pruned +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Sparse +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.ThinProduct +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Tiling public import Complexity.Algorithms.Reachability.BreadthFirst public import Complexity.Algorithms.Reachability.InductiveCounting public import Complexity.Algorithms.Reachability.InductiveCounting.Bounds @@ -1930,6 +1939,8 @@ public import Complexity.Crypto.Protocols.SigmaProtocol.Lyubashevsky.NativeExtra public import Complexity.Crypto.Protocols.SigmaProtocol.SchnorrIdentification public import Complexity.Crypto.Protocols.SigmaProtocol.SchnorrSignature public import Complexity.Crypto.Protocols.Signature.RandomOracle +public import Complexity.Foundations.Algebra.FiniteMagnitude +public import Complexity.Foundations.Algebra.GeometricTail public import Complexity.Foundations.Analysis.BernoulliInequality public import Complexity.Foundations.Analysis.HarmonicMean public import Complexity.Foundations.Analysis.InnerProductSpace.Basic @@ -2034,6 +2045,7 @@ public import Complexity.Foundations.Combinatorics.Expander.Walk public import Complexity.Foundations.Combinatorics.Expander.Zigzag public import Complexity.Foundations.Combinatorics.FiniteCover public import Complexity.Foundations.Combinatorics.FiniteCover.Basic +public import Complexity.Foundations.Combinatorics.FunctionPlaces public import Complexity.Foundations.Combinatorics.FunctionalGraph public import Complexity.Foundations.Combinatorics.FunctionalGraph.Active public import Complexity.Foundations.Combinatorics.FunctionalGraph.CheckpointForest @@ -2145,6 +2157,7 @@ public import Complexity.Foundations.Data.Finset.Average public import Complexity.Foundations.Data.Finset.Expect public import Complexity.Foundations.Data.Finset.Finpartition public import Complexity.Foundations.Data.Finset.Markov +public import Complexity.Foundations.Data.Finset.OrderGlue public import Complexity.Foundations.Data.Finset.Quantile public import Complexity.Foundations.Data.Finset.Quotient public import Complexity.Foundations.Data.Finset.Sampling @@ -2152,6 +2165,7 @@ public import Complexity.Foundations.Data.Finset.SubsetSum public import Complexity.Foundations.Data.Finset.SymmDiff public import Complexity.Foundations.Data.Hamming public import Complexity.Foundations.Data.Hamming.Expansion +public import Complexity.Foundations.Data.Int.LogBounds public import Complexity.Foundations.Data.Int.Order public import Complexity.Foundations.Data.Int.Residue public import Complexity.Foundations.Data.List.Append @@ -2284,6 +2298,7 @@ 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.NumberTheory.SqrtPrimeWindow public import Complexity.Foundations.Order.BooleanSubalgebra public import Complexity.Foundations.Order.Changes public import Complexity.Foundations.Order.UpperLower diff --git a/Complexity/Algorithms/ExactTriangle/Hashing.lean b/Complexity/Algorithms/ExactTriangle/Hashing.lean index 17ee536865..e8d4c7064b 100644 --- a/Complexity/Algorithms/ExactTriangle/Hashing.lean +++ b/Complexity/Algorithms/ExactTriangle/Hashing.lean @@ -9,15 +9,14 @@ module public import Complexity.Algorithms.ExactTriangle.Weights public import Complexity.Foundations.Data.Int.Residue -public import Complexity.Foundations.NumberTheory.PrimeWindow +public import Complexity.Foundations.NumberTheory.SqrtPrimeWindow 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 +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`. @@ -33,7 +32,7 @@ by `C * κ * n^3 * log n / sqrt D`. The prime count used here is proved for ever * `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. +* `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. @@ -53,71 +52,19 @@ The source pin is 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 +/-- A prime in `Nat.sqrtPrimeWindow 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, + p ∈ Nat.sqrtPrimeWindow D ∧ ∀ q ∈ Nat.sqrtPrimeWindow 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) - -/-! -### 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`. -/ @@ -131,7 +78,7 @@ private theorem eq_add_toNat_emod_iff {p r : ℕ} (hp : p ≠ 0) (hr : r < p) (u ← Int.add_emod] /-! -### Counting the triples with `S(a,b,c) ≡ 0 (mod p)` -/ +### Counting triples with a vanishing modular weight sum -/ section Counting @@ -174,7 +121,8 @@ theorem countZeroMod_eq (p : ℕ) : (countZeroMod wAB wBC wAC) p = (falsePositiv 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. + -- A triple whose `weightSum` is divisible by `p` 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] @@ -231,8 +179,8 @@ theorem falsePositiveCount_add_zeroCount_eq_sum_coeff {p : ℕ} (hp : p ≠ 0) : 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. + -- A zero `weightSum` modulo `p` means `wAC a c + wBC b c ≡ -wAB a b`, and the right side + -- 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 @@ -246,24 +194,6 @@ theorem falsePositiveCount_add_zeroCount_eq_sum_coeff {p : ℕ} (hp : p ≠ 0) : 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 -/ @@ -275,18 +205,19 @@ variable {n D p : ℕ} {κ : ℝ} (wAB wBC wAC : Fin n → Fin n → ℤ) minimum of the modular-zero count. -/ theorem exists_isSelectedPrime (D : ℕ) (hD : 16 ≤ D) : ∃ p, IsSelectedPrime wAB wBC wAC D p := by + have hsqrt : (4 : ℝ) ≤ Real.sqrt D := Real.four_le_sqrt_natCast hD 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⟩ + Nat.exists_prime_half_le_and_lt (Real.sqrt D) hsqrt + exact Finset.exists_min_image (Nat.sqrtPrimeWindow D) (countZeroMod wAB wBC wAC) ⟨q, + Nat.mem_sqrtPrimeWindow.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 + ∀ q ∈ Nat.sqrtPrimeWindow 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 @@ -301,17 +232,17 @@ 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 : ℝ) + (((Nat.sqrtPrimeWindow 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 hsqrt := Real.four_le_sqrt_natCast hD + set divisors := (Nat.sqrtPrimeWindow 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 + obtain ⟨hprime, hge, -⟩ := Nat.mem_sqrtPrimeWindow.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^κ`. + -- Their product divides the nonzero `weightSum`, whose absolute value is at most `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 @@ -335,18 +266,18 @@ edge weights bounded by `n^κ`, the sum of false-positive counts is at most `n^3 -/ theorem sum_falsePositiveCount_le (hD : 16 ≤ D) (hκ : 1 ≤ κ) (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) : - ((∑ p ∈ primesInRange D, (falsePositiveCount wAB wBC wAC) p : ℕ) : ℝ) + ((∑ p ∈ Nat.sqrtPrimeWindow D, (falsePositiveCount wAB wBC wAC) p : ℕ) : ℝ) ≤ (n : ℝ) ^ 3 * Real.logb (Real.sqrt D / 2) (3 * (n : ℝ) ^ κ) := by classical - have hsqrt := four_le_sqrt hD + have hsqrt := Real.four_le_sqrt_natCast 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 + have hswap : ∑ p ∈ Nat.sqrtPrimeWindow D, (falsePositiveCount wAB wBC wAC) p = + ∑ t : Fin n × Fin n × Fin n, + ((Nat.sqrtPrimeWindow 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 : ℝ) + (((Nat.sqrtPrimeWindow 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 @@ -359,7 +290,7 @@ theorem sum_falsePositiveCount_le (hD : 16 ≤ D) (hκ : 1 ≤ κ) (hT : Weights 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 : ℝ) + (((Nat.sqrtPrimeWindow 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 @@ -370,34 +301,24 @@ theorem sum_falsePositiveCount_le (hD : 16 ≤ D) (hκ : 1 ≤ κ) (hT : Weights 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) p * (Nat.sqrtPrimeWindow D).card ≤ ∑ q ∈ Nat.sqrtPrimeWindow 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) : +private 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] + have hsqrt : 0 < Real.sqrt D := Real.sqrt_pos.mpr + (by exact_mod_cast (show 0 < D by omega)) -- `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) @@ -410,9 +331,9 @@ theorem log_le_four_mul_log_sqrt_div_two (D : ℕ) (hD : 16 ≤ D) : 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 : ℝ} +theorem falsePositiveCount_le_of_le_card_sqrtPrimeWindow {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) + (hcard : c * (Real.sqrt D / Real.log D) ≤ ((Nat.sqrtPrimeWindow 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 : ℝ) ^ κ) @@ -421,8 +342,8 @@ theorem falsePositiveCount_le_of_le_card_primesInRange {n D p : ℕ} {κ c : ℝ -- 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 + set numPrimes := (Nat.sqrtPrimeWindow D).card + have hr : 4 ≤ r := Real.four_le_sqrt_natCast 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) @@ -459,7 +380,7 @@ theorem exists_falsePositiveCount_le : 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 + obtain ⟨c, hc, hcard⟩ := Nat.exists_pos_mul_sqrt_div_log_le_card_sqrtPrimeWindow 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) @@ -471,7 +392,8 @@ theorem exists_falsePositiveCount_le : 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 + (falsePositiveCount_le_of_le_card_sqrtPrimeWindow 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 @@ -481,20 +403,20 @@ theorem exists_falsePositiveCount_le : /-- 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 +noncomputable def falsePositiveConst : ℝ := Classical.choose exists_falsePositiveCount_le -/-- The universal false-positive constant `Hashing.falsePositiveConst` is at least one. +/-- The universal false-positive constant `falsePositiveConst` is at least one. -/ -theorem Hashing.one_le_falsePositiveConst : 1 ≤ Hashing.falsePositiveConst := +theorem one_le_falsePositiveConst : 1 ≤ 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`. +false-positive count by `falsePositiveConst * κ * n^3 * log n / sqrt D`. -/ -theorem Hashing.falsePositiveCount_le_const_mul {n D p : ℕ} {κ : ℝ} (hD : 16 ≤ D) (hDn : D ≤ n) +theorem 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 * + ((falsePositiveCount wAB wBC wAC) p : ℝ) ≤ falsePositiveConst * (κ * (n : ℝ) ^ 3 * Real.log n / Real.sqrt D) := (Classical.choose_spec exists_falsePositiveCount_le).2 hD hDn hκ wAB wBC wAC hT hp diff --git a/Complexity/Algorithms/ExactTriangle/Lopsided.lean b/Complexity/Algorithms/ExactTriangle/Lopsided.lean index d52e779437..3b2dbc0b28 100644 --- a/Complexity/Algorithms/ExactTriangle/Lopsided.lean +++ b/Complexity/Algorithms/ExactTriangle/Lopsided.lean @@ -31,11 +31,12 @@ fast solver or subquadratic running time. * `LopsidedGraph.matX_mul_matY_apply`: entries count common middle neighbors. * `LopsidedGraph.matX_mul_matY_apply_ne_zero_iff`: nonzero entries detect common neighbors. -* `LopsidedGraph.zero_one`: the converse for zero/one matrices. -* `LopsidedGraph.entries_from_counts`: bounded signed products from `4 * β²` counts. +* `LopsidedGraph.mul_apply_eq_commonNeighborCount_ofMatrices`: the converse for zero/one matrices. +* `LopsidedGraph.mul_apply_eq_sum_signedBitInstance_counts`: signed products from `4 * β²` counts. ## References +The lopsided matrix interface is used in [AlmanWilliams2026, Section 3]. Adapted from [Anthropic's MatrixLanguage.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/MatrixLanguage.lean) and the lopsided definitions in [PaperStatements.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/PaperStatements.lean#L798-L852). The adaptation uses `Set.EqOn`, pointwise signed parts, and direct finite cardinalities. @@ -73,7 +74,7 @@ noncomputable def commonNeighborCount {n : ℕ} (I : LopsidedGraph n) (a b : Fin (Finset.univ.filter fun v : I.M => I.adjA a v ∧ I.adjB v b).card /-- Two thin matrices describe adjacency where their entries equal one. -The counting theorem `zero_one` requires all entries to be zero or one. -/ +The theorem `mul_apply_eq_commonNeighborCount_ofMatrices` requires zero/one entries. -/ def ofMatrices {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D) (Fin n) ℤ) (W : Finset (Fin n × Fin n)) : LopsidedGraph n where M := Fin D @@ -180,16 +181,15 @@ theorem counting_eqOn_iff (ans : Fin n × Fin n → ℕ) : simp only [Set.EqOn, matX_mul_matY_apply, Nat.cast_inj] /-- Counts on queried edges provide their detection answers. -/ -theorem detection_eqOn_of_counting (ans : Fin n × Fin n → ℕ) +theorem detection_iff_of_counting (ans : Fin n × Fin n → ℕ) (h : Set.EqOn ans (fun q => I.commonNeighborCount q.1 q.2) (I.W : Set _)) : - Set.EqOn (fun q => ans q ≠ 0) (fun q => I.HasCommonNeighbor q.1 q.2) (I.W : Set _) := by + ∀ q ∈ I.W, ans q ≠ 0 ↔ I.HasCommonNeighbor q.1 q.2 := by intro q hq - change (ans q ≠ 0) = I.HasCommonNeighbor q.1 q.2 - apply propext rw [h hq, commonNeighborCount_ne_zero_iff] /-- For zero/one matrices, each product entry is the count in the graph they describe. -/ -theorem zero_one {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D) (Fin n) ℤ) +theorem mul_apply_eq_commonNeighborCount_ofMatrices {n D : ℕ} + (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D) (Fin n) ℤ) (W : Finset (Fin n × Fin n)) (hX : ∀ a j, X a j = 0 ∨ X a j = 1) (hY : ∀ j b, Y j b = 0 ∨ Y j b = 1) (a b : Fin n) : (X * Y) a b = ((ofMatrices X Y W).commonNeighborCount a b : ℤ) := by @@ -203,8 +203,8 @@ theorem zero_one {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D Every instance has exactly `D` middle vertices and the original queried edges `W`. -/ def signedBitInstance (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D) (Fin n) ℤ) (W : Finset (Fin n × Fin n)) (negativeX negativeY : Bool) (i j : ℕ) : LopsidedGraph n := - ofMatrices (IntMatrix.bit (if negativeX then (X.map negPart) else (X.map posPart)) i) - (IntMatrix.bit (if negativeY then (Y.map negPart) else (Y.map posPart)) j) W + ofMatrices (Matrix.bit (if negativeX then (X.map negPart) else (X.map posPart)) i) + (Matrix.bit (if negativeY then (Y.map negPart) else (Y.map posPart)) j) W /-- The finite family of sign/digit instances for a magnitude width `β`. Different slots may describe the same graph. -/ @@ -223,18 +223,19 @@ private theorem mul_apply_eq_sum_counts {n D : ℕ} (X : Matrix (Fin n) (Fin D) (Y : Matrix (Fin D) (Fin n) ℤ) (W : Finset (Fin n × Fin n)) (β : ℕ) (hX : ∀ a j, 0 ≤ X a j ∧ X a j < 2 ^ β) (hY : ∀ j b, 0 ≤ Y j b ∧ Y j b < 2 ^ β) (a b : Fin n) : (X * Y) a b = ∑ i ∈ Finset.range β, ∑ j ∈ Finset.range β, (2 : ℤ) ^ (i + j) * - ((ofMatrices (IntMatrix.bit X i) (IntMatrix.bit Y j) W).commonNeighborCount a b : ℤ) := by - rw [IntMatrix.product X Y β hX hY, Matrix.sum_apply] + ((ofMatrices (Matrix.bit X i) (Matrix.bit Y j) W).commonNeighborCount a b : ℤ) := by + rw [Matrix.mul_eq_sum_two_pow_smul_bit_mul_bit X Y β hX hY, Matrix.sum_apply] refine Finset.sum_congr rfl fun i _ => ?_ rw [Matrix.sum_apply] refine Finset.sum_congr rfl fun j _ => ?_ - rw [Matrix.smul_apply, smul_eq_mul, zero_one _ _ W (IntMatrix.bit_eq_zero_or_one X i) - (IntMatrix.bit_eq_zero_or_one Y j)] + rw [Matrix.smul_apply, smul_eq_mul, mul_apply_eq_commonNeighborCount_ofMatrices _ _ W + (Matrix.bit_eq_zero_or_one X i) + (Matrix.bit_eq_zero_or_one Y j)] /-- Signed bounded product entries are weighted sums of four count instances per digit pair, with signs `+, −, −, +`. Only counts at the same queried edge are needed. The absolute bounds are strict; `β = 0` therefore permits only zero entries. -/ -theorem entries_from_counts {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) +theorem mul_apply_eq_sum_signedBitInstance_counts {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) (Y : Matrix (Fin D) (Fin n) ℤ) (W : Finset (Fin n × Fin n)) (β : ℕ) (hX : ∀ a j, |X a j| < 2 ^ β) (hY : ∀ j b, |Y j b| < 2 ^ β) (a b : Fin n) : (X * Y) a b = ∑ i ∈ Finset.range β, ∑ j ∈ Finset.range β, (2 : ℤ) ^ (i + j) * @@ -247,10 +248,10 @@ theorem entries_from_counts {n D : ℕ} (X : Matrix (Fin n) (Fin D) ℤ) + ((signedBitInstance X Y W true true i j).commonNeighborCount a b : ℤ)) := by simp only [signedBitInstance, Bool.false_eq_true, ↓reduceIte] - obtain ⟨hXpos, hXneg⟩ := IntMatrix.parts_lt X hX - obtain ⟨hYpos, hYneg⟩ := IntMatrix.parts_lt Y hY + obtain ⟨hXpos, hXneg⟩ := Matrix.map_posPart_lt_and_map_negPart_lt X hX + obtain ⟨hYpos, hYneg⟩ := Matrix.map_posPart_lt_and_map_negPart_lt Y hY -- `XY` is a signed sum of the four products of a part of `X` and a part of `Y`. - conv_lhs => rw [IntMatrix.parts X, IntMatrix.parts Y] + conv_lhs => rw [Matrix.eq_map_posPart_sub_map_negPart X, Matrix.eq_map_posPart_sub_map_negPart Y] rw [Matrix.sub_mul, Matrix.mul_sub, Matrix.mul_sub] simp only [Matrix.sub_apply] rw [mul_apply_eq_sum_counts _ _ W β hXpos hYpos, diff --git a/Complexity/Algorithms/ExactTriangle/Splitting.lean b/Complexity/Algorithms/ExactTriangle/Splitting.lean index 4f6f435c3a..1433c41827 100644 --- a/Complexity/Algorithms/ExactTriangle/Splitting.lean +++ b/Complexity/Algorithms/ExactTriangle/Splitting.lean @@ -131,23 +131,15 @@ def totalChunks {n : ℕ} (wAB : Fin n → Fin n → ℤ) (D p : ℕ) : ℕ := variable {n : ℕ} -/-- `√D ≤ D ≤ n`. -/ -theorem sqrt_le_natCast {D : ℕ} (hD : 1 ≤ D) (hDn : D ≤ n) : Real.sqrt D ≤ n := - calc Real.sqrt D ≤ Real.sqrt D * Real.sqrt D := - le_mul_of_one_le_left (Real.sqrt_nonneg _) (Real.one_le_sqrt.mpr (by exact_mod_cast hD)) - _ = D := Real.mul_self_sqrt (Nat.cast_nonneg D) - _ ≤ n := by exact_mod_cast hDn - /-- The bound `⌊n²/√D⌋` on the size of a chunk is at least 1, because `√D ≤ n ≤ n²`. -/ theorem one_le_queryCap {D : ℕ} (hD : 1 ≤ D) (hDn : D ≤ n) : 1 ≤ queryCap n D := by rw [queryCap_eq_floor] have hn : (1 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn refine Nat.le_floor ?_ rw [Nat.cast_one, le_div_iff₀ (Real.sqrt_pos.mpr (by exact_mod_cast hD)), one_mul] - exact (sqrt_le_natCast hD hDn).trans (le_self_pow₀ hn two_ne_zero) - -private theorem four_le_sqrt {D : ℕ} (hD : 16 ≤ D) : 4 ≤ Real.sqrt (D : ℝ) := - Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) + have hsqrt : Real.sqrt (D : ℝ) ≤ D := + Real.sqrt_le_self_iff.mpr (Or.inr (by exact_mod_cast hD)) + exact (hsqrt.trans (by exact_mod_cast hDn)).trans (le_self_pow₀ hn two_ne_zero) variable {n D g p : ℕ} @@ -163,7 +155,7 @@ theorem sqrt_lt_sOf_add_one (D : ℕ) : Real.sqrt D < (sOf D : ℝ) + 1 := /-- `s > 0` because `D ≥ 16`. -/ private theorem sOf_pos (hD : 16 ≤ D) : (0 : ℝ) < sOf D := by - linarith [four_le_sqrt hD, sqrt_lt_sOf_add_one D] + linarith [Real.four_le_sqrt_natCast hD, sqrt_lt_sOf_add_one D] /-- The bound `⌈s/g⌉` on the size of a piece is at least 1. -/ theorem pieceSize_pos (hD : 16 ≤ D) (hg1 : 1 ≤ g) : 0 < pieceSize D g := by @@ -270,7 +262,7 @@ private theorem sq_lt_mul_queryCap_add (hD : 16 ≤ D) (hDn : D ≤ n) set r := Real.sqrt D with hr set δ := (sOf D : ℝ) + 1 - r with hδ set x := (n : ℝ) ^ 2 / r with hx - have hr4 : 4 ≤ r := four_le_sqrt hD + have hr4 : 4 ≤ r := Real.four_le_sqrt_natCast hD have hxr : x * r = (n : ℝ) ^ 2 := div_mul_cancel₀ _ (by linarith) have hx3 : r ^ 3 ≤ x := by rw [hx, le_div_iff₀ (by linarith)] @@ -304,8 +296,8 @@ private theorem sq_lt_mul_queryCap_add (hD : 16 ≤ D) (hDn : D ≤ n) because D ≥ 16, and √D ≤ √n ≤ 2n/3 because D ≤ n and we may assume n ≥ 3". -/ theorem le_sOf_and_sqrt_le (hD : 16 ≤ D) (hDn : D ≤ n) : 3 * Real.sqrt D / 4 ≤ (sOf D : ℝ) ∧ Real.sqrt D ≤ 2 * (n : ℝ) / 3 := by - have hsqrtD := four_le_sqrt hD - have hsqrtn := four_le_sqrt (hD.trans hDn) + have hsqrtD := Real.four_le_sqrt_natCast hD + have hsqrtn := Real.four_le_sqrt_natCast (hD.trans hDn) constructor · calc 3 * Real.sqrt D / 4 ≤ Real.sqrt D - 1 := by linarith _ ≤ (sOf D : ℝ) := by linarith [sqrt_lt_sOf_add_one D] @@ -385,10 +377,10 @@ private theorem totalChunks_le_add_sOf (hD : 16 ≤ D) (hDn : D ≤ n) (hp0 : p NOTE. A chunk holds a whole number of pairs, at most `⌊n²/√D⌋`, so the count is at most `p + (n² − p)/⌊n²/√D⌋`, and the second term can exceed `√D` slightly. The printed bound is still true, because the count is an integer, `D` is an integer and `D ≤ n` (`sq_lt_mul_queryCap_add`). -/ -theorem totalChunks_le (hD : 16 ≤ D) (hDn : D ≤ n) (hp : p ∈ primesInRange D) : +theorem totalChunks_le (hD : 16 ≤ D) (hDn : D ≤ n) (hp : p ∈ Nat.sqrtPrimeWindow D) : (totalChunks wAB D p : ℝ) ≤ (p : ℝ) + Real.sqrt D ∧ (p : ℝ) + Real.sqrt D ≤ 2 * Real.sqrt D := by - obtain ⟨hprime, hge, hlt⟩ := mem_primesInRange.mp hp + obtain ⟨hprime, hge, hlt⟩ := Nat.mem_sqrtPrimeWindow.mp hp have hnat : (totalChunks wAB D p : ℝ) ≤ (p : ℝ) + (sOf D : ℝ) := by exact_mod_cast totalChunks_le_add_sOf wAB hD hDn hprime.ne_zero hge exact ⟨by linarith [sOf_le_sqrt D], by linarith⟩ @@ -397,9 +389,9 @@ theorem totalChunks_le (hD : 16 ≤ D) (hDn : D ≤ n) (hp : p ∈ primesInRange /-- Proof of Theorem 17: "middle part C_k × ℤ_p, of size at most sp ≤ D". So each instance is an instance of Lop-AE-SparseTri(n, D). -/ -theorem card_middle_lopInstance (hD : 16 ≤ D) (hg1 : 1 ≤ g) (hp : p ∈ primesInRange D) +theorem card_middle_lopInstance (hD : 16 ≤ D) (hg1 : 1 ≤ g) (hp : p ∈ Nat.sqrtPrimeWindow D) (ι : InstanceIndex p) : Fintype.card (lopInstance wAB wBC wAC D g p ι).M ≤ D := by - obtain ⟨-, -, hlt⟩ := mem_primesInRange.mp hp + obtain ⟨-, -, hlt⟩ := Nat.mem_sqrtPrimeWindow.mp hp have hsp : (sOf D : ℝ) * (p : ℝ) ≤ D := calc (sOf D : ℝ) * (p : ℝ) ≤ Real.sqrt D * Real.sqrt D := mul_le_mul (sOf_le_sqrt D) hlt.le (Nat.cast_nonneg p) (Real.sqrt_nonneg _) @@ -475,7 +467,7 @@ theorem card_instanceIndices (D g p : ℕ) : /-- Proof of Theorem 17: "There are at most 2√D h ≤ 2√D(ng/s + 1) ≤ 4ng instances". -/ theorem card_instanceIndices_le (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) - (hp : p ∈ primesInRange D) : + (hp : p ∈ Nat.sqrtPrimeWindow D) : ((instanceIndices wAB D g p).card : ℝ) ≤ 4 * (n : ℝ) * (g : ℝ) := by obtain ⟨hs, hsqrt⟩ := le_sOf_and_sqrt_le hD hDn obtain ⟨hchunks, hchunks'⟩ := totalChunks_le wAB hD hDn hp diff --git a/Complexity/Algorithms/ExactTriangle/Witnesses.lean b/Complexity/Algorithms/ExactTriangle/Witnesses.lean index dfaeebb6aa..02c8cba7a2 100644 --- a/Complexity/Algorithms/ExactTriangle/Witnesses.lean +++ b/Complexity/Algorithms/ExactTriangle/Witnesses.lean @@ -275,26 +275,29 @@ theorem pieceSize_le_two_mul_sqrt_div (hg1 : 1 ≤ g) (hg : (g : ℝ) ≤ Real.s /-- For `16 ≤ D ≤ n`, `1 ≤ g ≤ sqrt D`, `κ ≥ 1` and edge weights bounded by `n^κ`, bound `(falsePositiveCount + 1) * pieceSize` at any selected prime by -`(2 * Hashing.falsePositiveConst + 2) * κ * n^3 * log n / g`. +`(2 * falsePositiveConst + 2) * κ * n^3 * log n / g`. This arithmetic estimate counts candidate vertices rather than native instructions. -/ theorem falsePositiveCount_add_one_mul_pieceSize_le {κ : ℝ} (wAB wBC wAC : Fin n → Fin n → ℤ) (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) (hg : (g : ℝ) ≤ Real.sqrt D) (hκ : 1 ≤ κ) (hT : WeightsBoundedBy wAB wBC wAC ((n : ℝ) ^ κ)) (hp : IsSelectedPrime wAB wBC wAC D p) : ((falsePositiveCount wAB wBC wAC p : ℝ) + 1) * (pieceSize D g : ℝ) - ≤ (2 * Hashing.falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by - have hsqrt : 4 ≤ Real.sqrt (D : ℝ) := - Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) + ≤ (2 * falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by + have hsqrt : 4 ≤ Real.sqrt (D : ℝ) := Real.four_le_sqrt_natCast hD have hn : (16 : ℝ) ≤ n := by exact_mod_cast hD.trans hDn have hlog : 1 ≤ Real.log n := (Real.le_log_iff_exp_le (by linarith)).2 (Real.exp_one_lt_three.le.trans (by linarith)) set B := κ * (n : ℝ) ^ 3 * Real.log n with hB -- The two terms: `F(p) √D = O(B)` by the bound on `F(p)`, and `√D ≤ n ≤ n³ ≤ B`. have hfalse : (falsePositiveCount wAB wBC wAC p : ℝ) * Real.sqrt D ≤ - Hashing.falsePositiveConst * B := by - have hF := Hashing.falsePositiveCount_le_const_mul hD hDn hκ wAB wBC wAC hT hp + falsePositiveConst * B := by + have hF := falsePositiveCount_le_const_mul hD hDn hκ wAB wBC wAC hT hp rwa [← mul_div_assoc, le_div_iff₀ (by linarith)] at hF have hone : Real.sqrt D ≤ B := - calc Real.sqrt D ≤ n := sqrt_le_natCast (by omega) hDn + calc + Real.sqrt D ≤ D := by + apply Real.sqrt_le_self_iff.mpr + exact Or.inr (by exact_mod_cast (show 1 ≤ D by omega)) + _ ≤ n := by exact_mod_cast hDn _ ≤ (n : ℝ) ^ 3 := le_self_pow₀ (by linarith) (by norm_num) _ ≤ κ * (n : ℝ) ^ 3 := le_mul_of_one_le_left (by positivity) hκ _ ≤ B := le_mul_of_one_le_right (by positivity) hlog @@ -304,8 +307,8 @@ theorem falsePositiveCount_add_one_mul_pieceSize_le {κ : ℝ} (wAB wBC wAC : Fi exact pieceSize_le_two_mul_sqrt_div hg1 hg _ = (2 * ((falsePositiveCount wAB wBC wAC p : ℝ) * Real.sqrt D) + 2 * Real.sqrt D) / (g : ℝ) := by ring - _ ≤ (2 * (Hashing.falsePositiveConst * B) + 2 * B) / (g : ℝ) := by gcongr - _ = (2 * Hashing.falsePositiveConst + 2) * (B / (g : ℝ)) := by ring + _ ≤ (2 * (falsePositiveConst * B) + 2 * B) / (g : ℝ) := by gcongr + _ = (2 * falsePositiveConst + 2) * (B / (g : ℝ)) := by ring /-! ### The output of the reduction -/ @@ -319,7 +322,7 @@ theorem recovery_correct (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 ≤ g) {L : List (InstanceIndex p × (Fin n × Fin n))} (hL : IsScanOrder wAB wBC wAC D g p ans L) : IsSearchAnswer wAB wBC wAC (runScans wAB wBC wAC D g L).1 ∧ (runScans wAB wBC wAC D g L).2 ≤ falsePositiveCount wAB wBC wAC p + 1 := by - have hp0 : p ≠ 0 := (mem_primesInRange.mp hp.1).1.ne_zero + have hp0 : p ≠ 0 := (Nat.mem_sqrtPrimeWindow.mp hp.1).1.ne_zero obtain ⟨hzero, hnone, hcount⟩ := runScans_spec wAB wBC wAC D g L refine ⟨⟨hzero, fun hout ⟨a, b, c, h0⟩ => ?_⟩, ?_⟩ · -- "A zero triangle is always found": its scan cannot have failed. @@ -363,8 +366,8 @@ theorem answersCorrect_of_counting (wAB wBC wAC : Fin n → Fin n → ℤ) (D g intro ι hι q hq apply propext simp only [decide_eq_true_eq] - exact Iff.of_eq ((lopInstance wAB wBC wAC D g p ι).detection_eqOn_of_counting - (counts ι) (hcounts ι hι) hq) + exact (lopInstance wAB wBC wAC D g p ι).detection_iff_of_counting + (counts ι) (hcounts ι hι) q hq /-- The recovery scan examines at most the uniform false-positive bound times the piece size. This counts candidate vertices only; it is not a charged native running-time theorem. -/ @@ -377,7 +380,7 @@ theorem recovery_vertexBound {κ : ℝ} (hD : 16 ≤ D) (hDn : D ≤ n) (hg1 : 1 {L : List (InstanceIndex p × (Fin n × Fin n))} (hL : IsScanOrder wAB wBC wAC D g p ans L) : IsSearchAnswer wAB wBC wAC (runScans wAB wBC wAC D g L).1 ∧ ((runScans wAB wBC wAC D g L).2 : ℝ) * pieceSize D g ≤ - (2 * Hashing.falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by + (2 * falsePositiveConst + 2) * (κ * (n : ℝ) ^ 3 * Real.log n / (g : ℝ)) := by obtain ⟨hcorrect, hcount⟩ := recovery_correct hD hDn hg1 wAB wBC wAC hp hans hL refine ⟨hcorrect, le_trans ?_ (falsePositiveCount_add_one_mul_pieceSize_le wAB wBC wAC hD hDn hg1 hg hκ hw hp)⟩ diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Batch.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Batch.lean new file mode 100644 index 0000000000..262f3f864e --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Batch.lean @@ -0,0 +1,194 @@ +/- +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.MatrixMultiplication.Schonhage.Levels + +/-! +# Batch correctness of the tensor product + +## Main results + +* `bilinearProduct_arrays_eq_mul`: exact integer matrix multiplication at selected strings. +* `batch_product_spec`: arbitrary output strings with `m` inner levels. + +The batch arrays vanish on strings with the wrong number of inner positions. +The coefficient rules select equal inner indices and the specified outer row and column, +yielding the usual product of integer matrices. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-! ### Cardinalities of index types -/ + +/-- The outer-index type has `N0 L m` elements. -/ +theorem card_outerStr (L m : ℕ) : Fintype.card (OuterStr L m) = N0 L m := by + simp [N0] + +/-- The inner-index type has `D m` elements. -/ +theorem card_innerStr (m : ℕ) : Fintype.card (InnerStr m) = D m := by + simp [D] + +/-- There are `K L m` subsets of `Fin L` with cardinality `m`. -/ +theorem card_subsets_eq_K (L m : ℕ) : + (univ.filter fun Q : Finset (Fin L) => Q.card = m).card = K L m := by + rw [univ_filter_card_eq, card_powersetCard, card_univ, Fintype.card_fin, K] + +/-! ### The input arrays -/ + +/-- The entry of the left input array at the string with inner set `Q`, outer part `r` and inner +part `π` is `X_Q[r, π]`. -/ +@[simp] +theorem arrayL_leftStrOf {L m : ℕ} (X : Finset (Fin L) → LeftMat L m) (Q : Finset (Fin L)) + (hQ : Q.card = m) (r : OuterStr L m) (π : InnerStr m) : + arrayL m X (leftStrOf Q hQ r π) = X Q r π := by + have h : (innerSetL (leftStrOf Q hQ r π)).card = m := by rw [innerSetL_leftStrOf, hQ] + rw [arrayL, dite_eq_left h, outerPartL_leftStrOf, innerPartL_leftStrOf, innerSetL_leftStrOf] + +/-- The entry of the right input array at the string with inner set `Q`, inner part `π` and outer +part `c` is `Y_Q[π, c]`. -/ +@[simp] +theorem arrayR_rightStrOf {L m : ℕ} (Y : Finset (Fin L) → RightMat L m) (Q : Finset (Fin L)) + (hQ : Q.card = m) (π : InnerStr m) (c : OuterStr L m) : + arrayR m Y (rightStrOf Q hQ π c) = Y Q π c := by + have h : (innerSetR (rightStrOf Q hQ π c)).card = m := by rw [innerSetR_rightStrOf, hQ] + rw [arrayR, dite_eq_left h, innerPartR_rightStrOf, outerPartR_rightStrOf, innerSetR_rightStrOf] + +/-! ### Coefficients and batch correctness -/ + +/-- The inner-output coefficient selects equal inner index pairs (Section 2.3.3). +Checked directly on `coefficient`; no polynomial coefficient bridge is assumed. -/ +theorem coefficient_z0 (s : LeftVar) (t : RightVar) : + coefficient s t .z0 = if ∃ i j, s = .p i j ∧ t = .q i j then 1 else 0 := by + decide +revert + +/-- With two outer inputs, the outer-output coefficient selects exactly the matching row +and column (Section 2.3.3). -/ +theorem coefficient_x_y_z (i' j' i j : Fin 3) : + coefficient (.x i') (.y j') (.z i j) = if i' = i ∧ j' = j then 1 else 0 := by + decide +revert + +/-- The statement on `γ(x_i', y_j', z_ij)` as it is used: if `s` and `t` are outer and +`γ(s, t, z_ij) ≠ 0`, then `(s, t) = (x_i, y_j)`. -/ +private theorem eq_x_y_of_gamma_ne_zero {s : LeftVar} {t : RightVar} {i j : Fin 3} + (hs : ¬ s.IsInner) (ht : ¬ t.IsInner) (h : coefficient s t (.z i j) ≠ 0) : s = .x i ∧ t = .y j + := by + rcases s with i' | _ + · rcases t with j' | _ + · by_contra hne + exact h ((coefficient_x_y_z i' j' i j).trans + (ite_eq_right fun ⟨hi, hj⟩ => hne ⟨hi ▸ rfl, hj ▸ rfl⟩)) + · exact (ht trivial).elim + · exact (hs trivial).elim + +/-- The summands of Lemma 8 that can be nonzero. Let `w` be the output string with inner +set `Q`, row `r` and column `c`. If `γ(u_ℓ, v_ℓ, w_ℓ) ≠ 0` at every level and the inner sets of `u` +and `v` have `m` elements, then `u` has inner set `Q`, outer part `r` and some inner part `π`, and +`v` has inner set `Q`, outer part `c` and the inner part with the same indices. -/ +private theorem exists_of_summand_ne_zero {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card = m) + (r c : OuterStr L m) {u : LeftStr L} {v : RightStr L} + (hγ : ∀ ℓ, coefficient (u ℓ) (v ℓ) (outStrOf Q hQ r c ℓ) ≠ 0) (hu : (innerSetL u).card = m) + (hv : (innerSetR v).card = m) : ∃ π, u = leftStrOf Q hQ r π ∧ v = rightStrOf Q hQ π c := by + simp only [outStrOf, forall_level_iff Q hQ, glue_innerLevel, glue_outerLevel] at hγ + obtain ⟨hγin, hγout⟩ := hγ + -- at a level ℓ ∈ Q we have w_ℓ = z₀ […] So (u_ℓ, v_ℓ) = (p_ij, q_ij) for some i, j. + have hin : ∀ k, ∃ ij : Fin 2 × Fin 2, + u (Finset.orderEmbOfFin Q hQ k) = .p ij.1 ij.2 ∧ v (Finset.orderEmbOfFin Q hQ k) = .q ij.1 + ij.2 := fun k => by + by_contra hne + exact hγin k ((coefficient_z0 _ _).trans (ite_eq_right fun ⟨i, j, h⟩ => hne ⟨(i, j), h⟩)) + choose π hπ using hin + -- Thus u and v have an inner variable at every level of Q, and their inner sets are therefore + -- exactly Q. So u_ℓ and v_ℓ are outer at the levels outside Q. + have huQ : innerSetL u = Q := + filter_eq_of_forall_innerLevel Q hQ (fun k => by rw [(hπ k).1]; trivial) hu + have hvQ : innerSetR v = Q := + filter_eq_of_forall_innerLevel Q hQ (fun k => by rw [(hπ k).2]; trivial) hv + have huout := ((filter_eq_iff_forall_level Q hQ).mp huQ).2 + have hvout := ((filter_eq_iff_forall_level Q hQ).mp hvQ).2 + -- at a level ℓ ∉ Q we have w_ℓ = z_ij for some i, j […] so (u_ℓ, v_ℓ) = (x_i, y_j). + have hout : ∀ k, u ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k) = .x (r k) ∧ v + ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k) = .y (c k) := + fun k => eq_x_y_of_gamma_ne_zero (huout k) (hvout k) (hγout k) + exact ⟨π, (eq_glue_iff Q hQ).mpr ⟨fun k => (hπ k).1, fun k => (hout k).1⟩, + (eq_glue_iff Q hQ).mpr ⟨fun k => (hπ k).2, fun k => (hout k).2⟩⟩ + +/-- In the proof of Lemma 9: the coefficient of the summand of `π` is 1. -/ +private theorem prod_gamma_eq_one {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card = m) + (r c : OuterStr L m) (π : InnerStr m) : + ∏ ℓ, coefficient (leftStrOf Q hQ r π ℓ) (rightStrOf Q hQ π c ℓ) (outStrOf Q hQ r c ℓ) = 1 := by + have hone : ∀ ℓ, + coefficient (leftStrOf Q hQ r π ℓ) (rightStrOf Q hQ π c ℓ) (outStrOf Q hQ r c ℓ) = 1 := by + simp only [forall_level_iff Q hQ, leftStrOf, rightStrOf, outStrOf, glue_innerLevel, + glue_outerLevel] + exact ⟨fun _ => (coefficient_z0 _ _).trans (ite_eq_left ⟨_, _, rfl, rfl⟩), + fun _ => (coefficient_x_y_z _ _ _ _).trans (ite_eq_left ⟨rfl, rfl⟩)⟩ + exact prod_eq_one fun ℓ _ => hone ℓ + +/-- At the output string for a size-`m` set `Q` and indices `r`, `c`, the bilinear product of +the batch arrays equals `(X Q * Y Q) r c`. The matrices have arbitrary integer entries +(cf. [AlmanWilliams2026, Lemma 9]). -/ +theorem bilinearProduct_arrays_eq_mul {L m : ℕ} (X : Finset (Fin L) → LeftMat L m) + (Y : Finset (Fin L) → RightMat L m) (Q : Finset (Fin L)) (hQ : Q.card = m) + (r c : OuterStr L m) : + bilinearProduct (arrayL m X) (arrayR m Y) (outStrOf Q hQ r c) = (X Q * Y Q) r c := by + -- To show: `∑_{(u,v)} ∏_ℓ γ(u_ℓ, v_ℓ, w_ℓ) a[u] b[v] = ∑_π X_Q[r, π] Y_Q[π, c]`. The right-hand + -- side is the part of the left-hand side over the pairs `(u, v)` that belong to some `π`. + rw [tensor_expansion, Matrix.mul_apply, ← Fintype.sum_prod_type'] + symm + refine Fintype.sum_of_injective (fun π => (leftStrOf Q hQ r π, rightStrOf Q hQ π c)) ?_ _ _ ?_ ?_ + · -- There is one such summand for each of the D strings π. + exact fun π π' h => (leftStrOf_inj.mp (congrArg Prod.fst h)).2 + · -- Every other summand is zero. a[u] b[v] ≠ 0 requires the inner sets of u and v to have + -- exactly m elements. + rintro ⟨u, v⟩ hnot + by_contra hne + obtain ⟨⟨hγ, ha⟩, hb⟩ : ((∏ ℓ, coefficient (u ℓ) (v ℓ) (outStrOf Q hQ r c ℓ)) ≠ 0 + ∧ arrayL m X u ≠ 0) ∧ arrayR m Y v ≠ 0 := by + simpa only [mul_ne_zero_iff] using hne + obtain ⟨π, rfl, rfl⟩ := exists_of_summand_ne_zero Q hQ r c + (fun ℓ => prod_ne_zero_iff.mp hγ ℓ (mem_univ ℓ)) + (by_contra fun h => ha (dite_eq_right h)) (by_contra fun h => hb (dite_eq_right h)) + exact hnot ⟨π, rfl⟩ + · -- The summand of `π` is `1 · X_Q[r, π] · Y_Q[π, c]`. + intro π + rw [prod_gamma_eq_one, one_mul, arrayL_leftStrOf, arrayR_rightStrOf] + +/-- An output string with exactly `m` inner levels selects the product of the two matrices +indexed by its inner set, at its extracted row and column. This is an equality of integer +values for the arrays `arrayL m X` and `arrayR m Y` (cf. [AlmanWilliams2026, Lemma 9]). -/ +theorem batch_product_spec {L m : ℕ} (X : Finset (Fin L) → LeftMat L m) (Y : Finset (Fin L) → + RightMat L m) + (w : OutStr L) (h : (innerSetO w).card = m) : + bilinearProduct (arrayL m X) (arrayR m Y) w + = (X (innerSetO w) * Y (innerSetO w)) (rowO w h) (colO w h) := by + -- `w` is the output string with the inner set, the row and the column of `w`. + conv_lhs => rw [outStr_eq_outStrOf h rfl h] + exact bilinearProduct_arrays_eq_mul X Y (innerSetO w) h (rowO w h) (colO w h) + +/-- `Full.result` on the batch arrays agrees with the matrix product for every size-`m` set +`Q` and every row and column index. -/ +theorem Full.arrays_eq_mul {L m : ℕ} (X : Finset (Fin L) → LeftMat L m) + (Y : Finset (Fin L) → RightMat L m) (Q : Finset (Fin L)) (hQ : Q.card = m) + (r c : OuterStr L m) : + Full.result L (arrayL m X) (arrayR m Y) (outStrOf Q hQ r c) = (X Q * Y Q) r c := by + rw [Full.result_eq_bilinearProduct] + exact bilinearProduct_arrays_eq_mul X Y Q hQ r c + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Layout.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Layout.lean new file mode 100644 index 0000000000..684266bc9a --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Layout.lean @@ -0,0 +1,145 @@ +/- +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.MatrixMultiplication.Schonhage.Levels +public import Mathlib.Data.Nat.Sqrt + +/-! +# Tiling layout data + +## Main definitions + +* `Layout`: index equivalences and an injective assignment of level sets to block pairs. +* `wantedStrings`, `tiles`: requested outputs grouped by tile. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The natural square root of the number `K L m` of size-`m` level sets. -/ +def K0 (L m : ℕ) : ℕ := Nat.sqrt (K L m) + +/-- `N₀ = 3^{L-m}` is positive. -/ +theorem N0_pos (L m : ℕ) : 0 < N0 L m := Nat.pow_pos (by decide) + +/-- `K₀ = ⌊√binom(L, m)⌋` is positive when `m ≤ L`. -/ +theorem K0_pos {L m : ℕ} (h : m ≤ L) : 0 < K0 L m := Nat.sqrt_pos.mpr (Nat.choose_pos h) + +/-- Indexing data for the block decomposition, under `m ≤ L`. + +The three equivalences identify block rows, block columns and the thin dimension with +their string indices. The injective `table` assigns a size-`m` subset of `Fin L` to each +pair of block indices in a tile. All tiles use this same table. The structure contains no +matrix values or solver certificate; the correctness results hold for any such indexing +choices. -/ +structure Layout (L m : ℕ) where + hmL : m ≤ L + /-- Index a block's rows by outer strings. -/ + rowIdx : Fin (N0 L m) ≃ OuterStr L m + /-- Index a block's columns by outer strings. -/ + colIdx : Fin (N0 L m) ≃ OuterStr L m + /-- Index the thin matrix dimension by strings of inner index pairs. -/ + innerIdx : Fin (D m) ≃ InnerStr m + /-- Assign a distinct size-`m` level set to each block product in a tile. -/ + table : Fin (K0 L m) × Fin (K0 L m) → Finset (Fin L) + table_card : ∀ gh, (table gh).card = m + table_injective : Function.Injective table + +/-- The common band dimension `K0 L m * N0 L m`: each band contains `K0 L m` consecutive +blocks of size `N0 L m`. -/ +def bandSize (L m : ℕ) : ℕ := K0 L m * N0 L m + +/-- The quotient `(N + bandSize L m - 1) / bandSize L m`. When `m ≤ L`, the band size is +positive and this is the number of bands needed to cover `N` positions, including `N = 0`. -/ +def numBands (L m N : ℕ) : ℕ := (N + bandSize L m - 1) / bandSize L m + +/-- The dimension covered by `numBands L m N` complete bands. -/ +def padN (L m N : ℕ) : ℕ := numBands L m N * bandSize L m + +/-- Extend `X` to natural-number row indices by returning zero at indices `I ≥ N`. -/ +def padRows {N D' : ℕ} (X : Matrix (Fin N) (Fin D') ℤ) (I : ℕ) (k : Fin D') : ℤ := + if h : I < N then X ⟨I, h⟩ k else 0 + +/-- Extend `Y` to natural-number column indices by returning zero at indices `J ≥ N`. -/ +def padCols {N D' : ℕ} (Y : Matrix (Fin D') (Fin N) ℤ) (k : Fin D') (J : ℕ) : ℤ := + if h : J < N then Y k ⟨J, h⟩ else 0 + +/-- Extract the `N0 L m` rows starting at `B * N0 L m` from `padRows X`, using `lay.rowIdx` +and `lay.innerIdx` to index the resulting integer matrix by strings. -/ +def rowBlock {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) (B : ℕ) : + LeftMat L m := + fun r π => padRows X (B * N0 L m + lay.rowIdx.symm r) (lay.innerIdx.symm π) + +/-- Extract the `N0 L m` columns starting at `B * N0 L m` from `padCols Y`, using +`lay.innerIdx` and `lay.colIdx` to index the resulting integer matrix by strings. -/ +def colBlock {L m N : ℕ} (lay : Layout L m) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (B : ℕ) : + RightMat L m := + fun π c => padCols Y (lay.innerIdx.symm π) (B * N0 L m + lay.colIdx.symm c) + +/-- For a table entry `Q = lay.table (g, h)`, select row block `β * K0 L m + g` of `X`; return +the zero matrix for sets outside the table's image. Injectivity of the table ensures that +the defining sum selects at most one matrix. -/ +def bandFamilyL {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) (β : ℕ) : + Finset (Fin L) → LeftMat L m := + fun Q => ∑ gh, if lay.table gh = Q then rowBlock lay X (β * K0 L m + gh.1) else 0 + +/-- For a table entry `Q = lay.table (g, h)`, select column block `β * K0 L m + h` of `Y`; +return the zero matrix for sets outside the table's image. -/ +def bandFamilyR {L m N : ℕ} (lay : Layout L m) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (β : ℕ) : + Finset (Fin L) → RightMat L m := + fun Q => ∑ gh, if lay.table gh = Q then colBlock lay Y (β * K0 L m + gh.2) else 0 + +/-- The left batch array formed from the block family for row band `β`. -/ +def bandArrayL {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) (β : ℕ) : + LeftStr L → ℤ := + arrayL m (bandFamilyL lay X β) + +/-- The right batch array formed from the block family for column band `β`. -/ +def bandArrayR {L m N : ℕ} (lay : Layout L m) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (β : ℕ) : + RightStr L → ℤ := + arrayR m (bandFamilyR lay Y β) + +/-- The band that contains row (or column) `I`. -/ +def bandOf (L m I : ℕ) : ℕ := I / bandSize L m + +/-- The position, within its band, of the block that contains row (or column) `I`. -/ +def blockOf {L m : ℕ} (lay : Layout L m) (I : ℕ) : Fin (K0 L m) := + ⟨I / N0 L m % K0 L m, Nat.mod_lt _ (K0_pos lay.hmL)⟩ + +/-- The position of row (or column) `I` within its block. -/ +def offsetOf (L m I : ℕ) : Fin (N0 L m) := ⟨I % N0 L m, Nat.mod_lt _ (N0_pos L m)⟩ + +/-- Encode the within-tile position `(I, J)` as an output string. Its inner set comes from the +table at the two block indices, and its outer indices come from the two offsets within +those blocks. The tile itself is identified separately by `(bandOf L m I, bandOf L m J)`. -/ +def outStrOfPos {L m : ℕ} (lay : Layout L m) (I J : ℕ) : OutStr L := + outStrOf (lay.table (blockOf lay I, blockOf lay J)) (lay.table_card _) + (lay.rowIdx (offsetOf L m I)) (lay.colIdx (offsetOf L m J)) + +/-- Output strings for the positions of `W` whose row band is `β` and column band is `β'`. +This is the image of the filtered position set under `outStrOfPos lay`. -/ +def wantedStrings {L m N : ℕ} (lay : Layout L m) (W : Finset (Fin N × Fin N)) (β β' : ℕ) : + Finset (OutStr L) := + (W.filter fun IJ => bandOf L m IJ.1 = β ∧ bandOf L m IJ.2 = β').image + fun IJ => outStrOfPos lay IJ.1 IJ.2 + +/-- All the tiles: pairs of a row band and a column band. -/ +def tiles (L m N : ℕ) : Finset (ℕ × ℕ) := range (numBands L m N) ×ˢ range (numBands L m N) + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Levels.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Levels.lean new file mode 100644 index 0000000000..da18d521e0 --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Levels.lean @@ -0,0 +1,322 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Full +public import Complexity.Foundations.Data.Finset.OrderGlue + +/-! +# Matrix entries indexed by tensor strings + +## Main definitions + +* `arrayL`, `arrayR`: batch arrays with exactly `m` inner levels. +* `leftStrOf`, `rightStrOf`, `outStrOf`: entry strings. + +## Main results + +* `existsUnique_leftStr`, `existsUnique_rightStr`, `existsUnique_outStr`: unique entry strings. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- Recognize the constructors `.p i j` as inner left variables. -/ +def LeftVar.IsInner : LeftVar → Prop + | .x _ => False + | .p _ _ => True +/-- Recognize the constructors `.q i j` as inner right variables. -/ +def RightVar.IsInner : RightVar → Prop + | .y _ => False + | .q _ _ => True +instance LeftVar.instDecidablePredIsInner : DecidablePred LeftVar.IsInner + | .x _ => isFalse fun h => h + | .p _ _ => isTrue trivial +instance RightVar.instDecidablePredIsInner : DecidablePred RightVar.IsInner + | .y _ => isFalse fun h => h + | .q _ _ => isTrue trivial +/-- Positions of the inner left variables in `u`. -/ +def innerSetL {L : ℕ} (u : LeftStr L) : Finset (Fin L) := univ.filter fun ℓ => (u ℓ).IsInner +/-- Positions of the inner right variables in `v`. -/ +def innerSetR {L : ℕ} (v : RightStr L) : Finset (Fin L) := univ.filter fun ℓ => (v ℓ).IsInner +/-- Positions at which the output string is `.z0`. -/ +def innerSetO {L : ℕ} (w : OutStr L) : Finset (Fin L) := univ.filter fun ℓ => (w ℓ).IsInner + +/-- The outer dimension `3 ^ (L - m)`, with natural-number subtraction. In particular, this +definition gives `1` if `m > L`. -/ +def N0 (L m : ℕ) : ℕ := 3 ^ (L - m) +/-- The thin dimension `4 ^ m`. -/ +def D (m : ℕ) : ℕ := 4 ^ m +/-- The binomial coefficient `L.choose m`, counting size-`m` subsets of `Fin L`. -/ +def K (L m : ℕ) : ℕ := L.choose m + +/-- Extract `i` from `.x i`; return `0` on inner left variables. -/ +def LeftVar.outerIndex : LeftVar → Fin 3 + | .x i => i + | .p _ _ => 0 +/-- Extract `(i, j)` from `.p i j`; return `(0, 0)` on outer left variables. -/ +def LeftVar.innerIndex : LeftVar → Fin 2 × Fin 2 + | .x _ => (0, 0) + | .p i j => (i, j) +/-- Extract `j` from `.y j`; return `0` on inner right variables. -/ +def RightVar.outerIndex : RightVar → Fin 3 + | .y j => j + | .q _ _ => 0 +/-- Extract `(i, j)` from `.q i j`; return `(0, 0)` on outer right variables. -/ +def RightVar.innerIndex : RightVar → Fin 2 × Fin 2 + | .y _ => (0, 0) + | .q i j => (i, j) +/-- Extract the first index from `.z i j`; return `0` on `.z0`. -/ +def OutVar.rowIndex : OutVar → Fin 3 + | .z i _ => i + | .z0 => 0 +/-- Extract the second index from `.z i j`; return `0` on `.z0`. -/ +def OutVar.colIndex : OutVar → Fin 3 + | .z _ j => j + | .z0 => 0 + +/-- Functions from `Fin (L - m)` to `Fin 3`, used as outer row and column indices. Each value +specifies the index of an outer left or right variable. -/ +abbrev OuterStr (L m : ℕ) : Type := Fin (L - m) → Fin 3 + +/-- Functions from `Fin m` to pairs of `Fin 2` indices, used for the common inner dimension of +the two matrices. Left and right inner variables with the same pair of indices share one +element of this type, so contraction over `InnerStr m` is ordinary matrix multiplication. -/ +abbrev InnerStr (m : ℕ) : Type := Fin m → Fin 2 × Fin 2 + +/-- Extract the outer left indices in increasing order of their positions, assuming `u` has +exactly `m` inner positions. -/ +def outerPartL {L m : ℕ} (u : LeftStr L) (h : (innerSetL u).card = m) : OuterStr L m := + fun k => (u (((innerSetL u)ᶜ.orderEmbOfFin (card_compl_of_card_eq (innerSetL u) h)) k)).outerIndex +/-- Extract the inner left index pairs in increasing order of their positions, assuming `u` +has exactly `m` inner positions. -/ +def innerPartL {L m : ℕ} (u : LeftStr L) (h : (innerSetL u).card = m) : InnerStr m := + fun k => (u (Finset.orderEmbOfFin (innerSetL u) h k)).innerIndex +/-- The outer part of a right string whose inner set has `m` elements. -/ +def outerPartR {L m : ℕ} (v : RightStr L) (h : (innerSetR v).card = m) : OuterStr L m := + fun k => (v (((innerSetR v)ᶜ.orderEmbOfFin (card_compl_of_card_eq (innerSetR v) h)) k)).outerIndex +/-- The inner part of a right string whose inner set has `m` elements. -/ +def innerPartR {L m : ℕ} (v : RightStr L) (h : (innerSetR v).card = m) : InnerStr m := + fun k => (v (Finset.orderEmbOfFin (innerSetR v) h k)).innerIndex + +/-- Read the first indices of the outer output variables in increasing position order, +assuming `w` has exactly `m` inner positions. -/ +def rowO {L m : ℕ} (w : OutStr L) (h : (innerSetO w).card = m) : OuterStr L m := + fun k => (w (((innerSetO w)ᶜ.orderEmbOfFin (card_compl_of_card_eq (innerSetO w) h)) k)).rowIndex +/-- Read the second indices of the outer output variables in increasing position order, +assuming `w` has exactly `m` inner positions. -/ +def colO {L m : ℕ} (w : OutStr L) (h : (innerSetO w).card = m) : OuterStr L m := + fun k => (w (((innerSetO w)ᶜ.orderEmbOfFin (card_compl_of_card_eq (innerSetO w) h)) k)).colIndex + +/-- Integer matrices with outer row indices and inner column indices. -/ +abbrev LeftMat (L m : ℕ) : Type := Matrix (OuterStr L m) (InnerStr m) ℤ +/-- Integer matrices with inner row indices and outer column indices. -/ +abbrev RightMat (L m : ℕ) : Type := Matrix (InnerStr m) (OuterStr L m) ℤ + +/-- Encode a family of integer matrices as a left array. A string with exactly `m` inner +positions selects `X` at that set of positions, with row `outerPartL` and column +`innerPartL`; every other string has value zero. -/ +def arrayL {L : ℕ} (m : ℕ) (X : Finset (Fin L) → LeftMat L m) : LeftStr L → ℤ := + fun u => + if h : (innerSetL u).card = m then X (innerSetL u) (outerPartL u h) (innerPartL u h) else 0 + +/-- Encode a family of integer matrices as a right array. A string with exactly `m` inner +positions selects `Y` at that set of positions, with row `innerPartR` and column +`outerPartR`; every other string has value zero. -/ +def arrayR {L : ℕ} (m : ℕ) (Y : Finset (Fin L) → RightMat L m) : RightStr L → ℤ := + fun v => + if h : (innerSetR v).card = m then Y (innerSetR v) (innerPartR v h) (outerPartR v h) else 0 + +/-! ### The block-count parameter `M` -/ + +/-- The quantity `K L m * N0 L m ^ 2` used in the leaf and tile bounds. -/ +def M (L m : ℕ) : ℕ := K L m * N0 L m ^ 2 + + +variable {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card = m) + +/-! ### Left, right and output strings with a given inner set, outer part and inner part -/ + +/-- The left string with inner set `Q`, outer part `r` and inner part `π`. -/ +def leftStrOf (r : OuterStr L m) (π : InnerStr m) : LeftStr L := + glue Q hQ (fun k => .p (π k).1 (π k).2) (fun k => .x (r k)) + +/-- The right string with inner set `Q`, inner part `π` and outer part `c`. -/ +def rightStrOf (π : InnerStr m) (c : OuterStr L m) : RightStr L := + glue Q hQ (fun k => .q (π k).1 (π k).2) (fun k => .y (c k)) + +/-- Construct an output string with `.z0` at `Q` and `.z (r k) (c k)` at the `k`-th position +of its complement. Both sets of positions use increasing order. -/ +def outStrOf (r c : OuterStr L m) : OutStr L := + glue Q hQ (fun _ => .z0) (fun k => .z (r k) (c k)) + +/-- The inner set of `leftStrOf Q hQ r π` is `Q`. -/ +@[simp] +theorem innerSetL_leftStrOf (r : OuterStr L m) (π : InnerStr m) : + innerSetL (leftStrOf Q hQ r π) = Q := + filter_glue Q hQ LeftVar.IsInner (fun _ => trivial) fun _ => id + +/-- The inner set of `rightStrOf Q hQ π c` is `Q`. -/ +@[simp] +theorem innerSetR_rightStrOf (π : InnerStr m) (c : OuterStr L m) : + innerSetR (rightStrOf Q hQ π c) = Q := + filter_glue Q hQ RightVar.IsInner (fun _ => trivial) fun _ => id + +/-- The inner set of `outStrOf Q hQ r c` is `Q`. -/ +@[simp] +theorem innerSetO_outStrOf (r c : OuterStr L m) : innerSetO (outStrOf Q hQ r c) = Q := + filter_glue Q hQ OutVar.IsInner (fun _ => trivial) fun _ => id + +/-- The outer part of `leftStrOf Q hQ r π` is `r`. -/ +@[simp] +theorem outerPartL_leftStrOf (r : OuterStr L m) (π : InnerStr m) + (h : (innerSetL (leftStrOf Q hQ r π)).card = m) : outerPartL (leftStrOf Q hQ r π) h = r := + funext fun k => congrArg LeftVar.outerIndex + (glue_outerLevel_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k) + +/-- The inner part of `leftStrOf Q hQ r π` is `π`. -/ +@[simp] +theorem innerPartL_leftStrOf (r : OuterStr L m) (π : InnerStr m) + (h : (innerSetL (leftStrOf Q hQ r π)).card = m) : innerPartL (leftStrOf Q hQ r π) h = π := + funext fun k => congrArg LeftVar.innerIndex + (glue_innerLevel_of_eq Q hQ _ _ (innerSetL_leftStrOf Q hQ r π) h k) + +/-- The inner part of `rightStrOf Q hQ π c` is `π`. -/ +@[simp] +theorem innerPartR_rightStrOf (π : InnerStr m) (c : OuterStr L m) + (h : (innerSetR (rightStrOf Q hQ π c)).card = m) : innerPartR (rightStrOf Q hQ π c) h = π := + funext fun k => congrArg RightVar.innerIndex + (glue_innerLevel_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k) + +/-- The outer part of `rightStrOf Q hQ π c` is `c`. -/ +@[simp] +theorem outerPartR_rightStrOf (π : InnerStr m) (c : OuterStr L m) + (h : (innerSetR (rightStrOf Q hQ π c)).card = m) : outerPartR (rightStrOf Q hQ π c) h = c := + funext fun k => congrArg RightVar.outerIndex + (glue_outerLevel_of_eq Q hQ _ _ (innerSetR_rightStrOf Q hQ π c) h k) + +/-- The variable of `outStrOf Q hQ r c` at the `k`-th lowest level outside its inner set. -/ +private theorem outStrOf_outerLevel (r c : OuterStr L m) + (h : (innerSetO (outStrOf Q hQ r c)).card = m) (k : Fin (L - m)) : + outStrOf Q hQ r c (((innerSetO (outStrOf Q hQ r c))ᶜ.orderEmbOfFin (card_compl_of_card_eq + (innerSetO (outStrOf Q hQ r c)) h)) k) = .z (r k) (c k) := by + rw [outerLevel_congr Q hQ (innerSetO_outStrOf Q hQ r c) h, outStrOf, glue_outerLevel] + +/-- The row of `outStrOf Q hQ r c` is `r`. -/ +@[simp] +theorem rowO_outStrOf (r c : OuterStr L m) (h : (innerSetO (outStrOf Q hQ r c)).card = m) : + rowO (outStrOf Q hQ r c) h = r := + funext fun k => congrArg OutVar.rowIndex (outStrOf_outerLevel Q hQ r c h k) + +/-- The column of `outStrOf Q hQ r c` is `c`. -/ +@[simp] +theorem colO_outStrOf (r c : OuterStr L m) (h : (innerSetO (outStrOf Q hQ r c)).card = m) : + colO (outStrOf Q hQ r c) h = c := + funext fun k => congrArg OutVar.colIndex (outStrOf_outerLevel Q hQ r c h k) + +variable {Q} + +/-- A left string with inner set `Q` is the string `leftStrOf` of its outer part and its inner +part. -/ +theorem leftStr_eq_leftStrOf {u : LeftStr L} (hu : innerSetL u = Q) (h : (innerSetL u).card = m) : + u = leftStrOf Q hQ (outerPartL u h) (innerPartL u h) := by + subst hu + refine glue_index _ h LeftVar.IsInner (inner := fun π : Fin 2 × Fin 2 => .p π.1 π.2) ?_ ?_ rfl + · rintro (i | ⟨i, j⟩) hs + exacts [hs.elim, rfl] + · rintro (i | ⟨i, j⟩) hs + exacts [rfl, (hs trivial).elim] + +/-- A right string with inner set `Q` is the string `rightStrOf` of its inner part and its outer +part. -/ +theorem rightStr_eq_rightStrOf {v : RightStr L} (hv : innerSetR v = Q) + (h : (innerSetR v).card = m) : v = rightStrOf Q hQ (innerPartR v h) (outerPartR v h) := by + subst hv + refine glue_index _ h RightVar.IsInner (inner := fun π : Fin 2 × Fin 2 => .q π.1 π.2) ?_ ?_ rfl + · rintro (j | ⟨i, j⟩) ht + exacts [ht.elim, rfl] + · rintro (j | ⟨i, j⟩) ht + exacts [rfl, (ht trivial).elim] + +/-- An output string with inner set `Q` is the string `outStrOf` of its row and its column. -/ +theorem outStr_eq_outStrOf {w : OutStr L} (hw : innerSetO w = Q) (h : (innerSetO w).card = m) : + w = outStrOf Q hQ (rowO w h) (colO w h) := by + subst hw + refine glue_index _ h OutVar.IsInner (inner := fun _ : Unit => .z0) + (outer := fun rc : Fin 3 × Fin 3 => .z rc.1 rc.2) (innerIndex := fun _ => ()) + (outerIndex := fun z => (z.rowIndex, z.colIndex)) ?_ ?_ rfl + · rintro (⟨i, j⟩ | _) hz + exacts [hz.elim, rfl] + · rintro (⟨i, j⟩ | _) hz + exacts [rfl, (hz trivial).elim] + +variable {hQ} + +/-- Left strings with the same inner set and different outer or inner parts are different. -/ +theorem leftStrOf_inj {r r' : OuterStr L m} {π π' : InnerStr m} : + leftStrOf Q hQ r π = leftStrOf Q hQ r' π' ↔ r = r' ∧ π = π' := by + refine ⟨fun h => ?_, by rintro ⟨rfl, rfl⟩; rfl⟩ + obtain ⟨hp, hx⟩ := (glue_inj Q hQ).mp h + exact ⟨funext fun k => LeftVar.x.inj (congrFun hx k), + funext fun k => Prod.ext (LeftVar.p.inj (congrFun hp k)).1 (LeftVar.p.inj (congrFun hp k)).2⟩ + +/-- Output strings with different inner sets, rows or columns are different. -/ +theorem outStrOf_inj {Q' : Finset (Fin L)} {hQ' : Q'.card = m} {r c r' c' : OuterStr L m} : + outStrOf Q hQ r c = outStrOf Q' hQ' r' c' ↔ Q = Q' ∧ r = r' ∧ c = c' := by + refine ⟨fun h => ?_, by rintro ⟨rfl, rfl, rfl⟩; rfl⟩ + obtain rfl : Q = Q' := by rw [← innerSetO_outStrOf Q hQ r c, h, innerSetO_outStrOf] + have hz := fun k => OutVar.z.inj (congrFun ((glue_inj Q hQ).mp h).2 k) + exact ⟨rfl, funext fun k => (hz k).1, funext fun k => (hz k).2⟩ + +variable (Q) (hQ) +include hQ + +/-- For a set `Q` of cardinality `m`, an outer index `r` and an inner index `π` determine a +unique left string with those extracted parts and inner set `Q`. -/ +theorem existsUnique_leftStr (r : OuterStr L m) (π : InnerStr m) : + ∃! u : LeftStr L, innerSetL u = Q + ∧ ∃ h : (innerSetL u).card = m, outerPartL u h = r ∧ innerPartL u h = π := by + have hset := innerSetL_leftStrOf Q hQ r π + refine ⟨leftStrOf Q hQ r π, ⟨hset, hset.symm ▸ hQ, outerPartL_leftStrOf Q hQ r π _, + innerPartL_leftStrOf Q hQ r π _⟩, ?_⟩ + rintro u ⟨hu, h, rfl, rfl⟩ + exact leftStr_eq_leftStrOf hQ hu h + +/-- For a set `Q` of cardinality `m`, an inner index `π` and an outer index `c` determine a +unique right string with those extracted parts and inner set `Q`. -/ +theorem existsUnique_rightStr (π : InnerStr m) (c : OuterStr L m) : + ∃! v : RightStr L, innerSetR v = Q + ∧ ∃ h : (innerSetR v).card = m, innerPartR v h = π ∧ outerPartR v h = c := by + have hset := innerSetR_rightStrOf Q hQ π c + refine ⟨rightStrOf Q hQ π c, ⟨hset, hset.symm ▸ hQ, + innerPartR_rightStrOf Q hQ π c _, outerPartR_rightStrOf Q hQ π c _⟩, ?_⟩ + rintro v ⟨hv, h, rfl, rfl⟩ + exact rightStr_eq_rightStrOf hQ hv h + +/-- For a set `Q` of cardinality `m`, a row `r` and a column `c` determine a unique output +string with those extracted indices and inner set `Q`. -/ +theorem existsUnique_outStr (r c : OuterStr L m) : + ∃! w : OutStr L, innerSetO w = Q + ∧ ∃ h : (innerSetO w).card = m, rowO w h = r ∧ colO w h = c := by + have hset := innerSetO_outStrOf Q hQ r c + refine ⟨outStrOf Q hQ r c, ⟨hset, hset.symm ▸ hQ, rowO_outStrOf Q hQ r c _, + colO_outStrOf Q hQ r c _⟩, ?_⟩ + rintro w ⟨hw, h, rfl, rfl⟩ + exact outStr_eq_outStrOf hQ hw h + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Magnitude.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Magnitude.lean new file mode 100644 index 0000000000..7a932b6803 --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Magnitude.lean @@ -0,0 +1,250 @@ +/- +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.MatrixMultiplication.Schonhage.ThinProduct +public import Complexity.Foundations.Algebra.FiniteMagnitude + +/-! +# Magnitudes of specified integer values and auxiliary counts + +## Main results + +* `abs_bandArrayL_le`, `abs_bandArrayR_le`: input entries are unchanged or zero. +* `abs_encodingL_le`, `abs_encodingR_le`, `Pruned.abs_result_le`: explicit integer magnitude bounds. +* `minutiae_total_le`: bound on layout subsets, tiles and wanted positions. + +The magnitude bounds concern band-array values, encoded values and returned pruned +values. They do not bound every intermediate state of a native implementation or certify +a machine word width. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-! ### Auxiliary finite counts -/ + +/-- `2^{m/9} ≤ 2^m`. -/ +private theorem two_rpow_ninth_le (m : ℕ) : (2 : ℝ) ^ ((m : ℝ) / 9) ≤ (2 : ℝ) ^ m := by + rw [← Real.rpow_natCast] + have hm : (0 : ℝ) ≤ m := Nat.cast_nonneg m + exact Real.rpow_le_rpow_of_exponent_le (by norm_num) (by linarith) + +/-- A number that is at least `2^{m/9}` is at least 1 after multiplication by `2^{-m/9}`. -/ +private theorem one_le_saving_mul {m : ℕ} {x : ℝ} (hx : (2 : ℝ) ^ ((m : ℝ) / 9) ≤ x) : + 1 ≤ saving m * x := by + rw [← saving_mul_two_rpow m] + exact mul_le_mul_of_nonneg_left hx (saving_pos m).le + +/-- When `D m ^ 18 ≤ N`, the number `K (19 * m) m` of size-`m` level sets is at most `N`. -/ +theorem K_le_of_D_pow_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : K (19 * m) m ≤ N := + (Nat.le_mul_of_pos_right _ (N0_pos _ _)).trans (K_mul_N0_le m N hN) + +/-- For `m ≤ L`, the tile count multiplied by `N0 L m` is at most four times the square of the +padded dimension `padN L m N`. -/ +theorem card_tiles_mul_N0_le {L m : ℕ} (hmL : m ≤ L) (N : ℕ) : + (tiles L m N).card * N0 L m ≤ 4 * padN L m N ^ 2 := by + -- `N₀ ≤ K N₀² = M`, and there are at most `4N²/M` tiles. + refine le_trans (Nat.mul_le_mul_left _ ?_) (card_tiles_mul_M_le L m N) + calc N0 L m = 1 * N0 L m ^ 1 := by ring + _ ≤ K L m * N0 L m ^ 2 := + Nat.mul_le_mul (Nat.choose_pos hmL) (Nat.pow_le_pow_right (N0_pos L m) (by norm_num)) + +/-- The size condition `D m ^ 18 ≤ N` implies the real bound `2 ^ (m / 9) ≤ N`. -/ +theorem two_rpow_le_of_D_pow_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : + (2 : ℝ) ^ ((m : ℝ) / 9) ≤ (N : ℝ) := by + -- `2^{m/9} ≤ 2^m ≤ 4^{18m} = D^18 ≤ N`. + have hpow : 2 ^ m ≤ N := + calc 2 ^ m ≤ 4 ^ m := Nat.pow_le_pow_left (by norm_num) m + _ ≤ (4 ^ m) ^ 18 := Nat.le_self_pow (by norm_num) _ + _ ≤ N := hN + exact (two_rpow_ninth_le m).trans (by exact_mod_cast hpow) + +/-- At `L = 19 * m`, the outer dimension is at least the real power `2 ^ (m / 9)`. -/ +theorem two_rpow_le_N0 (m : ℕ) : (2 : ℝ) ^ ((m : ℝ) / 9) ≤ (N0 (19 * m) m : ℝ) := by + -- `2^{m/9} ≤ 2^m ≤ 3^{18m} = N₀`. + have hpow : 2 ^ m ≤ N0 (19 * m) m := + calc 2 ^ m ≤ 3 ^ m := Nat.pow_le_pow_left (by norm_num) m + _ ≤ 3 ^ (19 * m - m) := Nat.pow_le_pow_right (by norm_num) (by omega) + exact (two_rpow_ninth_le m).trans (by exact_mod_cast hpow) + +/-- The three counts of the remaining minutiae (Section 2.4.4), one by one: the subsets in the +table, the tiles, and the positions of `W`. -/ +private theorem minutiae_parts (m N : ℕ) (W : Finset (Fin N × Fin N)) (hN : D m ^ 18 ≤ N) + (hW : W.card * 2 ^ m ≤ N ^ 2) : + (K (19 * m) m : ℝ) ≤ saving m * (N : ℝ) ^ 2 + ∧ ((tiles (19 * m) m N).card : ℝ) ≤ 4 * saving m * (padN (19 * m) m N : ℝ) ^ 2 + ∧ (W.card : ℝ) ≤ saving m * (N : ℝ) ^ 2 := by + have hc := saving_pos m + refine ⟨?_, ?_, ?_⟩ + · -- The subsets in the table: `K ≤ N ≤ 2^{-m/9} N²`. + calc (K (19 * m) m : ℝ) ≤ N := by exact_mod_cast K_le_of_D_pow_le m N hN + _ ≤ saving m * N * N := + le_mul_of_one_le_left N.cast_nonneg (one_le_saving_mul (two_rpow_le_of_D_pow_le m N hN)) + _ = saving m * (N : ℝ) ^ 2 := by ring + · -- The tiles: at most `4N²/N₀ ≤ 4 · 2^{-m/9} N²`. + have hT : ((tiles (19 * m) m N).card : ℝ) * N0 (19 * m) m + ≤ 4 * (padN (19 * m) m N : ℝ) ^ 2 := by + exact_mod_cast card_tiles_mul_N0_le (show m ≤ 19 * m by omega) N + calc ((tiles (19 * m) m N).card : ℝ) + ≤ (tiles (19 * m) m N).card * (saving m * N0 (19 * m) m) := + le_mul_of_one_le_right (Nat.cast_nonneg _) (one_le_saving_mul (two_rpow_le_N0 m)) + _ = saving m * ((tiles (19 * m) m N).card * N0 (19 * m) m) := by ring + _ ≤ saving m * (4 * (padN (19 * m) m N : ℝ) ^ 2) := by gcongr + _ = 4 * saving m * (padN (19 * m) m N : ℝ) ^ 2 := by ring + · -- The positions: `|W| ≤ 2^{-m} N² ≤ 2^{-m/9} N²`. + have hW' : (W.card : ℝ) * (2 : ℝ) ^ m ≤ (N : ℝ) ^ 2 := by exact_mod_cast hW + calc (W.card : ℝ) = saving m * (W.card * (2 : ℝ) ^ ((m : ℝ) / 9)) := by + rw [mul_left_comm, saving_mul_two_rpow, mul_one] + _ ≤ saving m * (W.card * (2 : ℝ) ^ m) := by gcongr; exact two_rpow_ninth_le m + _ ≤ saving m * (N : ℝ) ^ 2 := by gcongr + +/-- At `L = 19 * m`, the sum of `K L m`, the tile count and `W.card` is at most `18 * +D^(-1/18) * N²`, where `D = D m` and `N` is the original matrix dimension. The hypotheses +are `D m ^ 18 ≤ N` and `W.card * 2 ^ m ≤ N ^ 2`. This counts finite objects; it does not +count operations on them. -/ +theorem minutiae_total_le (m N : ℕ) (W : Finset (Fin N × Fin N)) (hN : D m ^ 18 ≤ N) + (hW : W.card * 2 ^ m ≤ N ^ 2) : + ((K (19 * m) m + (tiles (19 * m) m N).card + W.card : ℕ) : ℝ) + ≤ 18 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2 := by + obtain ⟨hK, hT, hW'⟩ := minutiae_parts m N W hN hW + have hc := saving_pos m + have hpad : (padN (19 * m) m N : ℝ) ≤ 2 * (N : ℝ) := by + exact_mod_cast padN_le_two_mul_of_D_pow_le m N hN + have hsq : saving m * (padN (19 * m) m N : ℝ) ^ 2 ≤ saving m * (2 * (N : ℝ)) ^ 2 := by gcongr + rw [← saving_eq_D_rpow] + push_cast + -- `K + |tiles| + |W| ≤ c N² + 4 c (2N)² + c N² = 18 c N²` with `c = 2^{-m/9}`. + linarith [hK, hT, hW', hsq] + +/-! ### Bounds on specified integer values -/ + +/-- The subsets in the table are distinct (Section 2.3.4), so a sum over the table that selects the +subset `Q` is one of its terms or 0. -/ +private theorem sum_table_eq_zero_or {L m : ℕ} {A : Type*} [AddCommMonoid A] (lay : Layout L m) + (F : Fin (K0 L m) × Fin (K0 L m) → A) (Q : Finset (Fin L)) : + (∑ gh, if lay.table gh = Q then F gh else 0) = 0 + ∨ ∃ gh, (∑ gh', if lay.table gh' = Q then F gh' else 0) = F gh := by + by_cases hQ : ∃ gh, lay.table gh = Q + · obtain ⟨gh, rfl⟩ := hQ + exact .inr ⟨gh, by simp [lay.table_injective.eq_iff]⟩ + · exact .inl (sum_eq_zero fun gh _ => ite_eq_right (not_exists.mp hQ gh)) + +/-- A nonnegative bound `A` on the absolute values of all entries of `X` also bounds every +value of its left band arrays, for every band index. Zero padding and unused strings are +included; `0 ≤ A` remains required when `N = 0`. -/ +theorem abs_bandArrayL_le {L m N : ℕ} (lay : Layout L m) {X : Matrix (Fin N) (Fin (D m)) ℤ} {A : ℤ} + (hA : 0 ≤ A) (hX : ∀ I k, |X I k| ≤ A) (β : ℕ) (u : LeftStr L) : + |bandArrayL lay X β u| ≤ A := by + unfold bandArrayL arrayL bandFamilyL + split_ifs with h + · rcases sum_table_eq_zero_or lay (fun gh => rowBlock lay X (β * K0 L m + gh.1)) (innerSetL u) + with hzero | ⟨gh, hgh⟩ + · rw [hzero] + simpa using hA + · rw [hgh] + unfold rowBlock padRows + split_ifs + · exact hX _ _ + · simpa using hA + · simpa using hA + +/-- A nonnegative bound `B` on the absolute values of all entries of `Y` also bounds every +value of its right band arrays, for every band index. Zero padding and unused strings are +included; `0 ≤ B` remains required when `N = 0`. -/ +theorem abs_bandArrayR_le {L m N : ℕ} (lay : Layout L m) {Y : Matrix (Fin (D m)) (Fin N) ℤ} {B : ℤ} + (hB : 0 ≤ B) (hY : ∀ k J, |Y k J| ≤ B) (β : ℕ) (v : RightStr L) : + |bandArrayR lay Y β v| ≤ B := by + unfold bandArrayR arrayR bandFamilyR + split_ifs with h + · rcases sum_table_eq_zero_or lay (fun gh => colBlock lay Y (β * K0 L m + gh.2)) (innerSetR v) + with hzero | ⟨gh, hgh⟩ + · rw [hzero] + simpa using hB + · rw [hgh] + unfold colBlock padCols + split_ifs + · exact hY _ _ + · simpa using hB + · simpa using hB + +/-- If every coefficient has absolute value at most one and every array entry has absolute +value at most `A`, then an encoded value has absolute value at most the array index type's +cardinality times `A`. This bounds the value of the finite sum defining `encodeWith`. -/ +theorem abs_encodeWith_le {α : Type} [Fintype α] {c : Term → α → ℤ} (hc : ∀ lam s, |c lam s| ≤ 1) + {L : ℕ} {a : (Fin L → α) → ℤ} {A : ℤ} (ha : ∀ u, |a u| ≤ A) (τ : Leaf L) : + |encodeWith c τ a| ≤ Fintype.card (Fin L → α) * A := by + have h := abs_sum_mul_le univ (c := fun u => ∏ ℓ, c (τ ℓ) (u ℓ)) + (fun u _ => abs_prod_le_one _ fun ℓ _ => hc _ _) fun u _ => ha u + simpa only [encodeWith, mul_comm, card_univ] using h + +/-- An absolute bound `A` on the left array entries gives the bound `7 ^ L * A` on each value +of `encodingL`. -/ +theorem abs_encodingL_le {L : ℕ} {a : LeftStr L → ℤ} {A : ℤ} (ha : ∀ u, |a u| ≤ A) (τ : Leaf L) : + |encodingL a τ| ≤ 7 ^ L * A := by + have h := abs_encodeWith_le abs_phi_le_one ha τ + rw [card_leftStr] at h + exact_mod_cast h + +/-- An absolute bound `B` on the right array entries gives the bound `7 ^ L * B` on each value +of `encodingR`. -/ +theorem abs_encodingR_le {L : ℕ} {b : RightStr L → ℤ} {B : ℤ} (hb : ∀ v, |b v| ≤ B) (τ : Leaf L) : + |encodingR b τ| ≤ 7 ^ L * B := by + have h := abs_encodeWith_le abs_psi_le_one hb τ + rw [card_rightStr] at h + exact_mod_cast h + +/-- If `encA` and `encB` are bounded in absolute value by `A` and `B`, every value returned by +`Pruned.result` is bounded by `10 ^ n * (A * B)`. This holds for every request set `S` and +output string `w`, including unrequested strings, whose returned value is zero. It bounds +the specified mathematical result, not native scratch states or machine word width. -/ +theorem Pruned.abs_result_le {n : ℕ} (encA encB : Leaf n → ℤ) (A B : ℤ) (hA : ∀ τ, |encA τ| ≤ A) + (hB : ∀ τ, |encB τ| ≤ B) (S : Finset (OutStr n)) (w : OutStr n) : + |Pruned.result n encA encB S w| ≤ 10 ^ n * (A * B) := by + have hA0 : 0 ≤ A := (abs_nonneg _).trans (hA fun _ => Term.P0) + have hAB : 0 ≤ A * B := mul_nonneg hA0 ((abs_nonneg _).trans (hB fun _ => Term.P0)) + -- The value is 0, or the sum over some of the `10^n` leaves of the product of the two numbers + -- looked up there. + rw [Pruned.result_eq_piecewise_sum] + simp only [Finset.piecewise] + split_ifs + · calc |∑ τ : Leaf n with Leaf.Contributes τ w, encA τ * encB τ| + ≤ (univ.filter fun τ : Leaf n => Leaf.Contributes τ w).card * (A * B) := + abs_sum_le_card_mul _ fun τ _ => by + rw [abs_mul] + exact mul_le_mul (hA τ) (hB τ) (abs_nonneg _) hA0 + _ ≤ 10 ^ n * (A * B) := by + gcongr + exact_mod_cast (card_le_univ _).trans_eq (card_vertex n) + · simp only [Pi.zero_apply, abs_zero] + positivity + +/-- Under `D m ^ 18 ≤ N`, the number `10 ^ (19 * m)` of leaf indices is at most `N ^ 2`. -/ +theorem ten_pow_le_sq (m N : ℕ) (hN : D m ^ 18 ≤ N) : 10 ^ (19 * m) ≤ N ^ 2 := by + -- `N₀ ≤ K N₀ ≤ N` and `2^{-m/9} ≤ 1`, so equation (6) gives `10^L ≤ N N₀ ≤ N²`. + have hN0 : (N0 (19 * m) m : ℝ) ≤ N := by + exact_mod_cast (Nat.le_mul_of_pos_left _ (Nat.choose_pos (by omega))).trans + (K_mul_N0_le m N hN) + have hc := saving_le_one m + have hreal : (10 : ℝ) ^ (19 * m) ≤ (N : ℝ) ^ 2 := + calc (10 : ℝ) ^ (19 * m) ≤ saving m * N * N0 (19 * m) m := encoding_size_le m N hN + _ ≤ 1 * N * N := by gcongr + _ = (N : ℝ) ^ 2 := by ring + exact_mod_cast hreal + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Orders.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Orders.lean new file mode 100644 index 0000000000..facc0f50c2 --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Orders.lean @@ -0,0 +1,311 @@ +/- +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.MatrixMultiplication.Schonhage.Levels +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Pruned +public import Complexity.Foundations.Combinatorics.FunctionPlaces +public import Mathlib.Algebra.Order.Ring.Rat +import Mathlib.Tactic.Linarith +import Mathlib.Tactic.Ring +import Mathlib.Tactic.Positivity + +/-! +# Orders and counts of contributing leaves + +## Main definitions + +* `order`, `alpha`, `beta`: leaf order and its two counting distributions. + +## Main results + +* `card_contributing_of_order`, `card_filter_order_eq`, `card_outStr_of_leaf`: exact counts. +* `two_pow_mul_beta_le`: geometric decay when `L = 19m`. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The set of levels at which a leaf chooses `P₀`. -/ +def P0Levels {L : ℕ} (τ : Leaf L) : Finset (Fin L) := univ.filter fun ℓ => τ ℓ = Term.P0 + +/-- The integer difference between `m` and the number of `P0` positions of a leaf. It can be +negative for arbitrary leaves; nonnegativity is proved for leaves contributing to an +output with `m` inner positions. -/ +def order {L : ℕ} (m : ℕ) (τ : Leaf L) : ℤ := (m : ℤ) - ((P0Levels τ).card : ℤ) + +/-- The quantity `m.choose d * 9 ^ d`, counting order-`d` leaves contributing to one output +with `m` inner positions. -/ +def alpha (m d : ℕ) : ℕ := m.choose d * 9 ^ d + +/-- The quantity `L.choose (m - d) * 9 ^ (L - m + d)`, using natural-number subtraction. Its +interpretation as the number of order-`d` leaves requires `d ≤ m ≤ L`. -/ +def beta (L m d : ℕ) : ℕ := L.choose (m - d) * 9 ^ (L - m + d) + +/-! ### One level: which terms contribute to which variables -/ + +/-- The private term of `z` is `P₀` exactly if `z = z₀`. -/ +theorem OutVar.privateTerm_eq_P0_iff (z : OutVar) : z.privateTerm = Term.P0 ↔ z.IsInner := by + cases z <;> simp [OutVar.IsInner, OutVar.privateTerm] + +/-- One term is `P₀`. -/ +private theorem card_term_P0 : (univ.filter fun lam : Term => lam = Term.P0).card = 1 := by decide + +/-- Nine terms are different from `P₀`. -/ +private theorem card_term_ne_P0 : (univ.filter fun lam : Term => ¬ lam = Term.P0).card = 9 := by + decide + +/-- Among the terms contributing to `z`, the number of those equal to `P₀`: one if `z = z₀`, none if +`z = z_ij`. -/ +private theorem card_contributing_P0 (z : OutVar) : + ((univ.filter fun lam : Term => lam.Contributes z).filter fun lam => lam = Term.P0).card + = if z.IsInner then 1 else 0 := by + decide +revert + +/-- Among the terms contributing to `z`, the number of those different from `P₀`: nine if `z = z₀`, +one if `z = z_ij`. -/ +private theorem card_contributing_ne_P0 (z : OutVar) : + ((univ.filter fun lam : Term => lam.Contributes z).filter fun lam => ¬ lam = Term.P0).card + = if z.IsInner then 9 else 1 := by + decide +revert + +/-- The number of terms contributing to `z`: ten if `z = z₀`, one if `z = z_ij`. -/ +private theorem card_contributing_term (z : OutVar) : + (univ.filter fun lam : Term => lam.Contributes z).card = if z.IsInner then 10 else 1 := by + decide +revert + +/-- Among the output variables to which `lam` contributes, exactly one is inner. -/ +private theorem card_contributed_inner (lam : Term) : + ((univ.filter fun z : OutVar => lam.Contributes z).filter fun z => z.IsInner).card = 1 := by + decide +revert + +/-- Among the output variables to which `lam` contributes, the number of outer ones: none if +`lam = P₀`, one if `lam = P_ij`. -/ +private theorem card_contributed_outer (lam : Term) : + ((univ.filter fun z : OutVar => lam.Contributes z).filter fun z => ¬ z.IsInner).card + = if lam = Term.P0 then 0 else 1 := by + decide +revert + +/-! ### The order of a leaf -/ + +/-- Membership in the set of levels at which a leaf chooses `P₀`. -/ +@[simp] +theorem mem_P0Levels {L : ℕ} {τ : Leaf L} {ℓ : Fin L} : ℓ ∈ P0Levels τ ↔ τ ℓ = Term.P0 := by + simp [P0Levels] + +/-- Membership in the inner set of an output string. -/ +@[simp] +theorem mem_innerSetO {L : ℕ} {η : OutStr L} {ℓ : Fin L} : ℓ ∈ innerSetO η ↔ (η ℓ).IsInner := by + simp [innerSetO] + +/-- A contributing leaf can use `P0` only at an inner position of the output string. -/ +theorem P0Levels_subset {L : ℕ} {τ : Leaf L} {w : OutStr L} (h : Leaf.Contributes τ w) : + P0Levels τ ⊆ innerSetO w := by + intro ℓ hℓ + rw [mem_P0Levels] at hℓ + rw [mem_innerSetO] + rcases (Term.contributes_iff _ _).mp (h ℓ) with hinner | hprivate + · exact hinner + · exact (OutVar.privateTerm_eq_P0_iff _).mp (hprivate.symm.trans hℓ) + +/-- A leaf contributing to a string with exactly `m` inner positions has nonnegative order +relative to `m`. -/ +theorem order_nonneg {L m : ℕ} {w : OutStr L} (hw : (innerSetO w).card = m) {τ : Leaf L} + (h : Leaf.Contributes τ w) : 0 ≤ order m τ := by + have hc := card_le_card (P0Levels_subset h) + unfold order + omega + +/-- The order relative to `m` of any leaf is at most `m`, regardless of which outputs it +contributes to. -/ +theorem order_le {L : ℕ} (m : ℕ) (τ : Leaf L) : order m τ ≤ m := by + unfold order + omega + +/-- A leaf has order `d` exactly if it chooses `P₀` at `m - d` levels. -/ +theorem order_eq_iff {L : ℕ} (m d : ℕ) (τ : Leaf L) : + order m τ = d ↔ (P0Levels τ).card + d = m := by + unfold order + omega + +/-! ### The counts -/ + +/-- Exactly `10 ^ m` leaves contribute to an output string with `m` inner positions. -/ +theorem card_filter_contributes {L m : ℕ} (w : OutStr L) (hw : (innerSetO w).card = m) : + (univ.filter fun τ : Leaf L => Leaf.Contributes τ w).card = 10 ^ m := by + -- Ten terms contribute to `z₀` and one to `z_ij`. + simp only [Leaf.filter_contributes_eq_piFinset, Fintype.card_piFinset, + card_contributing_term] + rw [prod_ite, prod_const, prod_const_one, mul_one, ← hw] + rfl + +/-- For a string with `m` inner positions, `alpha m d` leaves of order `d` contribute to it. +The natural parameter `d` is unrestricted; for `d > m` both sides are zero. -/ +theorem card_contributing_of_order {L m : ℕ} (w : OutStr L) (hw : (innerSetO w).card = m) + (d : ℕ) : + (univ.filter fun τ : Leaf L => Leaf.Contributes τ w ∧ order m τ = d).card = alpha m d := by + by_cases hd : d ≤ m + · -- Sort the strings of terms that contribute to `w` level by level and have `P₀` at exactly + -- `m - d` levels by the set `S` of these levels. + have hcount := Finset.card_pi_places_card + (fun ℓ => univ.filter fun lam : Term => lam.Contributes (w ℓ)) (fun lam => lam = Term.P0) + (m - d) + simp only [card_contributing_P0, card_contributing_ne_P0] at hcount + -- Given `S`, there is no such string unless `S ⊆ Q`. Then there are nine choices at each of + -- the `|Q| - |S| = d` levels of `Q` outside `S`, and one choice at every other level. + have hgiven : ∀ S ∈ powersetCard (m - d) (univ : Finset (Fin L)), + (∏ ℓ, if ℓ ∈ S then (if (w ℓ).IsInner then 1 else 0) else (if (w ℓ).IsInner then 9 else 1)) + = if S ⊆ innerSetO w then 9 ^ d else 0 := by + intro S hS + have hQS := prod_ite_ite_subset (innerSetO w) S 9 + rw [hw, (mem_powersetCard.mp hS).2, Nat.sub_sub_self hd] at hQS + rw [← hQS] + exact prod_congr rfl fun ℓ _ => by simp + -- There are `binom(m, m - d) = binom(m, d)` sets `S ⊆ Q`. + rw [sum_congr rfl hgiven, sum_powersetCard_ite_subset, hw, Nat.choose_symm hd] at hcount + -- These strings are the leaves of order `d` that contribute to `w`. + rw [alpha, ← hcount] + refine congrArg card (Finset.ext fun τ => ?_) + simp only [mem_filter, mem_univ, true_and, order_eq_iff] + exact and_congr Iff.rfl (by unfold P0Levels; omega) + · -- No leaf has order more than `m`. + rw [alpha, Nat.choose_eq_zero_of_lt (by omega), zero_mul, card_eq_zero] + refine filter_eq_empty_iff.mpr fun τ _ h => ?_ + have hle := order_le m τ + omega + +/-- If `d ≤ m ≤ L`, the full leaf index type contains exactly `beta L m d` leaves of order +`d`. -/ +theorem card_filter_order_eq {L m : ℕ} (hmL : m ≤ L) (d : ℕ) (hd : d ≤ m) : + (univ.filter fun τ : Leaf L => order m τ = d).card = beta L m d := by + -- The leaves of order `d` are the strings of terms with `P₀` at exactly `m - d` levels. + have hsort : (univ.filter fun τ : Leaf L => order m τ = d) + = univ.filter fun τ : Leaf L => (univ.filter fun ℓ => τ ℓ = Term.P0).card = m - d := by + ext τ + simp only [mem_filter, mem_univ, true_and, order_eq_iff] + unfold P0Levels + omega + rw [hsort, Finset.card_places_card (fun lam : Term => lam = Term.P0) (m - d), card_term_P0, + card_term_ne_P0, one_pow, one_mul, Fintype.card_fin, beta, + show L - (m - d) = L - m + d by omega] + +/-- Evaluate `beta` at order zero. -/ +theorem beta_zero (L m : ℕ) : beta L m 0 = L.choose m * 9 ^ (L - m) := by + simp [beta] + +/-- The order-zero count formula equals `M L m`. -/ +theorem beta_zero_eq_M (L m : ℕ) : beta L m 0 = M L m := by + rw [beta_zero, M, K, N0, ← pow_mul, mul_comm (L - m) 2, pow_mul] + norm_num + +/-- Assuming `m ≤ L`, a leaf of natural order `d` contributes to exactly `(L - m + d).choose +d` output strings with `m` inner positions. The order hypothesis forces `d ≤ m`. -/ +theorem card_outStr_of_leaf {L m : ℕ} (hmL : m ≤ L) (d : ℕ) (τ : Leaf L) + (hτ : order m τ = d) : + (univ.filter fun w : OutStr L => (innerSetO w).card = m ∧ Leaf.Contributes τ w).card + = (L - m + d).choose d := by + have hP0 : (P0Levels τ).card + d = m := (order_eq_iff m d τ).mp hτ + -- Sort the output strings to which `τ` contributes level by level, with `z₀` at exactly `m` + -- levels, by the set `S` of these levels, which is their inner set. + have hcount := Finset.card_pi_places_card + (fun ℓ => univ.filter fun z : OutVar => (τ ℓ).Contributes z) (fun z => z.IsInner) m + simp only [card_contributed_inner, card_contributed_outer] at hcount + -- Given `S`, there is one such string if `S` contains the levels at which `τ` chooses `P₀`, and + -- none otherwise. + have hgiven : ∀ S ∈ powersetCard m (univ : Finset (Fin L)), + (∏ ℓ, if ℓ ∈ S then 1 else (if τ ℓ = Term.P0 then 0 else 1)) + = if P0Levels τ ⊆ S then 1 else 0 := by + intro S _ + rw [← prod_ite_ite_superset] + exact prod_congr rfl fun ℓ _ => by simp + -- A set of `m` levels that contains these `m - d` levels is given by the `L - m` levels outside + -- it, among the other `L - m + d` levels. + rw [sum_congr rfl hgiven, ← card_filter, + card_powersetCard_superset _ (by simpa using hmL), Fintype.card_fin, + show L - (P0Levels τ).card = L - m + d by omega, Nat.choose_symm_add] at hcount + rw [← hcount] + refine congrArg card (Finset.ext fun w => ?_) + simp only [mem_filter, mem_univ, true_and] + exact and_comm + +/-- `β_d` is positive when `m ≤ L`. -/ +theorem beta_pos {L m : ℕ} (hmL : m ≤ L) (d : ℕ) : 0 < beta L m d := + Nat.mul_pos (Nat.choose_pos (by omega)) (by positivity) + +/-- For `1 ≤ d ≤ m ≤ L`, consecutive values of `beta` satisfy the displayed multiplication +identity, without division (cf. [AlmanWilliams2026, equation (5)]). -/ +theorem Equation5.ratio_nat {L m : ℕ} (hmL : m ≤ L) (d : ℕ) (hd1 : 1 ≤ d) (hd : d ≤ m) : + beta L m d * (L - m + d) = 9 * (m - d + 1) * beta L m (d - 1) := by + -- With `k = m - d` this is `binom(L, k) (L - k) = binom(L, k + 1) (k + 1)`, times a power of 9. + have hchoose : L.choose (m - d + 1) * (m - d + 1) = L.choose (m - d) * (L - m + d) := by + rw [Nat.choose_succ_right_eq, show L - (m - d) = L - m + d by omega] + have hpow : 9 ^ (L - m + d) = 9 * 9 ^ (L - m + (d - 1)) := by + rw [show L - m + d = L - m + (d - 1) + 1 by omega, pow_succ, mul_comm] + rw [beta, beta, show m - (d - 1) = m - d + 1 by omega, hpow] + calc L.choose (m - d) * (9 * 9 ^ (L - m + (d - 1))) * (L - m + d) + = 9 * 9 ^ (L - m + (d - 1)) * (L.choose (m - d) * (L - m + d)) := by ring + _ = 9 * 9 ^ (L - m + (d - 1)) * (L.choose (m - d + 1) * (m - d + 1)) := by rw [hchoose] + _ = 9 * (m - d + 1) * (L.choose (m - d + 1) * 9 ^ (L - m + (d - 1))) := by ring + +/-- For `1 ≤ d ≤ m ≤ L`, the ratio of consecutive `beta` values is `9 * (m - d + 1) / (L - m + +d)`, with the ratio and subtractions interpreted in `ℚ`. This extends the parameter ranges +used in [AlmanWilliams2026, equations (5) and (7)]. -/ +theorem beta_ratio {L m : ℕ} (hmL : m ≤ L) (d : ℕ) (hd1 : 1 ≤ d) (hd : d ≤ m) : + (beta L m d : ℚ) / (beta L m (d - 1) : ℚ) + = 9 * ((m : ℚ) - (d : ℚ) + 1) / ((L : ℚ) - (m : ℚ) + (d : ℚ)) := by + have hm : (m : ℚ) ≤ L := by exact_mod_cast hmL + have hd0 : (1 : ℚ) ≤ d := by exact_mod_cast hd1 + have hbeta : (0 : ℚ) < beta L m (d - 1) := by exact_mod_cast beta_pos hmL (d - 1) + have hnat := congrArg (Nat.cast : ℕ → ℚ) (Equation5.ratio_nat hmL d hd1 hd) + push_cast [Nat.cast_sub hmL, Nat.cast_sub hd] at hnat + rw [div_eq_div_iff hbeta.ne' (by linarith), hnat] + +/-- At `L = 19 * m` and `1 ≤ d ≤ m`, the ratio of consecutive `beta` values is bounded above +by `9 * m / (18 * m + 1)`, which is strictly below `1 / 2`. All divisions here are in `ℚ`. -/ +theorem beta_ratio_lt_half (m d : ℕ) (hd1 : 1 ≤ d) (hd : d ≤ m) : + (beta (19 * m) m d : ℚ) / (beta (19 * m) m (d - 1) : ℚ) ≤ 9 * (m : ℚ) / (18 * (m : ℚ) + 1) + ∧ 9 * (m : ℚ) / (18 * (m : ℚ) + 1) < 1 / 2 := by + have hm : (1 : ℚ) ≤ m := by exact_mod_cast hd1.trans hd + have hd0 : (1 : ℚ) ≤ d := by exact_mod_cast hd1 + refine ⟨?_, ?_⟩ + · -- The numerator is at most `9m` and the denominator at least `18m + 1`, because `d ≥ 1`. + rw [beta_ratio (by omega) d hd1 hd] + push_cast + exact div_le_div₀ (by positivity) (by linarith) (by positivity) (by linarith) + · rw [div_lt_div_iff₀ (by positivity) (by norm_num)] + linarith + +/-- For `d ≤ m` and `L = 19 * m`, the order count satisfies `2 ^ d * beta L m d ≤ M L m`. This +is the division-free geometric bound from [AlmanWilliams2026, equation (5)]. -/ +theorem two_pow_mul_beta_le (m d : ℕ) (hd : d ≤ m) : 2 ^ d * beta (19 * m) m d ≤ M (19 * m) m := by + induction d with + | zero => simp [beta_zero_eq_M] + | succ d ih => + -- By the two inequalities of (5), `β_{d+1} < β_d / 2`. + obtain ⟨hle, hlt⟩ := beta_ratio_lt_half m (d + 1) (by omega) hd + rw [Nat.add_sub_cancel] at hle + have hbeta : (0 : ℚ) < beta (19 * m) m d := by exact_mod_cast beta_pos (by omega) d + have hhalf : 2 * beta (19 * m) m (d + 1) ≤ beta (19 * m) m d := by + have hq := (div_lt_iff₀ hbeta).mp (hle.trans_lt hlt) + exact_mod_cast (by linarith : (2 * beta (19 * m) m (d + 1) : ℚ) ≤ beta (19 * m) m d) + calc 2 ^ (d + 1) * beta (19 * m) m (d + 1) = 2 ^ d * (2 * beta (19 * m) m (d + 1)) := by ring + _ ≤ 2 ^ d * beta (19 * m) m d := Nat.mul_le_mul_left _ hhalf + _ ≤ M (19 * m) m := ih (by omega) + + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Parameters.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Parameters.lean new file mode 100644 index 0000000000..38a936be5b --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Parameters.lean @@ -0,0 +1,131 @@ +/- +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.MatrixMultiplication.Schonhage.Sparse +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Tiling +public import Complexity.Foundations.Data.Nat.Choose.Real +public import Mathlib.Analysis.Complex.ExponentialBounds + +/-! +# Parameters of finite thin multiplication + +## Main results + +* `K_mul_N0_le_D_pow`: bound `K L m * N0 L m` by `D m ^ 18` at `L = 19 * m`. +* `encoding_size_le`: the shared encoding length bound. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- The binomial count at `L = 19 * m` is at most `(19 * Real.exp 1) ^ m`. This uses +`Nat.choose_le_exp_mul_div_pow` and includes `m = 0`, when both sides are one. -/ +theorem K_le_pow (m : ℕ) : (K (19 * m) m : ℝ) ≤ (19 * Real.exp 1) ^ m := by + rcases Nat.eq_zero_or_pos m with rfl | hm + · simp [K] + · have hm' : (m : ℝ) ≠ 0 := by exact_mod_cast hm.ne' + have hbase : Real.exp 1 * ((19 * m : ℕ) : ℝ) / (m : ℝ) = 19 * Real.exp 1 := by + push_cast + field_simp + exact hbase ▸ Nat.choose_le_exp_mul_div_pow (19 * m) m + +/-- The real constant `19 * exp 1 * 3 ^ 18` is strictly less than `4 ^ 18`. -/ +theorem nineteen_mul_exp_mul_three_pow_lt : 19 * Real.exp 1 * (3 : ℝ) ^ 18 < (4 : ℝ) ^ 18 := by + calc 19 * Real.exp 1 * (3 : ℝ) ^ 18 < 19 * 2.7182818286 * (3 : ℝ) ^ 18 := by + gcongr + exact Real.exp_one_lt_d9 + _ < (4 : ℝ) ^ 18 := by norm_num + +/-- `4^{18m} = D^18`. -/ +private theorem four_pow_eq (m : ℕ) : ((4 : ℝ) ^ 18) ^ m = (D m : ℝ) ^ 18 := by + rw [D, Nat.cast_pow, Nat.cast_ofNat, ← pow_mul, ← pow_mul, mul_comm] + +/-- At `L = 19 * m`, the product `K L m * N0 L m` is at most `D m ^ 18`. -/ +theorem K_mul_N0_le_D_pow (m : ℕ) : K (19 * m) m * N0 (19 * m) m ≤ D m ^ 18 := by + have hK := K_le_pow m + have key : (K (19 * m) m : ℝ) * (N0 (19 * m) m : ℝ) ≤ (D m : ℝ) ^ 18 := by + calc (K (19 * m) m : ℝ) * (N0 (19 * m) m : ℝ) + ≤ (19 * Real.exp 1) ^ m * ((3 : ℝ) ^ 18) ^ m := by + rw [N0_nineteen_mul m] + push_cast + rw [pow_mul] + exact mul_le_mul_of_nonneg_right hK (by positivity) + _ = (19 * Real.exp 1 * (3 : ℝ) ^ 18) ^ m := (mul_pow _ _ _).symm + _ ≤ ((4 : ℝ) ^ 18) ^ m := + pow_le_pow_left₀ (by positivity) nineteen_mul_exp_mul_three_pow_lt.le m + _ = (D m : ℝ) ^ 18 := four_pow_eq m + exact_mod_cast key + +/-- The size condition `D m ^ 18 ≤ N` bounds `K (19 * m) m * N0 (19 * m) m` by `N`. -/ +theorem K_mul_N0_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : K (19 * m) m * N0 (19 * m) m ≤ N := + (K_mul_N0_le_D_pow m).trans hN + +/-- The size condition `D m ^ 18 ≤ N` bounds the band dimension at `L = 19 * m` by `N`. -/ +theorem K0_mul_N0_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : K0 (19 * m) m * N0 (19 * m) m ≤ N := + (Nat.mul_le_mul_right _ (Nat.sqrt_le_self _)).trans (K_mul_N0_le m N hN) + +/-- Under `D m ^ 18 ≤ N`, padding to complete bands at `L = 19 * m` gives a dimension at most +`2 * N`. -/ +theorem padN_le_two_mul_of_D_pow_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : padN (19 * m) m N ≤ 2 * N := + padN_le_two_mul (19 * m) m N (K0_mul_N0_le m N hN) + +/-- The real constant `10 ^ 19 * 2 ^ (1 / 9) / 3 ^ 18` is strictly less than `4 ^ 18`. -/ +theorem ten_pow_mul_two_rpow_div_three_pow_lt : + (10 : ℝ) ^ 19 * (2 : ℝ) ^ (1 / 9 : ℝ) / (3 : ℝ) ^ 18 < (4 : ℝ) ^ 18 := by + have htwo : (2 : ℝ) ^ (1 / 9 : ℝ) ≤ 2 := by + calc (2 : ℝ) ^ (1 / 9 : ℝ) ≤ (2 : ℝ) ^ (1 : ℝ) := + Real.rpow_le_rpow_of_exponent_le (by norm_num) (by norm_num) + _ = 2 := Real.rpow_one 2 + rw [div_lt_iff₀ (by positivity)] + calc (10 : ℝ) ^ 19 * (2 : ℝ) ^ (1 / 9 : ℝ) ≤ (10 : ℝ) ^ 19 * 2 := + mul_le_mul_of_nonneg_left htwo (by positivity) + _ < (4 : ℝ) ^ 18 * (3 : ℝ) ^ 18 := by norm_num + +/-- At `L = 19 * m`, the real quotient `10 ^ L * 2 ^ (m / 9) / N0 L m` is at most `D m ^ 18`. -/ +theorem ten_pow_div_N0_le (m : ℕ) : + (10 : ℝ) ^ (19 * m) * (2 : ℝ) ^ ((m : ℝ) / 9) / (N0 (19 * m) m : ℝ) ≤ (D m : ℝ) ^ 18 := by + have hroot : (2 : ℝ) ^ ((m : ℝ) / 9) = ((2 : ℝ) ^ (1 / 9 : ℝ)) ^ m := by + rw [← Real.rpow_natCast _ m, ← Real.rpow_mul (by norm_num)] + congr 1 + ring + calc (10 : ℝ) ^ (19 * m) * (2 : ℝ) ^ ((m : ℝ) / 9) / (N0 (19 * m) m : ℝ) + = ((10 : ℝ) ^ 19 * (2 : ℝ) ^ (1 / 9 : ℝ) / (3 : ℝ) ^ 18) ^ m := by + rw [N0_nineteen_mul m, hroot] + push_cast + rw [pow_mul, pow_mul, div_pow, mul_pow] + _ ≤ ((4 : ℝ) ^ 18) ^ m := pow_le_pow_left₀ (by positivity) + ten_pow_mul_two_rpow_div_three_pow_lt.le m + _ = (D m : ℝ) ^ 18 := four_pow_eq m + +/-- For `L = 19 * m` and `D m ^ 18 ≤ N`, the encoding index count `10 ^ L` is at most `2 ^ (-m +/ 9) * N * N0 L m`, as a real inequality (cf. [AlmanWilliams2026, equation (6)]). -/ +theorem encoding_size_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : + (10 : ℝ) ^ (19 * m) ≤ (2 : ℝ) ^ (-(m : ℝ) / 9) * (N : ℝ) * (N0 (19 * m) m : ℝ) := by + have hN0 : (0 : ℝ) < (N0 (19 * m) m : ℝ) := by exact_mod_cast N0_pos _ _ + have hN' : (D m : ℝ) ^ 18 ≤ (N : ℝ) := by exact_mod_cast hN + have h := (ten_pow_div_N0_le m).trans hN' + rw [div_le_iff₀ hN0] at h + calc (10 : ℝ) ^ (19 * m) + = (10 : ℝ) ^ (19 * m) * (2 : ℝ) ^ ((m : ℝ) / 9) * saving m := by + rw [mul_assoc, mul_comm _ (saving m), saving_mul_two_rpow, mul_one] + _ ≤ (N : ℝ) * (N0 (19 * m) m : ℝ) * saving m := + mul_le_mul_of_nonneg_right h (saving_pos m).le + _ = saving m * (N : ℝ) * (N0 (19 * m) m : ℝ) := by ring + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Sparse.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Sparse.lean new file mode 100644 index 0000000000..1b2cab585a --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Sparse.lean @@ -0,0 +1,309 @@ +/- +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.MatrixMultiplication.Schonhage.Orders +public import Complexity.Foundations.Algebra.GeometricTail +public import Mathlib.Analysis.SpecialFunctions.Pow.Real + +/-! +# Sparse query bounds + +## Main definitions + +* `saving`: the factor `D^(-1/18)` for `D = 4^m`. + +## Main results + +* `card_leaves_le_sum_min`: sum the two order-wise bounds. +* `card_leaves_le_saving`: explicit sparse leaf bound, with no assumed sparse solver. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-! ### Specialization to `L = 19 * m` -/ + +/-- At `L = 19 * m`, the outer dimension is `3 ^ (18 * m)`. -/ +theorem N0_nineteen_mul (m : ℕ) : N0 (19 * m) m = 3 ^ (18 * m) := by + unfold N0 + congr 1 + omega + +/-- `M` is positive when `m ≤ L`. -/ +theorem M_pos {L m : ℕ} (hmL : m ≤ L) : 0 < M L m := + Nat.mul_pos (Nat.choose_pos hmL) (Nat.pow_pos (Nat.pow_pos (by decide))) + +/-- `D = 4^m = 2^{2m}`. -/ +private theorem D_eq_two_pow (m : ℕ) : (D m : ℝ) = (2 : ℝ) ^ (2 * m) := by + rw [D, pow_mul] + norm_num + +/-- The real square root of `D m` is `2 ^ m`. -/ +theorem two_pow_eq_sqrt_D (m : ℕ) : (2 : ℝ) ^ m = Real.sqrt (D m) := by + rw [D_eq_two_pow, pow_mul', Real.sqrt_sq (by positivity)] + +/-- The positive real factor `2 ^ (-m / 9)`, used in the finite leaf and encoding count +bounds. -/ +noncomputable def saving (m : ℕ) : ℝ := (2 : ℝ) ^ (-(m : ℝ) / 9) + +/-- `2^{-m/9}` is positive. -/ +theorem saving_pos (m : ℕ) : 0 < saving m := + Real.rpow_pos_of_pos two_pos _ + +/-- `2^{-m/9} ≤ 1`. -/ +theorem saving_le_one (m : ℕ) : saving m ≤ 1 := + Real.rpow_le_one_of_one_le_of_nonpos one_le_two + (div_nonpos_of_nonpos_of_nonneg (neg_nonpos.mpr m.cast_nonneg) (by norm_num)) + +/-- `2^{-m/9} · 2^{m/9} = 1`. -/ +theorem saving_mul_two_rpow (m : ℕ) : saving m * (2 : ℝ) ^ ((m : ℝ) / 9) = 1 := by + rw [saving, ← Real.rpow_add two_pos, show -(m : ℝ) / 9 + (m : ℝ) / 9 = 0 by ring, Real.rpow_zero] + +/-- Express `saving m` as the real power `D^(-1/18)` for `D = D m`. -/ +theorem saving_eq_D_rpow (m : ℕ) : saving m = (D m : ℝ) ^ (-(1 / 18 : ℝ)) := by + rw [saving, D_eq_two_pow, ← Real.rpow_natCast 2 (2 * m), ← Real.rpow_mul two_pos.le] + congr 1 + push_cast + ring + +/-! ### The first inequality -/ + +/-- If every string in `U` has `m` inner positions, its contributing leaves of natural order +`d` are bounded both by `U.card * alpha m d` and by `beta L m d`. The minimum combines +these bounds without counting a shared leaf more than once. -/ +theorem card_leaves_filter_order_le {L m : ℕ} (U : Finset (OutStr L)) + (hU : ∀ w ∈ U, (innerSetO w).card = m) (d : ℕ) : + ((leaves U).filter fun τ => order m τ = d).card ≤ min (U.card * alpha m d) (beta L m d) := by + refine le_min ?_ ?_ + · -- Charge each output string of U for its leaves of order d. + calc _ ≤ (U.biUnion fun w => + univ.filter fun τ : Leaf L => Leaf.Contributes τ w ∧ order m τ = d).card := by + refine card_le_card fun τ hτ => ?_ + obtain ⟨hτU, hτd⟩ := mem_filter.mp hτ + obtain ⟨w, hw, hc⟩ := mem_leaves.mp hτU + exact mem_biUnion.mpr ⟨w, hw, mem_filter.mpr ⟨mem_univ _, hc, hτd⟩⟩ + _ ≤ U.card * alpha m d := + card_biUnion_le_card_mul _ _ _ fun w hw => + (card_contributing_of_order w (hU w hw) d).le + · -- There are only β_d leaves of order d in all. + rcases U.eq_empty_or_nonempty with rfl | ⟨w₀, hw₀⟩ + · simp [leaves] + have hmL : m ≤ L := by + have h := card_le_univ (innerSetO w₀) + rwa [hU w₀ hw₀, Fintype.card_fin] at h + rcases le_or_gt d m with hd | hd + · rw [← card_filter_order_eq hmL d hd] + exact card_le_card (filter_subset_filter _ (subset_univ _)) + · -- No leaf has an order above m. + have hempty : ((leaves U).filter fun τ => order m τ = d) = ∅ := by + refine filter_false_of_mem fun τ _ h => ?_ + have hle := order_le m τ + omega + rw [hempty] + exact Nat.zero_le _ + +/-- Bound the cardinality of `leaves U` by summing the smaller of the per-output bound and the +full leaf count over natural orders `0` through `m`. Every string in `U` must have exactly +`m` inner positions; the empty set is permitted (cf. [AlmanWilliams2026, Lemma 11]). -/ +theorem card_leaves_le_sum_min {L m : ℕ} (U : Finset (OutStr L)) + (hU : ∀ w ∈ U, (innerSetO w).card = m) : + (leaves U).card ≤ ∑ d ∈ range (m + 1), min (U.card * alpha m d) (beta L m d) := by + -- Every leaf of `leaves(U)` has an order between 0 and m. + have hord : ∀ τ ∈ leaves U, 0 ≤ order m τ := by + intro τ hτ + obtain ⟨w, hw, hc⟩ := mem_leaves.mp hτ + exact order_nonneg (hU w hw) hc + have hmaps : ∀ τ ∈ leaves U, (order m τ).toNat ∈ range (m + 1) := by + intro τ _ + have hle := order_le m τ + rw [mem_range] + omega + -- Sum the smaller of the two bounds over the orders. + rw [card_eq_sum_card_fiberwise hmaps] + refine sum_le_sum fun d _ => le_trans ?_ (card_leaves_filter_order_le U hU d) + refine card_le_card fun τ hτ => ?_ + obtain ⟨hτU, hτd⟩ := mem_filter.mp hτ + have hnonneg := hord τ hτU + exact mem_filter.mpr ⟨hτU, by omega⟩ + +/-! ### The second inequality -/ + +/-- Equation (5) over the reals: `β_d ≤ M 2^{-d}` for `0 ≤ d ≤ m` and `L = 19m`. -/ +private theorem beta_le_M_mul (m d : ℕ) (hd : d ≤ m) : + (beta (19 * m) m d : ℝ) ≤ (M (19 * m) m : ℝ) * (1 / 2 : ℝ) ^ d := by + have h : (2 : ℝ) ^ d * (beta (19 * m) m d : ℝ) ≤ (M (19 * m) m : ℝ) := by + exact_mod_cast two_pow_mul_beta_le m d hd + rw [one_div, inv_pow, ← div_eq_mul_inv, le_div_iff₀ (by positivity)] + linarith + +/-- By (5), the orders `d ≥ a` contribute `∑ β_d ≤ M ∑_{d ≥ a} 2^{-d} < 2 · 2^{-a} M`. -/ +private theorem sum_beta_lt_of_subset (m a : ℕ) (s : Finset ℕ) (hs : s ⊆ Ico a (m + 1)) : + ((∑ d ∈ s, beta (19 * m) m d : ℕ) : ℝ) < 2 * (1 / 2 : ℝ) ^ a * (M (19 * m) m : ℝ) := by + have hM : (0 : ℝ) < (M (19 * m) m : ℝ) := by exact_mod_cast M_pos (show m ≤ 19 * m by omega) + push_cast + calc ∑ d ∈ s, (beta (19 * m) m d : ℝ) + ≤ ∑ d ∈ s, (M (19 * m) m : ℝ) * (1 / 2 : ℝ) ^ d := + sum_le_sum fun d hd => + beta_le_M_mul m d (by have := mem_Ico.mp (hs hd); omega) + _ ≤ ∑ d ∈ Ico a (m + 1), (M (19 * m) m : ℝ) * (1 / 2 : ℝ) ^ d := + sum_le_sum_of_subset_of_nonneg hs fun _ _ _ => by positivity + _ = (M (19 * m) m : ℝ) * ∑ d ∈ Ico a (m + 1), (1 / 2 : ℝ) ^ d := (mul_sum _ _ _).symm + _ < (M (19 * m) m : ℝ) * (2 * (1 / 2 : ℝ) ^ a) := + mul_lt_mul_of_pos_left + ((geom_sum_Ico_lt_of_lt_one (by norm_num) (by norm_num) a _).trans_eq (by ring)) hM + _ = 2 * (1 / 2 : ℝ) ^ a * (M (19 * m) m : ℝ) := by ring + +/-- At `L = 19 * m`, the sum of `beta L m d` over `0 ≤ d ≤ m` with `m < 9 * d` is strictly +below `2 * saving m * M L m`. This includes `m = 0`, when the sum is empty. -/ +theorem sum_beta_large_orders_lt (m : ℕ) : + ((∑ d ∈ (range (m + 1)).filter (fun d => m < 9 * d), beta (19 * m) m d : ℕ) : ℝ) + < 2 * saving m * (M (19 * m) m : ℝ) := by + -- The orders d > m/9 are the orders d ≥ ⌊m/9⌋ + 1. + have hs : (range (m + 1)).filter (fun d => m < 9 * d) ⊆ Ico (m / 9 + 1) (m + 1) := by + intro d hd + obtain ⟨hdm, hd9⟩ := mem_filter.mp hd + rw [mem_range] at hdm + rw [mem_Ico] + omega + refine (sum_beta_lt_of_subset m (m / 9 + 1) _ hs).trans_le ?_ + have hpow : (1 / 2 : ℝ) ^ (m / 9 + 1) ≤ saving m := by + have hfloor : (m : ℝ) ≤ 9 * ((m / 9 + 1 : ℕ) : ℝ) := by + exact_mod_cast (by omega : m ≤ 9 * (m / 9 + 1)) + rw [show (1 / 2 : ℝ) ^ (m / 9 + 1) = (2 : ℝ) ^ (-((m / 9 + 1 : ℕ) : ℝ)) by + rw [Real.rpow_neg (by norm_num), Real.rpow_natCast, one_div, inv_pow]] + exact Real.rpow_le_rpow_of_exponent_le (by norm_num) (by linarith) + -- Multiply `hpow` by `2` and by `M`. + gcongr + +/-- One term of the display in the proof of Lemma 11: `binom(m, d) 9^d = binom(m, d) 72^d 8^{-d}`, +which is at most `72^{m/9} binom(m, d) 8^{-d}` for `d ≤ m/9`. The factor `1^{m-d}` makes the right +side a term of the expansion of `(1/8 + 1)^m`. -/ +private theorem alpha_le (m d : ℕ) (hd : 9 * d ≤ m) : + (alpha m d : ℝ) + ≤ (72 : ℝ) ^ ((m : ℝ) / 9) * ((1 / 8 : ℝ) ^ d * 1 ^ (m - d) * (m.choose d : ℝ)) := by + have hd' : (9 : ℝ) * d ≤ m := by exact_mod_cast hd + have hpow : (72 : ℝ) ^ d ≤ (72 : ℝ) ^ ((m : ℝ) / 9) := by + rw [← Real.rpow_natCast (72 : ℝ) d] + exact Real.rpow_le_rpow_of_exponent_le (by norm_num) (by linarith) + calc (alpha m d : ℝ) = (72 : ℝ) ^ d * ((1 / 8 : ℝ) ^ d * 1 ^ (m - d) * (m.choose d : ℝ)) := by + rw [alpha, one_pow, mul_one, ← mul_assoc, ← mul_pow] + push_cast + norm_num [mul_comm] + _ ≤ _ := by gcongr + +/-- The first two steps of the display in the proof of Lemma 11: +`∑_{d ≤ m/9} binom(m, d) 9^d ≤ 72^{m/9} ∑_{d=0}^{m} binom(m, d) 8^{-d} = 72^{m/9} (9/8)^m`. -/ +private theorem sum_alpha_small_orders_le_aux (m : ℕ) : + ((∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), alpha m d : ℕ) : ℝ) + ≤ (72 : ℝ) ^ ((m : ℝ) / 9) * (9 / 8 : ℝ) ^ m := by + push_cast + calc ∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), (alpha m d : ℝ) + ≤ ∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), + (72 : ℝ) ^ ((m : ℝ) / 9) * ((1 / 8 : ℝ) ^ d * 1 ^ (m - d) * (m.choose d : ℝ)) := + sum_le_sum fun d hd => alpha_le m d (mem_filter.mp hd).2 + _ ≤ ∑ d ∈ range (m + 1), + (72 : ℝ) ^ ((m : ℝ) / 9) * ((1 / 8 : ℝ) ^ d * 1 ^ (m - d) * (m.choose d : ℝ)) := + sum_le_sum_of_subset_of_nonneg (filter_subset _ _) fun _ _ _ => by positivity + _ = (72 : ℝ) ^ ((m : ℝ) / 9) * (9 / 8 : ℝ) ^ m := by + rw [← mul_sum, ← add_pow] + norm_num + +/-- For `m ≥ 1`, the sum of `alpha m d` over `0 ≤ d ≤ m` with `9 * d ≤ m` is strictly below `2 +^ m * saving m`. For `m = 0`, the two quantities are equal, so the strict inequality needs +the stated positivity hypothesis. -/ +theorem sum_alpha_small_orders_lt (m : ℕ) (hm : 1 ≤ m) : + ((∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), alpha m d : ℕ) : ℝ) + < (2 : ℝ) ^ m * saving m := by + have hm9 : (0 : ℝ) < (m : ℝ) / 9 := by + have : (0 : ℝ) < (m : ℝ) := by exact_mod_cast hm + positivity + refine (sum_alpha_small_orders_le_aux m).trans_lt ?_ + -- 72^{m/9} (9/8)^m = (72 · (9/8)^9)^{m/9} + have hleft : (72 : ℝ) ^ ((m : ℝ) / 9) * (9 / 8 : ℝ) ^ m + = (72 * (9 / 8 : ℝ) ^ 9) ^ ((m : ℝ) / 9) := by + rw [Real.mul_rpow (by norm_num) (by positivity), ← Real.rpow_natCast (9 / 8 : ℝ) 9, + ← Real.rpow_mul (by norm_num), ← Real.rpow_natCast (9 / 8 : ℝ) m] + congr 2 + push_cast + ring + -- 256^{m/9} = 2^m · 2^{-m/9} + have hright : (256 : ℝ) ^ ((m : ℝ) / 9) = (2 : ℝ) ^ m * saving m := by + rw [saving, show (256 : ℝ) = (2 : ℝ) ^ ((8 : ℕ) : ℝ) by rw [Real.rpow_natCast]; norm_num, + ← Real.rpow_mul (by norm_num), ← Real.rpow_natCast (2 : ℝ) m, ← Real.rpow_add (by norm_num)] + congr 1 + push_cast + ring + rw [hleft, ← hright] + -- `72 · (9/8)^9 < 208 < 256`. + exact Real.rpow_lt_rpow (by positivity) (by norm_num) hm9 + +/-- The bound of `sum_alpha_small_orders_lt` with `≤`, which also holds for `m = 0`. -/ +private theorem sum_alpha_small_orders_le (m : ℕ) : + ((∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), alpha m d : ℕ) : ℝ) + ≤ (2 : ℝ) ^ m * saving m := by + rcases Nat.eq_zero_or_pos m with rfl | hm + · simp [alpha, saving] + · exact (sum_alpha_small_orders_lt m hm).le + +/-- Split the order sum at `m / 9`, applying the query-dependent estimate below the cutoff +and the layout-dependent estimate above it. -/ +private theorem card_leaves_le_split {L m : ℕ} (U : Finset (OutStr L)) + (hU : ∀ w ∈ U, (innerSetO w).card = m) : + (leaves U).card + ≤ U.card * ∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), alpha m d + + ∑ d ∈ (range (m + 1)).filter (fun d => m < 9 * d), beta L m d := by + refine (card_leaves_le_sum_min U hU).trans ?_ + rw [← sum_filter_add_sum_filter_not (range (m + 1)) (fun d => 9 * d ≤ m), mul_sum] + refine add_le_add (sum_le_sum fun d _ => min_le_left _ _) ?_ + rw [filter_congr fun d _ => (not_le (a := 9 * d) (b := m))] + exact sum_le_sum fun d _ => min_le_right _ _ + +/-- For strings of length `19 * m`, each with `m` inner positions, the contributing leaf set +has cardinality at most `saving m * (2 ^ m * U.card + 2 * M (19 * m) m)`. This real bound +permits `m = 0` and `U = ∅`; it makes no execution-cost claim (cf. [AlmanWilliams2026, +Lemma 11]). -/ +theorem card_leaves_le_saving {m : ℕ} (U : Finset (OutStr (19 * m))) + (hU : ∀ w ∈ U, (innerSetO w).card = m) : + ((leaves U).card : ℝ) ≤ saving m * ((2 : ℝ) ^ m * (U.card : ℝ) + 2 * (M (19 * m) m : ℝ)) := by + have hsplit : ((leaves U).card : ℝ) + ≤ (U.card : ℝ) * ((∑ d ∈ (range (m + 1)).filter (fun d => 9 * d ≤ m), alpha m d : ℕ) : ℝ) + + ((∑ d ∈ (range (m + 1)).filter (fun d => m < 9 * d), beta (19 * m) m d : ℕ) : ℝ) := by + exact_mod_cast card_leaves_le_split U hU + have hsmall := + mul_le_mul_of_nonneg_left (sum_alpha_small_orders_le m) (Nat.cast_nonneg (α := ℝ) U.card) + have hlarge := (sum_beta_large_orders_lt m).le + calc ((leaves U).card : ℝ) ≤ _ := hsplit + _ ≤ (U.card : ℝ) * ((2 : ℝ) ^ m * saving m) + 2 * saving m * (M (19 * m) m : ℝ) := + add_le_add hsmall hlarge + _ = _ := by ring + +/-- Combine the order-wise minimum bound on `leaves U` with the explicit real bound available +when `L = 19 * m`. Both conclusions assume only that each requested string has `m` inner +positions; the specialization is conditional on the displayed equality of parameters. -/ +theorem sparse_leaf_bounds {L m : ℕ} (U : Finset (OutStr L)) (hU : ∀ w ∈ U, (innerSetO w).card = + m) : + (leaves U).card ≤ ∑ d ∈ range (m + 1), min (U.card * alpha m d) (beta L m d) + ∧ (L = 19 * m → + ((leaves U).card : ℝ) + ≤ (2 : ℝ) ^ (-(m : ℝ) / 9) * ((2 : ℝ) ^ m * (U.card : ℝ) + 2 * (M L m : ℝ))) := by + refine ⟨card_leaves_le_sum_min U hU, ?_⟩ + rintro rfl + exact card_leaves_le_saving U hU + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/ThinProduct.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/ThinProduct.lean new file mode 100644 index 0000000000..cfe5b34c17 --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/ThinProduct.lean @@ -0,0 +1,205 @@ +/- +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.MatrixMultiplication.Schonhage.Parameters + +/-! +# Finite thin-product correctness and cardinality bounds + +## Main results + +* `thin_product_spec`: exact queried entries for arbitrary signed integer matrices. +* `encodings_total_le`, `total_leaves_le`, `total_passed_sets_le`: explicit finite totals. + +Correctness holds for arbitrary signed integer matrices and any layout. For `D = 4^m`, +`L = 19m`, `D^18 ≤ N` and `W.card * 2^m ≤ N²`, the encoding entry count is at most +`4 * D^(-1/18) * N²`, the sum of contributing-leaf cardinalities is at most +`33 * D^(-1/18) * N²`, and the sum of passed request-set cardinalities is at most +`(19m + 1) * 33 * D^(-1/18) * N²`. Here `N` is the original matrix dimension. +The definitions permit zero parameters and empty request sets; the bounds retain their +stated size and sparsity hypotheses. No native program or instruction bound is proved. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-- For `m ≤ L`, the band pair of any position `(I, J) : Fin N × Fin N` belongs to the finite +tile set. -/ +theorem tile_mem_tiles {L m N : ℕ} (hmL : m ≤ L) (I J : Fin N) : + (bandOf L m I, bandOf L m J) ∈ tiles L m N := + mem_product.mpr ⟨mem_range.mpr (bandOf_lt_numBands hmL N I I.2), + mem_range.mpr (bandOf_lt_numBands hmL N J J.2)⟩ + +/-- A requested position contributes its output string to `wantedStrings` for its band pair. -/ +theorem outStrOfPos_mem_wantedStrings {L m N : ℕ} (lay : Layout L m) (W : Finset (Fin N × Fin N)) + {I J : Fin N} (hIJ : (I, J) ∈ W) : + outStrOfPos lay I J ∈ wantedStrings lay W (bandOf L m I) (bandOf L m J) := + mem_image.mpr ⟨(I, J), mem_filter.mpr ⟨hIJ, rfl, rfl⟩, rfl⟩ + +/-- Every string requested in a tile has exactly `m` inner positions. -/ +theorem card_innerSetO_of_mem_wantedStrings {L m N : ℕ} (lay : Layout L m) + (W : Finset (Fin N × Fin N)) (β β' : ℕ) (w : OutStr L) (hw : w ∈ wantedStrings lay W β β') : + (innerSetO w).card = m := by + obtain ⟨IJ, -, rfl⟩ := mem_image.mp hw + exact card_innerSetO_outStrOfPos lay IJ.1 IJ.2 + +/-- For arbitrary signed integer matrices of shapes `N × D m` and `D m × N`, the pruned result +at any position in `W` equals the corresponding entry of `X * Y`. The encodings and +request set are those for that position's row and column bands. This holds for every +`Layout L m`, which includes `m ≤ L`, without the size or sparsity conditions used in the +counting bounds (cf. [AlmanWilliams2026, Theorem 5]). -/ +theorem thin_product_spec {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (W : Finset (Fin N × Fin N)) (I J : Fin N) + (hIJ : (I, J) ∈ W) : + Pruned.result L (encodingL (bandArrayL lay X (bandOf L m I))) + (encodingR (bandArrayR lay Y (bandOf L m J))) + (wantedStrings lay W (bandOf L m I) (bandOf L m J)) (outStrOfPos lay I J) + = (X * Y) I J := by + -- Lemma 10, with Lemma 7: `Pruned` returns the output of `Full` restricted to the set passed to + -- it. + rw [Pruned.result_eq_piecewise_full, Finset.piecewise_eq_of_mem _ _ _ + (outStrOfPos_mem_wantedStrings lay W hIJ)] + -- Lemma 9, through the tiling of Section 2.3.4. + exact Full.bandArray_eq_mul lay X Y I J + +/-- The real inequality `w ≤ N² / sqrt (D m)` is equivalent to the natural-number inequality +`w * 2 ^ m ≤ N ^ 2`, including `m = 0`. -/ +theorem le_sq_div_sqrt_D_iff (m N w : ℕ) : + (w : ℝ) ≤ (N : ℝ) ^ 2 / Real.sqrt (D m) ↔ w * 2 ^ m ≤ N ^ 2 := by + rw [← two_pow_eq_sqrt_D, le_div_iff₀ (by positivity)] + exact_mod_cast Iff.rfl + +/-- The number of bands times `N₀` is at most the padded size, because `K₀ ≥ 1`. -/ +private theorem numBands_mul_N0_le {L m : ℕ} (hmL : m ≤ L) (N : ℕ) : + numBands L m N * N0 L m ≤ padN L m N := + Nat.mul_le_mul_left _ (Nat.le_mul_of_pos_left _ (K0_pos hmL)) + +/-- At `L = 19 * m` with `D m ^ 18 ≤ N`, the encodings of all row-band and column-band arrays +have at most `4 * D^(-1/18) * N²` entries in total, where `D = D m`. The counted quantity +is `2 * numBands L m N * 10 ^ L`, and `N` is the original dimension. This is an entry +count, not a bound on instructions or encoding time. -/ +theorem encodings_total_le (m N : ℕ) (hN : D m ^ 18 ≤ N) : + ((2 * numBands (19 * m) m N * 10 ^ (19 * m) : ℕ) : ℝ) + ≤ 4 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2 := by + have hbands : (numBands (19 * m) m N : ℝ) * N0 (19 * m) m ≤ 2 * N := by + exact_mod_cast (numBands_mul_N0_le (by omega) N).trans (padN_le_two_mul_of_D_pow_le m N hN) + have heq6 : (10 : ℝ) ^ (19 * m) ≤ saving m * N * N0 (19 * m) m := encoding_size_le m N hN + have hc := saving_pos m + rw [← saving_eq_D_rpow] + push_cast + calc 2 * (numBands (19 * m) m N : ℝ) * 10 ^ (19 * m) + ≤ 2 * numBands (19 * m) m N * (saving m * N * N0 (19 * m) m) := by gcongr + _ = 2 * saving m * (N * (numBands (19 * m) m N * N0 (19 * m) m)) := by ring + _ ≤ 2 * saving m * (N * (2 * N)) := by gcongr + _ = 4 * saving m * (N : ℝ) ^ 2 := by ring + +/-- Summing the cardinalities of the requested-string sets over all tiles gives at most +`W.card`. -/ +theorem sum_card_wantedStrings_le {L m N : ℕ} (lay : Layout L m) (W : Finset (Fin N × Fin N)) : + ∑ T ∈ tiles L m N, (wantedStrings lay W T.1 T.2).card ≤ W.card := by + have hmap : ∀ IJ ∈ W, (bandOf L m IJ.1, bandOf L m IJ.2) ∈ tiles L m N := fun IJ _ => + tile_mem_tiles lay.hmL IJ.1 IJ.2 + rw [card_eq_sum_card_fiberwise hmap] + refine sum_le_sum fun T _ => ?_ + refine card_image_le.trans (le_of_eq ?_) + congr 1 + +/-- The sum of `M L m` over the finite tile set is at most `4 * padN L m N ^ 2`. -/ +theorem sum_M_le (L m N : ℕ) : ∑ _T ∈ tiles L m N, M L m ≤ 4 * padN L m N ^ 2 := by + rw [sum_const, smul_eq_mul] + exact card_tiles_mul_M_le L m N + +/-- At `L = 19 * m`, bound the sum of contributing-leaf cardinalities for the tile request +sets by `saving m * (2 ^ m * W.card + 2 * ∑ T ∈ tiles L m N, M L m)`. A leaf used in +distinct tiles is counted once for each such tile. -/ +theorem sum_card_leaves_le {m N : ℕ} (lay : Layout (19 * m) m) (W : Finset (Fin N × Fin N)) : + ((∑ T ∈ tiles (19 * m) m N, (leaves (wantedStrings lay W T.1 T.2)).card : ℕ) : ℝ) + ≤ saving m * ((2 : ℝ) ^ m * (W.card : ℝ) + + 2 * ((∑ _T ∈ tiles (19 * m) m N, M (19 * m) m : ℕ) : ℝ)) := by + have hc := saving_pos m + have hsum : ((∑ T ∈ tiles (19 * m) m N, (wantedStrings lay W T.1 T.2).card : ℕ) : ℝ) + ≤ (W.card : ℝ) := by + exact_mod_cast sum_card_wantedStrings_le lay W + push_cast at hsum ⊢ + calc ∑ T ∈ tiles (19 * m) m N, ((leaves (wantedStrings lay W T.1 T.2)).card : ℝ) + ≤ ∑ T ∈ tiles (19 * m) m N, saving m + * ((2 : ℝ) ^ m * ((wantedStrings lay W T.1 T.2).card : ℝ) + 2 * (M (19 * m) m : ℝ)) := + sum_le_sum fun T _ => + card_leaves_le_saving _ (card_innerSetO_of_mem_wantedStrings lay W T.1 T.2) + _ = saving m + * ((2 : ℝ) ^ m * ∑ T ∈ tiles (19 * m) m N, ((wantedStrings lay W T.1 T.2).card : ℝ) + + 2 * ∑ _T ∈ tiles (19 * m) m N, (M (19 * m) m : ℝ)) := by + rw [← mul_sum, sum_add_distrib, ← mul_sum, ← mul_sum] + _ ≤ _ := by gcongr + +/-- At `L = 19 * m`, assume `D m ^ 18 ≤ N` and `W.card * 2 ^ m ≤ N ^ 2`. The sum of +contributing-leaf cardinalities over all tiles is then at most `33 * D^(-1/18) * N²`, +where `D = D m` and `N` is the original dimension. This counts leaves with their tile +multiplicities, not executed instructions. -/ +theorem total_leaves_le {m N : ℕ} (lay : Layout (19 * m) m) (W : Finset (Fin N × Fin N)) + (hN : D m ^ 18 ≤ N) (hW : W.card * 2 ^ m ≤ N ^ 2) : + ((∑ T ∈ tiles (19 * m) m N, (leaves (wantedStrings lay W T.1 T.2)).card : ℕ) : ℝ) + ≤ 33 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2 := by + have hc := saving_pos m + have hWN : (2 : ℝ) ^ m * (W.card : ℝ) ≤ (N : ℝ) ^ 2 := by + rw [mul_comm] + exact_mod_cast hW + have hM : ((∑ _T ∈ tiles (19 * m) m N, M (19 * m) m : ℕ) : ℝ) ≤ 16 * (N : ℝ) ^ 2 := by + have hpad := Nat.pow_le_pow_left (padN_le_two_mul_of_D_pow_le m N hN) 2 + exact_mod_cast (sum_M_le (19 * m) m N).trans (by linarith) + refine (sum_card_leaves_le lay W).trans ?_ + rw [← saving_eq_D_rpow] + calc saving m + * ((2 : ℝ) ^ m * (W.card : ℝ) + 2 * ((∑ _T ∈ tiles (19 * m) m N, M (19 * m) m : ℕ) : ℝ)) + ≤ saving m * ((N : ℝ) ^ 2 + 2 * (16 * (N : ℝ) ^ 2)) := by gcongr + _ = 33 * saving m * (N : ℝ) ^ 2 := by ring + +/-- Lemma 10 summed over the tiles: the total size of the sets handled by the pruned recursions is +at most `L + 1` times the total number of leaves. -/ +private theorem sum_totalSize_le {m N : ℕ} (lay : Layout (19 * m) m) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (W : Finset (Fin N × Fin N)) : + ((∑ T ∈ tiles (19 * m) m N, + Pruned.totalSize (encodingL (bandArrayL lay X T.1)) (encodingR (bandArrayR lay Y T.2)) + (wantedStrings lay W T.1 T.2) : ℕ) : ℝ) + ≤ ((19 * m + 1 : ℕ) : ℝ) + * ((∑ T ∈ tiles (19 * m) m N, (leaves (wantedStrings lay W T.1 T.2)).card : ℕ) : ℝ) := by + have h : ∑ T ∈ tiles (19 * m) m N, + Pruned.totalSize (encodingL (bandArrayL lay X T.1)) (encodingR (bandArrayR lay Y T.2)) + (wantedStrings lay W T.1 T.2) + ≤ (19 * m + 1) * ∑ T ∈ tiles (19 * m) m N, (leaves (wantedStrings lay W T.1 T.2)).card := by + rw [mul_sum] + exact sum_le_sum fun T _ => Pruned.totalSize_le _ _ _ + exact_mod_cast h + +/-- At `L = 19 * m`, assume `D m ^ 18 ≤ N` and `W.card * 2 ^ m ≤ N ^ 2`. Summing +`Pruned.totalSize` over all tiles gives at most `(19 * m + 1) * (33 * D^(-1/18) * N²)`, +where `D = D m` and `N` is the original dimension. The quantity sums request-set +cardinalities through the pruned recursions; it is not an instruction count. -/ +theorem total_passed_sets_le {m N : ℕ} (lay : Layout (19 * m) m) (X : Matrix (Fin N) (Fin (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (W : Finset (Fin N × Fin N)) (hN : D m ^ 18 ≤ N) + (hW : W.card * 2 ^ m ≤ N ^ 2) : + ((∑ T ∈ tiles (19 * m) m N, + Pruned.totalSize (encodingL (bandArrayL lay X T.1)) (encodingR (bandArrayR lay Y T.2)) + (wantedStrings lay W T.1 T.2) : ℕ) : ℝ) + ≤ ((19 * m + 1 : ℕ) : ℝ) * (33 * (D m : ℝ) ^ (-(1 / 18 : ℝ)) * (N : ℝ) ^ 2) := + (sum_totalSize_le lay X Y W).trans + (mul_le_mul_of_nonneg_left (total_leaves_le lay W hN hW) (Nat.cast_nonneg _)) + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Algorithms/MatrixMultiplication/Schonhage/Tiling.lean b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Tiling.lean new file mode 100644 index 0000000000..980e8f750d --- /dev/null +++ b/Complexity/Algorithms/MatrixMultiplication/Schonhage/Tiling.lean @@ -0,0 +1,237 @@ +/- +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.MatrixMultiplication.Schonhage.Batch +public import Complexity.Algorithms.MatrixMultiplication.Schonhage.Layout +import Mathlib.Tactic.Ring + +/-! +# Padding and tiling correctness + +## Main results + +* `nonempty_layout`: valid layouts exist whenever `m ≤ L`. +* `Full.bandArray_eq_mul`: product entries as `Full.result` values of band arrays. +* `outStrOfPos_injective`: positions in a tile have distinct output strings. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Cslib.Algorithms.MatrixMultiplication.Schonhage + +/-! ### The numbers -/ + +/-- The square of the natural square root `K0 L m` is at most `K L m`. -/ +theorem K0_sq_le_K (L m : ℕ) : K0 L m ^ 2 ≤ K L m := + Nat.sqrt_le' _ + +/-- The binomial count `K L m` is at most four times `K0 L m ^ 2`. -/ +theorem K_le_four_mul_K0_sq (L m : ℕ) : K L m ≤ 4 * K0 L m ^ 2 := by + have hlt : K L m < (K0 L m + 1) ^ 2 := Nat.lt_succ_sqrt' _ + rcases Nat.eq_zero_or_pos (K0 L m) with h0 | h0 + · rw [h0] at hlt ⊢ + omega + · calc K L m ≤ (K0 L m + 1) ^ 2 := hlt.le + _ ≤ (2 * K0 L m) ^ 2 := Nat.pow_le_pow_left (by omega) 2 + _ = 4 * K0 L m ^ 2 := by ring + +/-- A band has at least one row. -/ +theorem bandSize_pos {L m : ℕ} (hmL : m ≤ L) : 0 < bandSize L m := + Nat.mul_pos (K0_pos hmL) (N0_pos L m) + +/-- The padded dimension is a multiple of the band size. -/ +theorem bandSize_dvd_padN (L m N : ℕ) : bandSize L m ∣ padN L m N := + Dvd.intro_left _ rfl + +/-- For `m ≤ L`, padding covers all `N` original positions. The parameter condition ensures a +positive band size. -/ +theorem le_padN {L m : ℕ} (hmL : m ≤ L) (N : ℕ) : N ≤ padN L m N := by + have h := Nat.lt_div_mul_add (a := N + bandSize L m - 1) (bandSize_pos hmL) + unfold padN numBands + omega + +/-- For `m ≤ L`, the padded dimension is strictly less than `N + bandSize L m`, including `N = +0`. -/ +theorem padN_lt_add_bandSize {L m : ℕ} (hmL : m ≤ L) (N : ℕ) : + padN L m N < N + bandSize L m := by + have h := Nat.div_mul_le_self (N + bandSize L m - 1) (bandSize L m) + have hb := bandSize_pos hmL + unfold padN numBands + omega + +/-- If the band size is at most `N`, the padded dimension is at most `2 * N`. -/ +theorem padN_le_two_mul (L m N : ℕ) (hN : K0 L m * N0 L m ≤ N) : padN L m N ≤ 2 * N := by + rcases Nat.lt_or_ge L m with hLm | hmL + · -- for `m > L` there are no subsets, `K₀ = 0`, and the padded size is 0 + have hK : K0 L m = 0 := by simp [K0, K, Nat.choose_eq_zero_of_lt hLm] + simp [padN, bandSize, hK] + · have h := padN_lt_add_bandSize hmL N + have hN' : bandSize L m ≤ N := hN + omega + +/-- The number of tiles is the square of the number of bands. -/ +theorem card_tiles (L m N : ℕ) : (tiles L m N).card = numBands L m N ^ 2 := by + simp [tiles, sq] + +/-- The tile count multiplied by `M L m` is at most `4 * padN L m N ^ 2`. The bound is stated +without division and requires no positivity hypothesis. -/ +theorem card_tiles_mul_M_le (L m N : ℕ) : (tiles L m N).card * M L m ≤ 4 * padN L m N ^ 2 := by + rw [card_tiles] + calc numBands L m N ^ 2 * M L m + = numBands L m N ^ 2 * N0 L m ^ 2 * K L m := by unfold M; ring + _ ≤ numBands L m N ^ 2 * N0 L m ^ 2 * (4 * K0 L m ^ 2) := + Nat.mul_le_mul_left _ (K_le_four_mul_K0_sq L m) + _ = 4 * padN L m N ^ 2 := by unfold padN bandSize; ring + +/-! ### The table of subsets and the matrices of a tile -/ + +/-- Indexing data `Layout L m` exist for every `m ≤ L`, including `m = 0`. -/ +theorem nonempty_layout {L m : ℕ} (hmL : m ≤ L) : Nonempty (Layout L m) := by + have hcard : Fintype.card (Fin (K0 L m) × Fin (K0 L m)) + ≤ Fintype.card {Q : Finset (Fin L) // Q.card = m} := by + rw [Fintype.card_prod, Fintype.card_fin, Fintype.card_subtype, card_subsets_eq_K, ← sq] + exact K0_sq_le_K L m + obtain ⟨e⟩ := Function.Embedding.nonempty_of_card_le hcard + exact ⟨{ hmL := hmL + rowIdx := (Fintype.equivFinOfCardEq (card_outerStr L m)).symm + colIdx := (Fintype.equivFinOfCardEq (card_outerStr L m)).symm + innerIdx := (Fintype.equivFinOfCardEq (card_innerStr m)).symm + table := fun gh => (e gh).1 + table_card := fun gh => (e gh).2 + table_injective := fun _ _ h => e.injective (Subtype.ext h) }⟩ + +/-- The subsets in the table are distinct, so a sum over the table that selects the subset of the +block product `gh` has one term. -/ +private theorem sum_table_eq {L m : ℕ} {A : Type*} [AddCommMonoid A] (lay : Layout L m) + (F : Fin (K0 L m) × Fin (K0 L m) → A) (gh : Fin (K0 L m) × Fin (K0 L m)) : + (∑ gh', if lay.table gh' = lay.table gh then F gh' else 0) = F gh := by + simp [lay.table_injective.eq_iff] + +/-- At the table entry `(g, h)`, the left family selects row block `β * K0 L m + g`. -/ +@[simp] +theorem bandFamilyL_table {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) (β : ℕ) + (g h : Fin (K0 L m)) : + bandFamilyL lay X β (lay.table (g, h)) = rowBlock lay X (β * K0 L m + g) := + sum_table_eq lay _ (g, h) + +/-- The matrix `Y_Q` of the block product in position `(g, h)` is the `h`-th column block of the +column band. -/ +@[simp] +theorem bandFamilyR_table {L m N : ℕ} (lay : Layout L m) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (β : ℕ) + (g h : Fin (K0 L m)) : + bandFamilyR lay Y β (lay.table (g, h)) = colBlock lay Y (β * K0 L m + h) := + sum_table_eq lay _ (g, h) + +/-- The left family is zero at any level set outside the table image. -/ +theorem bandFamilyL_eq_zero {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin (D m)) ℤ) + (β : ℕ) {Q : Finset (Fin L)} (hQ : ∀ gh, lay.table gh ≠ Q) : bandFamilyL lay X β Q = 0 := + sum_eq_zero fun gh _ => ite_eq_right (hQ gh) + +/-- The right family is zero at any level set outside the table image. -/ +theorem bandFamilyR_eq_zero {L m N : ℕ} (lay : Layout L m) (Y : Matrix (Fin (D m)) (Fin N) ℤ) + (β : ℕ) {Q : Finset (Fin L)} (hQ : ∀ gh, lay.table gh ≠ Q) : bandFamilyR lay Y β Q = 0 := + sum_eq_zero fun gh _ => ite_eq_right (hQ gh) + +/-! ### Product identities for tile arrays -/ + +/-- For the band arrays of tile `(β, β')`, `Full.result` at the output string for table entry +`(g, h)` equals the product entry of row block `β * K0 L m + g` and column block `β' * K0 +L m + h`. -/ +theorem Full.bandArray_eq_block_mul {L m N : ℕ} (lay : Layout L m) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (β β' : ℕ) + (g h : Fin (K0 L m)) (r c : OuterStr L m) : + Full.result L (bandArrayL lay X β) (bandArrayR lay Y β') + (outStrOf (lay.table (g, h)) (lay.table_card (g, h)) r c) + = (rowBlock lay X (β * K0 L m + g) * colBlock lay Y (β' * K0 L m + h)) r c := by + rw [bandArrayL, bandArrayR, Full.arrays_eq_mul, bandFamilyL_table, bandFamilyR_table] + +/-- A row (or column) is recovered from its band, the position of its block in the band and its +position in the block. -/ +private theorem position_decomp {L m : ℕ} (lay : Layout L m) (I : ℕ) : + (bandOf L m I * K0 L m + (blockOf lay I : ℕ)) * N0 L m + (offsetOf L m I : ℕ) = I := by + have hband : bandOf L m I = I / N0 L m / K0 L m := by + rw [bandOf, bandSize, Nat.mul_comm, Nat.div_div_eq_div_mul] + have hblock : bandOf L m I * K0 L m + (blockOf lay I : ℕ) = I / N0 L m := by + rw [hband] + exact Nat.div_add_mod' _ _ + rw [hblock] + exact Nat.div_add_mod' _ _ + +/-- The band arrays for the tile containing an original position `(I, J)` give `Full.result` +equal to `(X * Y) I J` at `outStrOfPos lay I J`. The matrices may have arbitrary signed +integer entries. -/ +theorem Full.bandArray_eq_mul {L m N : ℕ} (lay : Layout L m) + (X : Matrix (Fin N) (Fin (D m)) ℤ) (Y : Matrix (Fin (D m)) (Fin N) ℤ) (I J : Fin N) : + Full.result L (bandArrayL lay X (bandOf L m I)) (bandArrayR lay Y (bandOf L m J)) (outStrOfPos + lay I J) + = (X * Y) I J := by + unfold outStrOfPos + rw [Full.bandArray_eq_block_mul, Matrix.mul_apply, Matrix.mul_apply, + ← Equiv.sum_comp lay.innerIdx] + -- Both sides are sums over the `D` inner indices; the summands agree because row `I` is recovered + -- from its band, block and place in the block, and likewise column `J`. + refine sum_congr rfl fun k _ => ?_ + simp only [rowBlock, colBlock, Equiv.symm_apply_apply, position_decomp, padRows, padCols, + I.isLt, J.isLt, dite_eq_left, Fin.eta] + +/-- At the output string of an original position `(I, J)`, the bilinear product of its two +band arrays equals `(X * Y) I J`. -/ +theorem bilinearProduct_bandArray_eq_mul {L m N : ℕ} (lay : Layout L m) (X : Matrix (Fin N) (Fin + (D m)) ℤ) + (Y : Matrix (Fin (D m)) (Fin N) ℤ) (I J : Fin N) : + bilinearProduct (bandArrayL lay X (bandOf L m I)) (bandArrayR lay Y (bandOf L m J)) + (outStrOfPos lay I J) + = (X * Y) I J := by + rw [← Full.result_eq_bilinearProduct] + exact Full.bandArray_eq_mul lay X Y I J + +/-- For `m ≤ L`, any index `I < N` has band index below `numBands L m N`. -/ +theorem bandOf_lt_numBands {L m : ℕ} (hmL : m ≤ L) (N I : ℕ) (hI : I < N) : + bandOf L m I < numBands L m N := by + rw [bandOf, Nat.div_lt_iff_lt_mul (bandSize_pos hmL)] + exact hI.trans_le (le_padN hmL N) + +/-! ### The output string of a position -/ + +/-- The inner set of a position's output string is the table entry indexed by its two +within-band block indices. -/ +theorem innerSetO_outStrOfPos {L m : ℕ} (lay : Layout L m) (I J : ℕ) : + innerSetO (outStrOfPos lay I J) = lay.table (blockOf lay I, blockOf lay J) := + innerSetO_outStrOf _ _ _ _ + +/-- The inner set of the output string of a position has exactly `m` elements. -/ +theorem card_innerSetO_outStrOfPos {L m : ℕ} (lay : Layout L m) (I J : ℕ) : + (innerSetO (outStrOfPos lay I J)).card = m := + innerSetO_outStrOfPos lay I J ▸ lay.table_card _ + +/-- Within a fixed pair of row and column bands, equality of output strings implies equality +of both natural-number position indices. -/ +theorem outStrOfPos_injective {L m : ℕ} (lay : Layout L m) (I J I' J' : ℕ) + (hI : bandOf L m I = bandOf L m I') (hJ : bandOf L m J = bandOf L m J') + (h : outStrOfPos lay I J = outStrOfPos lay I' J') : I = I' ∧ J = J' := by + -- The two strings have the same inner set, row and column; so the two positions lie in the same + -- block product and at the same place in it. + obtain ⟨hQ, hr, hc⟩ := outStrOf_inj.mp h + obtain ⟨hblockI, hblockJ⟩ := Prod.mk.inj (lay.table_injective hQ) + constructor + · rw [← position_decomp lay I, ← position_decomp lay I', hI, hblockI, + lay.rowIdx.injective hr] + · rw [← position_decomp lay J, ← position_decomp lay J', hJ, hblockJ, + lay.colIdx.injective hc] + +end Cslib.Algorithms.MatrixMultiplication.Schonhage diff --git a/Complexity/Foundations/Algebra/FiniteMagnitude.lean b/Complexity/Foundations/Algebra/FiniteMagnitude.lean new file mode 100644 index 0000000000..4bf60b86e2 --- /dev/null +++ b/Complexity/Foundations/Algebra/FiniteMagnitude.lean @@ -0,0 +1,53 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Algebra.Order.BigOperators.Ring.Finset + +/-! +# Magnitudes of finite sums and products + +## Main results + +* `Finset.abs_sum_mul_le`: bounded signed linear combinations. +* `Finset.abs_prod_le_one`: products of coefficients bounded by one. + +## References + +Adapted selectively from +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +namespace Finset + +variable {ι R : Type*} + +/-- If `|f i| ≤ A` for all `i ∈ s`, then `|∑ i ∈ s, f i| ≤ s.card * A`. -/ +theorem abs_sum_le_card_mul [Ring R] [LinearOrder R] [IsOrderedRing R] (s : Finset ι) {f : ι → R} + {A : R} (h : ∀ i ∈ s, |f i| ≤ A) : |∑ i ∈ s, f i| ≤ s.card * A := + (abs_sum_le_sum_abs f s).trans ((sum_le_card_nsmul s _ A h).trans_eq (nsmul_eq_mul _ _)) + +/-- If `|c i| ≤ 1` and `|a i| ≤ A` for all `i ∈ s`, then `|∑ i ∈ s, c i * a i| ≤ s.card * A`. -/ +theorem abs_sum_mul_le [Ring R] [LinearOrder R] [IsOrderedRing R] (s : Finset ι) {c a : ι → R} + {A : R} (hc : ∀ i ∈ s, |c i| ≤ 1) (ha : ∀ i ∈ s, |a i| ≤ A) : + |∑ i ∈ s, c i * a i| ≤ s.card * A := + abs_sum_le_card_mul s fun i hi => + (abs_mul (c i) (a i)).trans_le + ((mul_le_of_le_one_left (abs_nonneg _) (hc i hi)).trans (ha i hi)) + +/-- If `|f i| ≤ 1` for all `i ∈ s`, then `|∏ i ∈ s, f i| ≤ 1`. -/ +theorem abs_prod_le_one [CommRing R] [LinearOrder R] [IsStrictOrderedRing R] (s : Finset ι) + {f : ι → R} (h : ∀ i ∈ s, |f i| ≤ 1) : |∏ i ∈ s, f i| ≤ 1 := + (abs_prod s f).trans_le (prod_le_one₀ (fun _ _ => abs_nonneg _) h) + + +end Finset diff --git a/Complexity/Foundations/Algebra/GeometricTail.lean b/Complexity/Foundations/Algebra/GeometricTail.lean new file mode 100644 index 0000000000..ca762bc791 --- /dev/null +++ b/Complexity/Foundations/Algebra/GeometricTail.lean @@ -0,0 +1,40 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Algebra.Field.GeomSum +public import Mathlib.Algebra.Order.BigOperators.Ring.Finset + +/-! +# Strict bounds for finite geometric tails + +## Main results + +* `geom_sum_Ico_lt_of_lt_one`: the strict bound for `∑ i ∈ Ico m n, x ^ i` when `0 < x < 1`. + +## References + +Adapted selectively from +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +/-- For `0 < x < 1`, `∑ i ∈ Ico m n, x ^ i < x ^ m / (1 - x)`. Mathlib's +`geom_sum_Ico_le_of_lt_one` has `≤`, for `0 ≤ x`. -/ +theorem geom_sum_Ico_lt_of_lt_one {K : Type*} [Field K] [LinearOrder K] [IsStrictOrderedRing K] + {x : K} (hx : 0 < x) (hx1 : x < 1) (m n : ℕ) : + ∑ i ∈ Finset.Ico m n, x ^ i < x ^ m / (1 - x) := by + have hpos : 0 < 1 - x := sub_pos.2 hx1 + rcases le_total m n with hmn | hnm + · rw [geom_sum_Ico' hx1.ne hmn] + exact div_lt_div_of_pos_right (sub_lt_self _ (pow_pos hx n)) hpos + · rw [Finset.Ico_eq_empty_of_le hnm, Finset.sum_empty] + exact div_pos (pow_pos hx m) hpos diff --git a/Complexity/Foundations/Combinatorics/FunctionPlaces.lean b/Complexity/Foundations/Combinatorics/FunctionPlaces.lean new file mode 100644 index 0000000000..206c61f3da --- /dev/null +++ b/Complexity/Foundations/Combinatorics/FunctionPlaces.lean @@ -0,0 +1,140 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Algebra.Group.Action.Defs +public import Mathlib.Algebra.BigOperators.GroupWithZero.Finset +public import Mathlib.Data.Finset.Powerset +public import Mathlib.Data.Fintype.BigOperators + +/-! +# Counting functions by constrained positions + +## Main results + +* `Finset.card_pi_places_card`: count functions by the positions with a property. +* `Finset.card_places_card`: binomial count without pointwise constraints. + +## References + +Adapted selectively from +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. +-/ + +@[expose] public section + +namespace Finset + +variable {ι : Type*} [Fintype ι] [DecidableEq ι] + +/-- A product that takes the value `a` on `S` and `b` elsewhere. -/ +theorem prod_ite_mem_const {M : Type*} [CommMonoid M] (S : Finset ι) (a b : M) : + (∏ i, if i ∈ S then a else b) = a ^ #S * b ^ (Fintype.card ι - #S) := by + rw [prod_ite, prod_const, prod_const, filter_mem_eq_inter, univ_inter, filter_not, + filter_mem_eq_inter, univ_inter, ← compl_eq_univ_sdiff, card_compl] + +/-- A product with the factors `1` on `S ∩ Q`, `0` on `S \ Q`, `c` on `Q \ S` and `1` elsewhere: it +is `0` unless `S ⊆ Q`. -/ +theorem prod_ite_ite_subset {M : Type*} [CommMonoidWithZero M] (Q S : Finset ι) (c : M) : + (∏ i, if i ∈ S then (if i ∈ Q then 1 else 0) else (if i ∈ Q then c else 1)) + = if S ⊆ Q then c ^ (#Q - #S) else 0 := by + split_ifs with h + · have hterm : ∀ i, (if i ∈ S then (if i ∈ Q then 1 else 0) else (if i ∈ Q then c else 1)) + = if i ∈ Q \ S then c else 1 := fun i => by + by_cases hS : i ∈ S + · simp [hS, h hS] + · simp [hS] + simp_rw [hterm] + rw [prod_ite_mem_const, one_pow, mul_one, card_sdiff_of_subset h] + · obtain ⟨i, hiS, hiQ⟩ := not_subset.mp h + exact prod_eq_zero (mem_univ i) (by simp [hiS, hiQ]) + +/-- A product with the factor `0` on `T \ S` and `1` elsewhere: it is `0` unless `T ⊆ S`. -/ +theorem prod_ite_ite_superset {M : Type*} [CommMonoidWithZero M] (T S : Finset ι) : + (∏ i, if i ∈ S then 1 else (if i ∈ T then 0 else 1)) = if T ⊆ S then (1 : M) else 0 := by + split_ifs with h + · refine prod_eq_one fun i _ => ?_ + by_cases hS : i ∈ S + · simp [hS] + · have hT : i ∉ T := fun hT => hS (h hT) + simp [hS, hT] + · obtain ⟨i, hiT, hiS⟩ := not_subset.mp h + exact prod_eq_zero (mem_univ i) (by simp [hiS, hiT]) + +/-- The number of sets of `k` elements inside `Q`, each counted `c` times. -/ +theorem sum_powersetCard_ite_subset (Q : Finset ι) (k c : ℕ) : + (∑ S ∈ powersetCard k (univ : Finset ι), if S ⊆ Q then c else 0) = (#Q).choose k * c := by + rw [← sum_filter, sum_const, smul_eq_mul, ← card_powersetCard] + congr 2 + ext S + simp [mem_powersetCard, and_comm] + +/-- The number of sets of `k` elements that contain `T`: choose the complement outside `T`. -/ +theorem card_powersetCard_superset (T : Finset ι) {k : ℕ} (hk : k ≤ Fintype.card ι) : + #{S ∈ powersetCard k (univ : Finset ι) | T ⊆ S} = + (Fintype.card ι - #T).choose (Fintype.card ι - k) := by + rw [← card_compl, ← card_powersetCard] + refine card_bij' (fun S _ => Sᶜ) (fun S _ => Sᶜ) ?_ ?_ (fun _ _ => compl_compl _) + (fun _ _ => compl_compl _) + · intro S hS + simp only [mem_filter, mem_powersetCard, subset_univ, true_and] at hS + simp [mem_powersetCard, hS.1, hS.2, card_compl] + · intro S hS + simp only [mem_powersetCard] at hS + simp only [mem_filter, mem_powersetCard, subset_univ, true_and, card_compl, hS.2] + exact ⟨by omega, by rw [← compl_subset_compl, compl_compl]; exact hS.1⟩ + +/-! ## Functions by the set of places with a property -/ + +variable {α : Type*} [Fintype α] + +/-- Count the functions satisfying `f i ∈ A i` and having property `p` exactly on `S` by +multiplying the permitted choices at each position. -/ +theorem card_pi_places_eq [DecidableEq α] (A : ι → Finset α) (p : α → Prop) [DecidablePred p] + (S : Finset ι) : + #{f : ι → α | (∀ i, f i ∈ A i) ∧ ({i | p (f i)} : Finset ι) = S} = + ∏ i, if i ∈ S then #{a ∈ A i | p a} else #{a ∈ A i | ¬ p a} := by + have hmem (i : ι) (x : α) : (x ∈ if i ∈ S then {a ∈ A i | p a} else {a ∈ A i | ¬ p a}) ↔ + x ∈ A i ∧ (p x ↔ i ∈ S) := by + split_ifs with hi <;> simp [hi] + simp_rw [← apply_ite Finset.card, ← Fintype.card_piFinset] + congr 1 + ext f + simp only [mem_filter, mem_univ, true_and, Fintype.mem_piFinset, hmem, forall_and, + Finset.ext_iff] + +/-- Partition the functions satisfying `f i ∈ A i` and having property `p` at exactly `k` +positions by their size-`k` set of such positions, and sum the corresponding products of +choice counts. -/ +theorem card_pi_places_card [DecidableEq α] (A : ι → Finset α) (p : α → Prop) [DecidablePred p] + (k : ℕ) : + #{f : ι → α | (∀ i, f i ∈ A i) ∧ #{i | p (f i)} = k} = + ∑ S ∈ powersetCard k (univ : Finset ι), + ∏ i, if i ∈ S then #{a ∈ A i | p a} else #{a ∈ A i | ¬ p a} := by + rw [card_eq_sum_card_fiberwise (f := fun f : ι → α => ({i | p (f i)} : Finset ι)) + (t := powersetCard k (univ : Finset ι)) + fun f hf => mem_coe.2 (mem_powersetCard.2 ⟨subset_univ _, (mem_filter.1 hf).2.2⟩)] + refine sum_congr rfl fun S hS => ?_ + rw [← card_pi_places_eq, filter_filter] + exact congrArg card (filter_congr fun f _ => + ⟨fun h => ⟨h.1.1, h.2⟩, fun h => ⟨⟨h.1, h.2 ▸ (mem_powersetCard.1 hS).2⟩, h.2⟩⟩) + +/-- The number of functions `ι → α` that have `p` at exactly `k` places. -/ +theorem card_places_card (p : α → Prop) [DecidablePred p] (k : ℕ) : + #{f : ι → α | #{i | p (f i)} = k} = + (Fintype.card ι).choose k * (#{a | p a} ^ k * #{a | ¬ p a} ^ (Fintype.card ι - k)) := by + classical + have h := card_pi_places_card (ι := ι) (fun _ => (univ : Finset α)) p k + simp only [mem_univ, implies_true, true_and] at h + rw [h, sum_congr rfl fun S hS => by + rw [prod_ite_mem_const, (mem_powersetCard.1 hS).2], sum_const, card_powersetCard, card_univ, + smul_eq_mul] + +end Finset diff --git a/Complexity/Foundations/Data/Finset/OrderGlue.lean b/Complexity/Foundations/Data/Finset/OrderGlue.lean new file mode 100644 index 0000000000..9a61f690c7 --- /dev/null +++ b/Complexity/Foundations/Data/Finset/OrderGlue.lean @@ -0,0 +1,187 @@ +/- +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.Data.Finset.Sort + +/-! +# Gluing functions along ordered finite subsets + +## Main definitions + +* `Finset.glue`: join two strings along a subset and its complement. + +## Main results + +* `Finset.eq_glue_iff`, `Finset.glue_inj`: reconstruction and uniqueness. + +## References + +Adapted selectively from [AlmanWilliams2026, Section 2] and +[Anthropic's pinned formal source](https://github.com/anthropics/formal-math/tree/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec2). +Source pin: `16a5ce41304fcdaea4e0be7c132e1c361dc26da5`. These are finite mathematical +counts and identities; they give no native program or instruction-cost certificate. +-/ + +@[expose] public section + +open Finset + +namespace Finset + +/-- The complement of a size-`m` subset of `Fin L` has size `L - m`. -/ +theorem card_compl_of_card_eq {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card = m) : + Qᶜ.card = L - m := by + rw [Finset.card_compl, Fintype.card_fin, hQ] + +variable {L m : ℕ} (Q : Finset (Fin L)) (hQ : Q.card = m) + +/-! ### The levels of a set and of its complement, in order -/ + +/-- The `k`-th lowest level outside `Q` is not a level of `Q`. -/ +theorem outerLevel_notMem (k : Fin (L - m)) : (Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k ∉ + Q := + mem_compl.mp (Qᶜ.orderEmbOfFin_mem (card_compl_of_card_eq Q hQ) k) + +/-- Every level of `Q` is the `k`-th lowest level of `Q` for some `k`. -/ +theorem exists_innerLevel {ℓ : Fin L} (hℓ : ℓ ∈ Q) : ∃ k, Finset.orderEmbOfFin Q hQ k = ℓ := + (Q.range_orderEmbOfFin hQ).ge hℓ + +/-- Every level outside `Q` is the `k`-th lowest level outside `Q` for some `k`. -/ +theorem exists_outerLevel {ℓ : Fin L} (hℓ : ℓ ∉ Q) : ∃ k, (Qᶜ.orderEmbOfFin (card_compl_of_card_eq + Q hQ)) k = ℓ := + (Qᶜ.range_orderEmbOfFin (card_compl_of_card_eq Q hQ)).ge (mem_coe.mpr (mem_compl.mpr hℓ)) + +/-- The levels of `Q` in order depend only on the set `Q`. -/ +theorem innerLevel_congr {Q' : Finset (Fin L)} (h : Q' = Q) (hQ' : Q'.card = m) : + Finset.orderEmbOfFin Q' hQ' = Finset.orderEmbOfFin Q hQ := by + subst h + rfl + +/-- The levels outside `Q` in order depend only on the set `Q`. -/ +theorem outerLevel_congr {Q' : Finset (Fin L)} (h : Q' = Q) (hQ' : Q'.card = m) : + (Q'ᶜ.orderEmbOfFin (card_compl_of_card_eq Q' hQ')) = (Qᶜ.orderEmbOfFin (card_compl_of_card_eq + Q hQ)) := by + subst h + rfl + +/-- A property holds at every level if and only if it holds at the levels of `Q` and at the levels +outside `Q`. -/ +theorem forall_level_iff {P : Fin L → Prop} : + (∀ ℓ, P ℓ) ↔ (∀ k, P (Finset.orderEmbOfFin Q hQ k)) ∧ ∀ k, P ((Qᶜ.orderEmbOfFin + (card_compl_of_card_eq Q hQ)) k) := by + refine ⟨fun h => ⟨fun _ => h _, fun _ => h _⟩, fun ⟨hin, hout⟩ ℓ => ?_⟩ + by_cases hℓ : ℓ ∈ Q + · obtain ⟨k, rfl⟩ := exists_innerLevel Q hQ hℓ + exact hin k + · obtain ⟨k, rfl⟩ := exists_outerLevel Q hQ hℓ + exact hout k + +/-- A property cuts out the set `Q` if and only if it holds at the levels of `Q` and fails at the +other levels. -/ +theorem filter_eq_iff_forall_level {P : Fin L → Prop} [DecidablePred P] : + univ.filter P = Q ↔ (∀ k, P (Finset.orderEmbOfFin Q hQ k)) ∧ ∀ k, ¬ P ((Qᶜ.orderEmbOfFin + (card_compl_of_card_eq Q hQ)) k) := by + simp only [Finset.ext_iff, forall_level_iff Q hQ, mem_filter, mem_univ, true_and, + Finset.orderEmbOfFin_mem, outerLevel_notMem Q hQ, iff_true, iff_false] + +/-- A property that holds at the `m` levels of `Q`, and at `m` levels in all, cuts out `Q`. -/ +theorem filter_eq_of_forall_innerLevel {P : Fin L → Prop} [DecidablePred P] + (hin : ∀ k, P (Finset.orderEmbOfFin Q hQ k)) (hcard : (univ.filter P).card = m) : univ.filter + P = Q := by + refine (eq_of_subset_of_card_le (fun ℓ hℓ => ?_) (hcard.trans hQ.symm).le).symm + obtain ⟨k, rfl⟩ := exists_innerLevel Q hQ hℓ + exact mem_filter.mpr ⟨mem_univ _, hin k⟩ + +/-! ### Gluing two strings along a set of levels -/ + +variable {α : Type*} + +/-- The string that has the letters of `f` at the levels of `Q` and the letters of `g` at the other +levels, both in the order of the levels. -/ +def glue (f : Fin m → α) (g : Fin (L - m) → α) : Fin L → α := + fun ℓ => + if hℓ : ℓ ∈ Q then f ((Q.orderIsoOfFin hQ).symm ⟨ℓ, hℓ⟩) + else g ((Qᶜ.orderIsoOfFin (card_compl_of_card_eq Q hQ)).symm ⟨ℓ, mem_compl.mpr hℓ⟩) + +/-- At the `k`-th lowest level of `Q`, the glued string has the `k`-th letter of `f`. -/ +@[simp] +theorem glue_innerLevel (f : Fin m → α) (g : Fin (L - m) → α) (k : Fin m) : + glue Q hQ f g (Finset.orderEmbOfFin Q hQ k) = f k := by + rw [glue, dite_eq_left (Finset.orderEmbOfFin_mem Q hQ k)] + exact congrArg f ((OrderIso.symm_apply_eq _).mpr (Subtype.ext rfl)) + +/-- At the `k`-th lowest level outside `Q`, the glued string has the `k`-th letter of `g`. -/ +@[simp] +theorem glue_outerLevel (f : Fin m → α) (g : Fin (L - m) → α) (k : Fin (L - m)) : + glue Q hQ f g ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k) = g k := by + rw [glue, dite_eq_right (outerLevel_notMem Q hQ k)] + exact congrArg g ((OrderIso.symm_apply_eq _).mpr (Subtype.ext rfl)) + +/-- At the `k`-th lowest level of a set `Q'` that is equal to `Q`, the string glued along `Q` has +the `k`-th letter of `f`. -/ +theorem glue_innerLevel_of_eq (f : Fin m → α) (g : Fin (L - m) → α) {Q' : Finset (Fin L)} + (h : Q' = Q) (hQ' : Q'.card = m) (k : Fin m) : glue Q hQ f g (Finset.orderEmbOfFin Q' hQ' k) = + f k := by + rw [innerLevel_congr Q hQ h, glue_innerLevel] + +/-- At the `k`-th lowest level outside a set `Q'` that is equal to `Q`, the string glued along `Q` +has the `k`-th letter of `g`. -/ +theorem glue_outerLevel_of_eq (f : Fin m → α) (g : Fin (L - m) → α) {Q' : Finset (Fin L)} + (h : Q' = Q) (hQ' : Q'.card = m) (k : Fin (L - m)) : + glue Q hQ f g ((Q'ᶜ.orderEmbOfFin (card_compl_of_card_eq Q' hQ')) k) = g k := by + subst h + apply glue_outerLevel + +variable {f f' : Fin m → α} {g g' : Fin (L - m) → α} + +/-- The glued string is the only string with the letters of `f` at the levels of `Q` and the letters +of `g` at the other levels. -/ +theorem eq_glue_iff {u : Fin L → α} : + u = glue Q hQ f g + ↔ (∀ k, u (Finset.orderEmbOfFin Q hQ k) = f k) ∧ ∀ k, u ((Qᶜ.orderEmbOfFin + (card_compl_of_card_eq Q hQ)) k) = g k := by + simp only [funext_iff, forall_level_iff Q hQ, glue_innerLevel, glue_outerLevel] + +/-- The two strings can be read off the glued string. -/ +theorem glue_inj : glue Q hQ f g = glue Q hQ f' g' ↔ f = f' ∧ g = g' := by + rw [eq_glue_iff] + simp only [glue_innerLevel, glue_outerLevel, funext_iff] + +/-- A map applied to a glued string, letter by letter. -/ +theorem map_glue {β : Type*} (φ : α → β) : + (fun ℓ => φ (glue Q hQ f g ℓ)) = glue Q hQ (fun k => φ (f k)) (fun k => φ (g k)) := by + simp only [eq_glue_iff, glue_innerLevel, glue_outerLevel, implies_true, and_self] + +/-! ### Alphabets with inner and outer letters -/ + +variable (IsInner : α → Prop) [DecidablePred IsInner] + +/-- If the letters of `f` are inner and those of `g` are outer, then the glued string has its inner +letters exactly at the levels of `Q`. -/ +theorem filter_glue (hf : ∀ k, IsInner (f k)) (hg : ∀ k, ¬ IsInner (g k)) : + (univ.filter fun ℓ => IsInner (glue Q hQ f g ℓ)) = Q := by + simpa only [filter_eq_iff_forall_level Q hQ, glue_innerLevel, glue_outerLevel] using + And.intro hf hg + +/-- Reconstruct `u` by extracting and reinserting its inner and outer indices. The hypotheses +require that both index projections have the stated reconstruction property on their +respective classes of letters, and that the inner letters of `u` occur exactly on the +size-`m` set `Q`. -/ +theorem glue_index {I O : Type*} {inner : I → α} {outer : O → α} {innerIndex : α → I} + {outerIndex : α → O} (hI : ∀ s, IsInner s → inner (innerIndex s) = s) + (hO : ∀ s, ¬ IsInner s → outer (outerIndex s) = s) {u : Fin L → α} + (hu : (univ.filter fun ℓ => IsInner (u ℓ)) = Q) : + u = glue Q hQ (fun k => inner (innerIndex (u (Finset.orderEmbOfFin Q hQ k)))) + (fun k => outer (outerIndex (u ((Qᶜ.orderEmbOfFin (card_compl_of_card_eq Q hQ)) k)))) := by + obtain ⟨hin, hout⟩ := (filter_eq_iff_forall_level Q hQ).mp hu + exact (eq_glue_iff Q hQ).mpr ⟨fun k => (hI _ (hin k)).symm, fun k => (hO _ (hout k)).symm⟩ + + +end Finset diff --git a/Complexity/Foundations/Data/Int/LogBounds.lean b/Complexity/Foundations/Data/Int/LogBounds.lean new file mode 100644 index 0000000000..57fd66785e --- /dev/null +++ b/Complexity/Foundations/Data/Int/LogBounds.lean @@ -0,0 +1,34 @@ +/- +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.Algebra.Order.Ring.Int +public import Mathlib.Data.Nat.Log + +/-! +# Integer magnitude and binary logarithms + +## Main results + +* `Int.abs_lt_two_pow_log_succ`: a natural magnitude bound supplies enough binary digits. +-/ + +@[expose] public section + +namespace Int + +/-- An integer bounded in magnitude by `U` has magnitude below `2^(log₂ U+1)`, including +`U = 0`. This is the scalar consequence of `Nat.lt_pow_succ_log_self`. -/ +theorem abs_lt_two_pow_log_succ {U : ℕ} (x : ℤ) (hx : |x| ≤ U) : + |x| < 2 ^ (Nat.log 2 U + 1) := by + have h := Nat.lt_pow_succ_log_self (by decide : 1 < 2) U + have hcast : (U : ℤ) < 2 ^ (Nat.log 2 U + 1) := by exact_mod_cast h + exact hx.trans_lt hcast + +end Int diff --git a/Complexity/Foundations/Data/Matrix/IntegerBits.lean b/Complexity/Foundations/Data/Matrix/IntegerBits.lean index 5083b77db3..dd5f54d29d 100644 --- a/Complexity/Foundations/Data/Matrix/IntegerBits.lean +++ b/Complexity/Foundations/Data/Matrix/IntegerBits.lean @@ -10,7 +10,6 @@ module public import Complexity.Init public import Mathlib.Algebra.Order.Group.PosPart public import Mathlib.Data.Matrix.Basic -public import Mathlib.Data.Nat.Log import Mathlib.Tactic.Positivity import Mathlib.Tactic.Ring @@ -22,24 +21,25 @@ parts. These are algebraic identities, with no algorithm or running-time asserti ## Main definitions -* `IntMatrix.bit`: the binary digit matrix of a given weight. +* `Matrix.bit`: the binary digit matrix of a given weight. ## Main results -* `IntMatrix.binary`: expansion of a nonnegative bounded matrix. -* `IntMatrix.product`: expansion of a product into products of digit matrices. -* `IntMatrix.abs_lt_two_pow_log`: an explicit integer bit bound, including zero. +* `Matrix.eq_sum_two_pow_smul_bit`: expansion of a nonnegative bounded matrix. +* `Matrix.mul_eq_sum_two_pow_smul_bit_mul_bit`: expansion of a product into digit products. +* `Matrix.eq_map_posPart_sub_map_negPart`: reconstruction from pointwise signed parts. +* `Matrix.map_posPart_lt_and_map_negPart_lt`: bounds for both pointwise signed parts. ## References Adapted from [Anthropic's MatrixLanguage.lean](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/MatrixLanguage.lean), -`matBit` and `Footnote8.binary/product`. General index types and the integer bit bound are original -adaptation code. +`matBit` and `Footnote8.binary/product`. The adaptation permits general index types and uses +Mathlib's existing pointwise signed parts. -/ @[expose] public section -namespace Cslib.IntMatrix +namespace Cslib.Matrix /-- The binary digit of weight `2^i`, entry by entry. For negative entries this is the Euclidean remainder of the quotient; reconstruction below assumes nonnegative entries. -/ @@ -58,7 +58,7 @@ private theorem emod_two_pow_eq_sum_digits (x : ℤ) (β : ℕ) : ring /-- A matrix with entries in `[0, 2^β)` is the weighted sum of its first `β` digit matrices. -/ -theorem binary {m n : Type*} (X : Matrix m n ℤ) (β : ℕ) +theorem eq_sum_two_pow_smul_bit {m n : Type*} (X : Matrix m n ℤ) (β : ℕ) (hX : ∀ a j, 0 ≤ X a j ∧ X a j < 2 ^ β) : X = ∑ i ∈ Finset.range β, (2 : ℤ) ^ i • bit X i := by ext a j @@ -73,11 +73,12 @@ theorem bit_eq_zero_or_one {m n : Type*} (X : Matrix m n ℤ) (i : ℕ) (a : m) omega /-- Expand a product of nonnegative bounded matrices into products of binary digits. -/ -theorem product {m d n : Type*} [Fintype d] (X : Matrix m d ℤ) (Y : Matrix d n ℤ) +theorem mul_eq_sum_two_pow_smul_bit_mul_bit {m d n : Type*} [Fintype d] + (X : Matrix m d ℤ) (Y : Matrix d n ℤ) (β : ℕ) (hX : ∀ a j, 0 ≤ X a j ∧ X a j < 2 ^ β) (hY : ∀ j b, 0 ≤ Y j b ∧ Y j b < 2 ^ β) : X * Y = ∑ i ∈ Finset.range β, ∑ j ∈ Finset.range β, (2 : ℤ) ^ (i + j) • (bit X i * bit Y j) := by - conv_lhs => rw [binary X β hX, binary Y β hY] + conv_lhs => rw [eq_sum_two_pow_smul_bit X β hX, eq_sum_two_pow_smul_bit Y β hY] rw [Matrix.sum_mul] refine Finset.sum_congr rfl fun i _ => ?_ rw [Matrix.mul_sum] @@ -86,14 +87,14 @@ theorem product {m d n : Type*} [Fintype d] (X : Matrix m d ℤ) (Y : Matrix d n /-- Entrywise positive and negative parts reconstruct an integer matrix. `Matrix.map` reuses scalar signed parts without adding an ordering instance to matrices. -/ -theorem parts {m n : Type*} (X : Matrix m n ℤ) : +theorem eq_map_posPart_sub_map_negPart {m n : Type*} (X : Matrix m n ℤ) : X = X.map posPart - X.map negPart := by ext a j exact (posPart_sub_negPart (X a j)).symm /-- Both pointwise signed parts of an absolutely bounded matrix satisfy the same nonnegative bound. This uses the existing ordered-group operations, without introducing matrix variants. -/ -theorem parts_lt {m n : Type*} (X : Matrix m n ℤ) {β : ℕ} +theorem map_posPart_lt_and_map_negPart_lt {m n : Type*} (X : Matrix m n ℤ) {β : ℕ} (hX : ∀ a j, |X a j| < 2 ^ β) : (∀ a j, 0 ≤ (X.map posPart) a j ∧ (X.map posPart) a j < 2 ^ β) ∧ ∀ a j, 0 ≤ (X.map negPart) a j ∧ (X.map negPart) a j < 2 ^ β := by @@ -101,12 +102,4 @@ theorem parts_lt {m n : Type*} (X : Matrix m n ℤ) {β : ℕ} exact ⟨fun a j => ⟨le_max_right _ _, max_lt ((le_abs_self _).trans_lt (hX a j)) hpow⟩, fun a j => ⟨le_max_right _ _, max_lt ((neg_le_abs _).trans_lt (hX a j)) hpow⟩⟩ -/-- An integer bounded in magnitude by `U` needs at most `log₂ U + 1` magnitude bits. -This also holds for `U = 0`; no real logarithms or asymptotic claim are involved. -/ -theorem abs_lt_two_pow_log {U : ℕ} (x : ℤ) (hx : |x| ≤ U) : - |x| < 2 ^ (Nat.log 2 U + 1) := by - have h := Nat.lt_pow_succ_log_self (by decide : 1 < 2) U - have hcast : (U : ℤ) < 2 ^ (Nat.log 2 U + 1) := by exact_mod_cast h - exact hx.trans_lt hcast - -end Cslib.IntMatrix +end Cslib.Matrix diff --git a/Complexity/Foundations/Data/Nat/Choose/Real.lean b/Complexity/Foundations/Data/Nat/Choose/Real.lean index 8265f9eab3..b8e5e8f6aa 100644 --- a/Complexity/Foundations/Data/Nat/Choose/Real.lean +++ b/Complexity/Foundations/Data/Nat/Choose/Real.lean @@ -1,7 +1,7 @@ /- -Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Copyright (c) 2026 Samuel Schlesinger, Anthropic, PBC, Erik Demaine. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. -Authors: Samuel Schlesinger, OpenAI +Authors: Samuel Schlesinger, OpenAI, Anthropic, PBC, Erik Demaine -/ module @@ -34,6 +34,7 @@ generating-function bound, motivated by the geometric thresholds in Appendix A.3 * `Cslib.sum_choose_le_exp_mul_pow`: a real-parameter generating-function bound. * `Cslib.sum_choose_le_exp_one_mul_div_pow`: the bound with parameter `q / e`. +* `Nat.choose_le_exp_mul_div_pow`: its single-binomial-coefficient consequence, for all `n, k`. * `Cslib.sum_choose_le_two_rpow_div`: the geometric tail under an exact integer budget test. * `Cslib.sum_choose_le_two_pow_div_add_one`: its natural-number rounded form. * `Cslib.le_mul_logb_of_two_pow_le_sum_choose_sq`: the logarithmic consequence. @@ -113,6 +114,33 @@ theorem sum_choose_le_exp_one_mul_div_pow {q e : ℕ} (heq : e ≤ q) : rw [mul_div_assoc, mul_pow, ← Real.exp_nat_mul] simp +end Cslib + +namespace Nat + +/-- A single binomial coefficient is at most `(exp 1 * n / k)^k`, including `k = 0`. +This follows from `Cslib.sum_choose_le_exp_one_mul_div_pow` when `k ≤ n`; otherwise the +coefficient vanishes. The statement is also present in +[Anthropic's formal source](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Util/Choose.lean#L52-L53); +this proof reuses the partial-sum bound. -/ +theorem choose_le_exp_mul_div_pow (n k : ℕ) : + (n.choose k : ℝ) ≤ (Real.exp 1 * (n : ℝ) / (k : ℝ)) ^ k := by + by_cases hkn : k ≤ n + · calc + (n.choose k : ℝ) ≤ ∑ j ∈ Finset.range (k + 1), (n.choose j : ℝ) := + Finset.single_le_sum (fun j _ => Nat.cast_nonneg (n.choose j)) + (Finset.mem_range.mpr (by omega)) + _ ≤ (Real.exp 1 * (n : ℝ) / (k : ℝ)) ^ k := + Cslib.sum_choose_le_exp_one_mul_div_pow hkn + · rw [choose_eq_zero_of_lt (by omega), cast_zero] + positivity + +end Nat + +namespace Cslib + +open Finset + /-- For any input size `N`, the tail with dimension `C * N` and radius `Γ * N / M` is at most `2 ^ (N / B)` when the exact integer test `B * Γ * ⌈log₂ (4 * C * M)⌉ ≤ M` holds. This is an original corollary of `sum_choose_le_exp_mul_pow`, motivated by Appendix A.3 of diff --git a/Complexity/Foundations/NumberTheory/SqrtPrimeWindow.lean b/Complexity/Foundations/NumberTheory/SqrtPrimeWindow.lean new file mode 100644 index 0000000000..7c00a95c8c --- /dev/null +++ b/Complexity/Foundations/NumberTheory/SqrtPrimeWindow.lean @@ -0,0 +1,120 @@ +/- +Copyright (c) 2026 Anthropic, PBC, Erik Demaine. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Anthropic, PBC, Erik Demaine +-/ + +module +-- upstream: candidate + +public import Complexity.Foundations.NumberTheory.PrimeWindow + +/-! +# Prime windows below a square root + +This finite window is independent of the application that chooses a prime from it. Membership +has both real-endpoint and integer-square characterizations, including the closed lower and +open upper endpoints. Counting reuses the general prime-window estimates. + +## Main definitions + +* `Nat.sqrtPrimeWindow`: the primes in `[sqrt D / 2, sqrt D)`. + +## Main results + +* `Real.four_le_sqrt_natCast`: the shared endpoint bound for natural parameters at least 16. +* `Nat.mem_sqrtPrimeWindow`, `Nat.mem_sqrtPrimeWindow_iff_sq`: membership tests. +* `Nat.sqrtPrimeWindow_eq_filter`: an integer-only enumeration. +* `Nat.card_sqrtPrimeWindow_lt_sqrt`: the upper cardinality bound. +* `Nat.exists_pos_mul_sqrt_div_log_le_card_sqrtPrimeWindow`: a uniform lower bound. + +## References + +Adapted from the prime-window definitions and estimates in +[Anthropic's pinned hashing source](https://github.com/anthropics/formal-math/blob/16a5ce41304fcdaea4e0be7c132e1c361dc26da5/3sum-apsp/ThreeSumApsp/Sec3/Theorem17/Hashing.lean). +-/ + +@[expose] public section + +/-- The real square root of a natural window parameter at least 16 is at least four. -/ +theorem Real.four_le_sqrt_natCast {D : ℕ} (hD : 16 ≤ D) : 4 ≤ Real.sqrt (D : ℝ) := + Real.le_sqrt_of_sq_le (by norm_num; exact_mod_cast hD) + +namespace Nat + +open Classical in +/-- Prime numbers in the half-open real interval `[sqrt D / 2, sqrt D)`. +-/ +noncomputable def sqrtPrimeWindow (D : ℕ) : Finset ℕ := + (Finset.range D).filter fun p => p.Prime ∧ Real.sqrt D / 2 ≤ (p : ℝ) ∧ (p : ℝ) < Real.sqrt D + +/-- Window membership is equivalent to primality and the two real endpoint inequalities. The +finite enumeration bound in `sqrtPrimeWindow` excludes no member of the interval. +-/ +theorem mem_sqrtPrimeWindow {D p : ℕ} : + p ∈ sqrtPrimeWindow D ↔ p.Prime ∧ Real.sqrt D / 2 ≤ (p : ℝ) ∧ (p : ℝ) < Real.sqrt D := by + unfold sqrtPrimeWindow + 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_sqrtPrimeWindow_iff_sq {D p : ℕ} : + p ∈ sqrtPrimeWindow D ↔ p.Prime ∧ D ≤ 4 * p ^ 2 ∧ p ^ 2 < D := by + rw [mem_sqrtPrimeWindow] + 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 sqrtPrimeWindow_eq_filter (D : ℕ) : + sqrtPrimeWindow D = (Finset.range D).filter fun p => p.Prime ∧ D ≤ 4 * p ^ 2 ∧ p ^ 2 < D := by + ext p + rw [mem_sqrtPrimeWindow_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⟩ + +/-- For `D ≥ 16`, the number of primes in `sqrtPrimeWindow D` is strictly less than `sqrt D`. +-/ +theorem card_sqrtPrimeWindow_lt_sqrt (D : ℕ) (hD : 16 ≤ D) : + ((sqrtPrimeWindow D).card : ℝ) < Real.sqrt D := by + have hpos : 0 < Real.sqrt D := Real.sqrt_pos.mpr (by exact_mod_cast (by omega : 0 < D)) + -- The primes in the range are among the numbers `1, …, ⌈√D⌉ - 1`. + have hsub : sqrtPrimeWindow D ⊆ Finset.Ico 1 ⌈Real.sqrt D⌉₊ := by + intro p hp + obtain ⟨hprime, -, hlt⟩ := mem_sqrtPrimeWindow.1 hp + exact Finset.mem_Ico.2 ⟨hprime.one_lt.le, Nat.lt_ceil.2 hlt⟩ + have hcard : (sqrtPrimeWindow 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 : ((sqrtPrimeWindow D).card : ℝ) + 1 ≤ (⌈Real.sqrt D⌉₊ : ℝ) := by exact_mod_cast hcard + linarith [Nat.ceil_lt_add_one hpos.le] + +/-- 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_pos_mul_sqrt_div_log_le_card_sqrtPrimeWindow : + ∃ c : ℝ, 0 < c ∧ ∀ D : ℕ, 16 ≤ D → + c * (Real.sqrt D / Real.log D) ≤ ((sqrtPrimeWindow 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_sqrtPrimeWindow.2 ⟨hp, hge, hlt⟩⟩ + +end Nat diff --git a/ComplexityTests/ExactTriangleHashing.lean b/ComplexityTests/ExactTriangleHashing.lean index 41df09a874..16acad7bb4 100644 --- a/ComplexityTests/ExactTriangleHashing.lean +++ b/ComplexityTests/ExactTriangleHashing.lean @@ -20,17 +20,17 @@ 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] +private theorem window16 : Nat.sqrtPrimeWindow 16 = {2, 3} := by + rw [Nat.sqrtPrimeWindow_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] +example : 2 ∈ Nat.sqrtPrimeWindow 16 ∧ 5 ∉ Nat.sqrtPrimeWindow 25 ∧ 3 ∈ Nat.sqrtPrimeWindow 25 := by + simp only [Nat.mem_sqrtPrimeWindow_iff_sq] decide -example : 2 ∉ primesInRange 17 := by rw [mem_primesInRange_iff_sq]; decide +example : 2 ∉ Nat.sqrtPrimeWindow 17 := by rw [Nat.mem_sqrtPrimeWindow_iff_sq]; decide example : countZeroMod ab z z 2 = 8 ∧ countZeroMod ab z z 3 = 4 := by unfold countZeroMod weightSum ab z @@ -104,8 +104,8 @@ example : ¬WeightsBoundedBy (fun _ _ : Fin 1 => (0 : ℤ)) (fun _ _ => 0) (fun -- 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) + falsePositiveConst * (1 * (16 : ℝ) ^ 3 * Real.log 16 / Real.sqrt 16) := by + apply 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, ?_⟩ diff --git a/ComplexityTests/ExactTriangleSplitting.lean b/ComplexityTests/ExactTriangleSplitting.lean index 1768d67f64..3125f48121 100644 --- a/ComplexityTests/ExactTriangleSplitting.lean +++ b/ComplexityTests/ExactTriangleSplitting.lean @@ -43,8 +43,8 @@ private def zeroWeight : Fin 16 → Fin 16 → ℤ := fun _ _ => 0 private def falseWeight : Fin 16 → Fin 16 → ℤ := fun _ _ => 6 -- The prime window admits p = 3 at the smallest supported D. -private theorem prime3 : 3 ∈ primesInRange 16 := by - rw [mem_primesInRange_iff_sq] +private theorem prime3 : 3 ∈ Nat.sqrtPrimeWindow 16 := by + rw [Nat.mem_sqrtPrimeWindow_iff_sq] decide -- All-zero matrices supply an admitted, nonempty splitting with the promised bounds. @@ -139,7 +139,7 @@ example : ∃ L, IsScanOrder zeroWeight zeroWeight zeroWeight 16 2 3 private theorem falseSelected : IsSelectedPrime falseWeight zeroWeight zeroWeight 16 3 := by refine ⟨prime3, fun q hq => ?_⟩ - obtain ⟨hprime, -, hsq⟩ := mem_primesInRange_iff_sq.mp hq + obtain ⟨hprime, -, hsq⟩ := Nat.mem_sqrtPrimeWindow_iff_sq.mp hq have hq : q = 2 ∨ q = 3 := by have := hprime.two_le have : q < 4 := by nlinarith diff --git a/ComplexityTests/Lopsided.lean b/ComplexityTests/Lopsided.lean index 737cbf386a..ecae7913f4 100644 --- a/ComplexityTests/Lopsided.lean +++ b/ComplexityTests/Lopsided.lean @@ -5,6 +5,7 @@ Authors: Erik Demaine -/ import Complexity.Algorithms.ExactTriangle.Lopsided +import Complexity.Foundations.Data.Int.LogBounds import Mathlib.LinearAlgebra.Matrix.Notation /-! Concrete tests for padding, counts, queried answers and signed reconstruction. -/ @@ -80,7 +81,7 @@ example : (signedBitInstance signedX signedY {(0, 0)} true false 0 1).commonNeig classical simp only [signedBitInstance, commonNeighborCount_ofMatrices, ↓reduceIte] norm_num [signedX, signedY, - Cslib.IntMatrix.bit, posPart_def, negPart_def, Finset.univ_fin2, + Cslib.Matrix.bit, posPart_def, negPart_def, Finset.univ_fin2, Finset.filter_insert, Finset.filter_singleton] example : (signedBitInstance signedX signedY {(0, 0)} false true 1 0).commonNeighborCount @@ -88,7 +89,7 @@ example : (signedBitInstance signedX signedY {(0, 0)} false true 1 0).commonNeig classical simp only [signedBitInstance, commonNeighborCount_ofMatrices, ↓reduceIte] norm_num [signedX, signedY, - Cslib.IntMatrix.bit, posPart_def, negPart_def, Finset.univ_fin2, + Cslib.Matrix.bit, posPart_def, negPart_def, Finset.univ_fin2, Finset.filter_insert, Finset.filter_singleton] -- Apply the signed reconstruction theorem to an actual product, including all four signs. @@ -98,26 +99,27 @@ example : (∑ i ∈ Finset.range 2, ∑ j ∈ Finset.range 2, (2 : ℤ) ^ (i + - ((signedBitInstance signedX signedY {(0, 0)} true false i j).commonNeighborCount 0 0 : ℤ) + ((signedBitInstance signedX signedY {(0, 0)} true true i j).commonNeighborCount 0 0 : ℤ))) = -8 := by - rw [← entries_from_counts signedX signedY {(0, 0)} 2 + rw [← mul_apply_eq_sum_signedBitInstance_counts signedX signedY {(0, 0)} 2 (by intro a j; fin_cases a; fin_cases j <;> norm_num [signedX]) (by intro j b; fin_cases j <;> fin_cases b <;> norm_num [signedY])] norm_num [signedX, signedY, Matrix.mul_apply, Fin.sum_univ_two] -example : Cslib.IntMatrix.bit (!![-1, 0, 3] : Matrix (Fin 1) (Fin 3) ℤ) 1 0 0 = 1 := by - norm_num [Cslib.IntMatrix.bit] +example : Cslib.Matrix.bit (!![-1, 0, 3] : Matrix (Fin 1) (Fin 3) ℤ) 1 0 0 = 1 := by + norm_num [Cslib.Matrix.bit] example : |(-7 : ℤ)| < 2 ^ (Nat.log 2 7 + 1) := by - exact Cslib.IntMatrix.abs_lt_two_pow_log (-7) (by norm_num) + exact Int.abs_lt_two_pow_log_succ (-7) (by norm_num) example : |(0 : ℤ)| < 2 ^ (Nat.log 2 0 + 1) := by - exact Cslib.IntMatrix.abs_lt_two_pow_log 0 (by norm_num) + exact Int.abs_lt_two_pow_log_succ 0 (by norm_num) example : (signedX * (-signedY)) 0 0 = 8 := by norm_num [signedX, signedY, Matrix.mul_apply, Fin.sum_univ_two] example : (ofMatrices (!![1, 1] : Matrix (Fin 1) (Fin 2) ℤ) !![1; 1] {(0, 0)}).commonNeighborCount 0 0 = 2 := by - have h := zero_one (!![1, 1] : Matrix (Fin 1) (Fin 2) ℤ) !![1; 1] {(0, 0)} + have h := mul_apply_eq_commonNeighborCount_ofMatrices + (!![1, 1] : Matrix (Fin 1) (Fin 2) ℤ) !![1; 1] {(0, 0)} (by intro a j; fin_cases a; fin_cases j <;> simp) (by intro j b; fin_cases j <;> fin_cases b <;> simp) 0 0 norm_num [Matrix.mul_apply, Fin.sum_univ_two] at h @@ -128,7 +130,7 @@ example : (signedBitInstance signedX (-signedY) {(0, 0)} true true 0 1).commonNe classical simp only [signedBitInstance, commonNeighborCount_ofMatrices, ↓reduceIte] norm_num [signedX, signedY, - Cslib.IntMatrix.bit, posPart_def, negPart_def, Finset.univ_fin2, + Cslib.Matrix.bit, posPart_def, negPart_def, Finset.univ_fin2, Finset.filter_insert, Finset.filter_singleton] example : (signedBitInstances signedX signedY {(0, 0)} 2 (true, false, 0, 1)).W = @@ -149,7 +151,7 @@ example : Set.EqOn (fun _ : Fin 0 × Fin 0 => (99 : ℕ)) -- With beta zero, strict magnitude bounds are satisfied by zero matrices of nonempty shape. example : ((0 : Matrix (Fin 1) (Fin 2) ℤ) * (0 : Matrix (Fin 2) (Fin 1) ℤ)) 0 0 = 0 := by - have h := entries_from_counts (0 : Matrix (Fin 1) (Fin 2) ℤ) 0 {(0, 0)} 0 + have h := mul_apply_eq_sum_signedBitInstance_counts (0 : Matrix (Fin 1) (Fin 2) ℤ) 0 {(0, 0)} 0 (by simp) (by simp) 0 0 simpa only [Finset.range_zero, Finset.sum_empty] using h diff --git a/ComplexityTests/SchonhageThinProduct.lean b/ComplexityTests/SchonhageThinProduct.lean new file mode 100644 index 0000000000..9ab0067ef8 --- /dev/null +++ b/ComplexityTests/SchonhageThinProduct.lean @@ -0,0 +1,159 @@ +/- +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.MatrixMultiplication.Schonhage.Magnitude +import Complexity.Algorithms.MatrixMultiplication.Schonhage.Full +import Complexity.Algorithms.MatrixMultiplication.Schonhage.Pruned + +/-! +Regression clients for finite thin multiplication. Small computable layouts exercise actual signed +products and padding. The positive-depth parameter instance uses symbolic matrices so its large +outer dimension is never enumerated. No example claims a native instruction bound. +-/ + +open Finset Cslib.Algorithms.MatrixMultiplication.Schonhage + +namespace SchonhageThinProductTests + +/-- One outer level: blocks have three rows, and the inner dimension is one. -/ +private def outerLayout : Layout 1 0 where + hmL := by decide + rowIdx := (Equiv.funUnique (Fin 1) (Fin 3)).symm + colIdx := (Equiv.funUnique (Fin 1) (Fin 3)).symm + innerIdx := by + change Fin 1 ≃ (Fin 0 → Fin 2 × Fin 2) + exact Equiv.ofUnique _ _ + table := fun _ => ∅ + table_card := fun _ => rfl + table_injective := by + change Function.Injective (fun _ : Fin 1 × Fin 1 => _) + exact fun _ _ _ => Subsingleton.elim _ _ + +/-- One inner level: four inner indices and one row per block. -/ +private def innerLayout : Layout 1 1 where + hmL := by decide + rowIdx := by + change Fin 1 ≃ (Fin 0 → Fin 3) + exact Equiv.ofUnique _ _ + colIdx := by + change Fin 1 ≃ (Fin 0 → Fin 3) + exact Equiv.ofUnique _ _ + innerIdx := (finProdFinEquiv (m := 2) (n := 2)).symm.trans + (Equiv.funUnique (Fin 1) (Fin 2 × Fin 2)).symm + table := fun _ => {0} + table_card := fun _ => by decide + table_injective := by + change Function.Injective (fun _ : Fin 1 × Fin 1 => _) + exact fun _ _ _ => Subsingleton.elim _ _ + +private def X : Matrix (Fin 4) (Fin 1) ℤ := !![-3; 2; 0; 5] +private def Y : Matrix (Fin 1) (Fin 4) ℤ := !![2, -1, 4, 0] +private def W : Finset (Fin 4 × Fin 4) := {(0, 1), (3, 2), (2, 3)} + +-- N is not a multiple of the band size. The last two padded rows and columns are zero. +example : padN 1 0 4 = 6 := by decide +example : numBands 1 0 4 = 2 := by decide +example : padN 1 0 0 = 0 := by decide +kernel +example : padRows X 4 0 = 0 ∧ padCols Y 0 5 = 0 := by decide +kernel +example : padRows X 3 0 = 5 ∧ padCols Y 0 2 = 4 := by decide + +example : outStrOfPos outerLayout 0 1 = ![.z 0 1] := by decide +example : bandArrayL outerLayout X 0 ![.x 0] = -3 := by decide +kernel +-- These compute the finite recursion itself, including encodings and the wanted-output mapping. +example : Pruned.result 1 (encodingL (bandArrayL outerLayout X 0)) + (encodingR (bandArrayR outerLayout Y 0)) (wantedStrings outerLayout W 0 0) + (outStrOfPos outerLayout 0 1) = 3 := by decide +kernel + +example : Pruned.result 1 (encodingL (bandArrayL outerLayout X 1)) + (encodingR (bandArrayR outerLayout Y 0)) (wantedStrings outerLayout W 1 0) + (outStrOfPos outerLayout 3 2) = 20 := by decide +kernel + +example : Pruned.result 1 (encodingL (bandArrayL outerLayout X 0)) + (encodingR (bandArrayR outerLayout Y 1)) (wantedStrings outerLayout W 0 1) + (outStrOfPos outerLayout 2 3) = 0 := by decide +kernel + +example : ∑ T ∈ tiles 1 0 4, (wantedStrings outerLayout W T.1 T.2).card = 3 := by decide +kernel +example : ∑ T ∈ tiles 1 0 4, (leaves (wantedStrings outerLayout W T.1 T.2)).card = 3 := by + decide +example : wantedStrings outerLayout W 1 1 = ∅ := by decide +example : (leaves (wantedStrings outerLayout W 1 1)).card = 0 := by decide +kernel + +-- Four signed summands with a single inner output; there is no outer error coefficient here. +private def A : Matrix (Fin 1) (Fin 4) ℤ := !![-3, 2, 0, 5] +private def B : Matrix (Fin 4) (Fin 1) ℤ := !![2; -1; 4; 0] + +example : Full.result 1 (bandArrayL innerLayout A 0) (bandArrayR innerLayout B 0) + (outStrOfPos innerLayout 0 0) = -8 := by decide +kernel + +-- A mixed inner/outer output checks the exact-inner-count zeroing that eliminates error terms. +private def mixedX (_Q : Finset (Fin 2)) : LeftMat 2 1 := + fun _ π => A 0 (finProdFinEquiv (π 0)) +private def mixedY (_Q : Finset (Fin 2)) : RightMat 2 1 := + fun π _ => B (finProdFinEquiv (π 0)) 0 +private def mixedOutput : OutStr 2 := ![.z0, .z 0 0] + +example : Pruned.result 2 (encodingL (arrayL 1 mixedX)) (encodingR (arrayR 1 mixedY)) + {mixedOutput} mixedOutput = -8 := by decide +kernel + +example : Pruned.result 1 (encodingL (bandArrayL innerLayout A 0)) + (encodingR (bandArrayR innerLayout B 0)) + (wantedStrings innerLayout ({(0, 0)} : Finset (Fin 1 × Fin 1)) 0 0) + (outStrOfPos innerLayout 0 0) = -8 := by decide +kernel + +example : (univ.filter fun τ : Leaf 1 => + Leaf.Contributes τ (fun _ => OutVar.z0) ∧ order 1 τ = 0).card = 1 := by decide +example : (univ.filter fun τ : Leaf 1 => + Leaf.Contributes τ (fun _ => OutVar.z0) ∧ order 1 τ = 1).card = 9 := by decide +example : (leaves ({fun _ => OutVar.z0} : Finset (OutStr 1))).card = 10 := by decide +example : order 0 (fun _ : Fin 1 => Term.P0) = -1 := by decide +example : (leaves (∅ : Finset (OutStr 0))).card = 0 := by decide +kernel + +example (β β' : ℕ) : wantedStrings outerLayout (∅ : Finset (Fin 4 × Fin 4)) β β' = ∅ := by + simp [wantedStrings] + +-- Existence at both the degenerate and positive-depth boundary, without assumptions on layouts. +example : Nonempty (Layout 0 0) := nonempty_layout (by decide) +example : Nonempty (Layout 19 1) := nonempty_layout (by decide) + +private abbrev admittedN : ℕ := 68719476736 +private noncomputable def admittedLayout : Layout 19 1 := + Classical.choice (nonempty_layout (by decide)) +private def singletonQuery : Finset (Fin admittedN × Fin admittedN) := {(0, 0)} + +-- All parameters of the complete composition are simultaneously satisfiable with m > 0. +example : D 1 ^ 18 ≤ admittedN ∧ singletonQuery.card * 2 ^ 1 ≤ admittedN ^ 2 := by + constructor + · rfl + · norm_num [singletonQuery, admittedN] + +example : + ((2 * numBands (19 * 1) 1 admittedN * 10 ^ (19 * 1) : ℕ) : ℝ) ≤ + 4 * (D 1 : ℝ) ^ (-(1 / 18 : ℝ)) * (admittedN : ℝ) ^ 2 := by + have hN : D 1 ^ 18 ≤ admittedN := by norm_num [D, admittedN] + exact encodings_total_le 1 admittedN hN + +example : + ((∑ T ∈ tiles (19 * 1) 1 admittedN, + (leaves (wantedStrings admittedLayout singletonQuery T.1 T.2)).card : ℕ) : ℝ) ≤ + 33 * (D 1 : ℝ) ^ (-(1 / 18 : ℝ)) * (admittedN : ℝ) ^ 2 := by + have hN : D 1 ^ 18 ≤ admittedN := by norm_num [D, admittedN] + have hW : singletonQuery.card * 2 ^ 1 ≤ admittedN ^ 2 := by + norm_num [singletonQuery, admittedN] + exact total_leaves_le (m := 1) (N := admittedN) admittedLayout singletonQuery hN hW + +example (X' : Matrix (Fin admittedN) (Fin (D 1)) ℤ) + (Y' : Matrix (Fin (D 1)) (Fin admittedN) ℤ) : + ((∑ T ∈ tiles (19 * 1) 1 admittedN, + Pruned.totalSize (encodingL (bandArrayL admittedLayout X' T.1)) + (encodingR (bandArrayR admittedLayout Y' T.2)) + (wantedStrings admittedLayout singletonQuery T.1 T.2) : ℕ) : ℝ) ≤ + ((19 * 1 + 1 : ℕ) : ℝ) * (33 * (D 1 : ℝ) ^ (-(1 / 18 : ℝ)) * (admittedN : ℝ) ^ 2) := by + have hN : D 1 ^ 18 ≤ admittedN := by norm_num [D, admittedN] + have hW : singletonQuery.card * 2 ^ 1 ≤ admittedN ^ 2 := by + norm_num [singletonQuery, admittedN] + exact total_passed_sets_le (m := 1) (N := admittedN) admittedLayout X' Y' singletonQuery hN hW + +end SchonhageThinProductTests diff --git a/NOTICE b/NOTICE index 5fb9b6c145..dcdc12ee98 100644 --- a/NOTICE +++ b/NOTICE @@ -3,15 +3,21 @@ Complexity/Algorithms/ExactTriangle/{Defs,Weights,Lopsided,Hashing,Splitting,Wit Complexity/Algorithms/MatrixMultiplication/Schonhage/, Complexity/Algorithms/ThreeSum/, Complexity/Algorithms/WordRAM/ThreeSum/, -Complexity/Foundations/Data/Finset/{Average,Quotient}.lean, +Complexity/Foundations/Algebra/{FiniteMagnitude,GeometricTail}.lean, +Complexity/Foundations/Combinatorics/FunctionPlaces.lean, +Complexity/Foundations/Data/Finset/{Average,OrderGlue,Quotient}.lean, Complexity/Foundations/Data/Int/Residue.lean, Complexity/Foundations/Data/Nat/CeilDiv.lean, Complexity/Foundations/Data/Matrix/{IntegerBits,PairChunks}.lean and -Complexity/Foundations/NumberTheory/{PrimeCounting,PrimeWindow}.lean and +Complexity/Foundations/NumberTheory/{PrimeCounting,PrimeWindow,SqrtPrimeWindow}.lean and Complexity/Foundations/RingTheory/CyclicIntRing.lean are adapted from anthropics/formal-math, 3sum-apsp, revision 16a5ce41304fcdaea4e0be7c132e1c361dc26da5. +The single-coefficient binomial bound in Complexity/Foundations/Data/Nat/Choose/Real.lean +retains its source credit from that adaptation; its proof now derives the statement from the +library's existing partial-binomial-sum bound. + ThreeSumApsp: a Lean 4 formalization of "Truly Subquadratic 3SUM and Truly Subcubic APSP via Triangles in Sparse Lopsided Graphs" Copyright 2026 Anthropic, PBC diff --git a/ROADMAP.md b/ROADMAP.md index ddae6f5ff5..edfbd97314 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -429,8 +429,10 @@ each statement against the numbered result. Organize modules by subject, not cha and the selected-prime false-positive bound ([#1442]), finite residue/chunk/piece splitting and witness recovery ([#1447]) are proved mathematically. The tensor identity, full recursion, shared encodings and pruned-evaluation core ([#1459]) include exact visited-leaf and passed-set - bounds, with the empty zero-depth root handled explicitly. It still needs tiling and sparse - leaf bounds, native tensor/sparse-query programs, native Chan–He preprocessing, the reduction + bounds, with the empty zero-depth root handled explicitly. The finite thin-product component + ([#1463]) proves valid layouts, arbitrary-size padding/tiling, exact queried signed products and + sparse encoding/leaf/passed-set totals. These counts are not native instruction costs. It still + needs native 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 diff --git a/UPSTREAM.md b/UPSTREAM.md index 153abe289b..1d86bbb455 100644 --- a/UPSTREAM.md +++ b/UPSTREAM.md @@ -177,14 +177,18 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he common-neighbor counts as product entries, detection, the zero/one converse, and signed products from `4 * β²` count instances with the original queried edges. Its digit identities live in `Foundations/Data/Matrix/IntegerBits.lean`, reuse Mathlib's pointwise positive/negative parts, - and give an explicit `Nat.log 2 U + 1` magnitude-bit bound including zero. Query answers use - `Set.EqOn`; no additional weighted-triangle predicate or source machine is introduced. + with descriptive lemmas in `Cslib.Matrix`. The scalar magnitude-bit bound, including zero, + lives in `Foundations/Data/Int/LogBounds.lean`. Count answers use `Set.EqOn`; detection + consumes a direct iff on queried edges. 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. `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. + predicate. The application-independent `Nat.sqrtPrimeWindow` API lives in + `Foundations/NumberTheory/SqrtPrimeWindow.lean`; its integer-square characterization supports + native enumeration, and splitting and witness recovery reuse it. 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 @@ -199,10 +203,15 @@ Do not add a ledger row for an ordinary candidate or repeat its source marker he error, full evaluation, shared encodings, pruned values, visited leaves and passed-set totals. Its pure finite recursion traces are not native machine semantics or charged instruction costs. The zero-depth empty-query root exception is explicit, with an unconditional root-or-leaf API. - See [#1459](https://github.com/cslib-community/complexitylib/issues/1459); tiling/sparse leaf bounds - remain in #1394. - The existing - signed and unsigned APIs were reused; no source machine or separate run semantics were copied. + See [#1459](https://github.com/cslib-community/complexitylib/issues/1459). The complete finite + thin-product component ([#1463](https://github.com/cslib-community/complexitylib/issues/1463)) + adapts `Sec2/{Levels,Lemma9,Tiling,Orders,Lemma11,Theorem5}` and needed finite magnitude bounds. + Layouts contain only indexing choices, whose existence is proved; sparse bounds are derived + from exact contributing-leaf counts. Generic function-position counts, ordered gluing, binomial + and finite-sum bounds live in small Foundations modules. Mathlib `Finset.orderEmbOfFin` indexes + both a level set and its complement, and `Finset.piecewise` restricts answers. No source machine, + problem-statement universe or solver/time hypothesis is imported. + The existing signed and unsigned APIs were reused; no source machine or separate run semantics were copied. Tensor/sparse query programs, the full reduction host, native Chan–He preprocessing, fast Exact Triangle and full input/time/width/asymptotic composition remain in [#1394](https://github.com/cslib-community/complexitylib/issues/1394), including the improved @@ -286,6 +295,7 @@ agreement during upstream review. | `Foundations/Combinatorics/{FiniteCover,Hoeffding,Hypergeometric,PowersetTail,Sunflower}` | Consider Mathlib combinatorics; probabilistic-method lemmas belong near `Combinatorics/Probabilistic`. | | `Foundations/Data/Finset/Expect` | Product-type expectation helpers extracted from `Foundations/Combinatorics/Expander/Products` (#317); Mathlib candidates. | | `Foundations/Data/Finset/Markov` | Mathlib candidate: `Algebra/Order/BigOperators/Group/Finset` (`Finset.card_filter_mul_le_sum`, `Finset.mul_card_filter_le_sum`). | +| `Foundations/Data/Finset/OrderGlue` | Follow-up to #1496: rename the generic `forall_level_iff`, `exists_innerLevel`, `outerLevel_notMem` and `filter_glue` using order-embedding/complement terminology (for example `forall_iff_orderEmb_and_compl`), retaining this topic module. Migrate the concrete `Schonhage/Levels` and `Schonhage/Batch` clients together with the active native band-preparation client; defer that API migration until the latter is stable. | | `Foundations/Data/Matrix/Permanent` | Mathlib candidate: Laplace expansion of the permanent along row and column 0, beside `LinearAlgebra/Matrix/Permanent` (mirroring `Matrix.det_succ_row_zero`). | | `Foundations/LinearAlgebra/{Rank,Matrix/CutRank,Matrix/Lindsey}` | Mathlib candidates: `LinearAlgebra/Dimension/LinearMap` and `LinearAlgebra/Matrix/Rank`; placement of Lindsey's lemma needs agreement. | | `Probability/PairwiseIndependentHash/Shift` | Consider integration into `PairwiseIndependentHash`. |