feat(wordram): build prime tables and collision minimum searches

Complexity.lean

@@ -118,6 +118,8 @@ public import Complexity.Algorithms.WordRAM.PrefixSums.Dynamic
 public import Complexity.Algorithms.WordRAM.PrefixSums.Query
 public import Complexity.Algorithms.WordRAM.PrefixSums.Tree
 public import Complexity.Algorithms.WordRAM.PrefixSums.Update
+public import Complexity.Algorithms.WordRAM.PrimeTable
+public import Complexity.Algorithms.WordRAM.PrimeTrial
 public import Complexity.Algorithms.WordRAM.PriorityQueue.Basic
 public import Complexity.Algorithms.WordRAM.PriorityQueue.BinaryHeap
 public import Complexity.Algorithms.WordRAM.PriorityQueue.BinaryHeap.Model
@@ -144,10 +146,16 @@ public import Complexity.Algorithms.WordRAM.RadixSort.Scan
 public import Complexity.Algorithms.WordRAM.Stack
 public import Complexity.Algorithms.WordRAM.Stack.Controller
 public import Complexity.Algorithms.WordRAM.ThreeSum.Cell
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.CollisionTest
 public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Collisions
 public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Count
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.CountCall
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.CountFrame
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.CountTriple
 public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Parameters
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Query
 public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Residues
+public import Complexity.Algorithms.WordRAM.ThreeSum.ChanHe.Search
 public import Complexity.Algorithms.WordRAM.ThreeSum.Fill
 public import Complexity.Algorithms.WordRAM.ThreeSum.Host.Loop
 public import Complexity.Algorithms.WordRAM.ThreeSum.Host.Prepare

Complexity/Algorithms/WordRAM/PrimeTable.lean

@@ -0,0 +1,224 @@
+/-
+Copyright (c) 2026 Erik Demaine. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Erik Demaine
+-/
+
+module
+-- upstream: candidate
+
+public import Complexity.Algorithms.WordRAM.PrimeTrial
+public import Mathlib.NumberTheory.PrimeCounting
+
+/-!
+# Native prime tables
+
+This fixed program tests every candidate through the supplied bound, writes its primality
+bit, and counts the primes. Trial division, table addresses, writes and loops are charged.
+There is no primality oracle or assumed sieve.
+
+## Main definitions
+
+* `PrimeTable.program`: the complete native table builder.
+* `PrimeTable.memory`: its output bitmap and memory frame.
+
+## Main results
+
+* `PrimeTable.exec`: exact bitmap and prime count in an explicit quadratic instruction bound.
+-/
+
+@[expose] public section
+
+namespace Cslib.Algorithms.WordRAM.PrimeTable
+
+open Cslib.WordRAM
+open Finset
+
+/-- Registers of the trial-division subroutine. -/
+def trialRegisters (i : Fin 6) : Fin 15 := ⟨i.val, by omega⟩
+
+private theorem trial_injective : Function.Injective trialRegisters := by
+  intro i j h
+  apply Fin.ext
+  exact congrArg (fun i : Fin 15 => i.val) h
+
+/-- Current trial-division candidate. -/
+abbrev arg : Fin 15 := 0
+/-- Current candidate's primality bit. -/
+abbrev flag : Fin 15 := 2
+/-- Largest candidate. -/
+abbrev limit : Fin 15 := 6
+/-- Output bitmap base. -/
+abbrev base : Fin 15 := 7
+/-- First candidate not yet written. -/
+abbrev cursor : Fin 15 := 8
+/-- Address scratch. -/
+abbrev addr : Fin 15 := 9
+/-- Number of primes written. -/
+abbrev count : Fin 15 := 10
+/-- One past the largest candidate. -/
+abbrev stop : Fin 15 := 11
+/-- Constant one. -/
+abbrev one : Fin 15 := 12
+/-- Constant zero. -/
+abbrev zero : Fin 15 := 13
+/-- Constant two. -/
+abbrev two : Fin 15 := 14
+
+/-- Handle zero and one directly; test all larger candidates by charged native trial division. -/
+def test : Code .standard 15 :=
+  .branch (.lt cursor two) (ram[ flag ← #0 ])
+    (.seq (ram[ arg ← cursor ]) (PrimeTrial.program.rename trialRegisters))
+
+/-- Test, write the bit, update the count, and advance the candidate. -/
+def body : Code .standard 15 :=
+  .seq test (ram[ addr ← base + cursor; *addr ← flag;
+    count ← count + flag; cursor ← cursor + one ])
+
+/-- A finite native prime-table program with its bound and base supplied in registers. -/
+def program : Code .standard 15 :=
+  .seq (ram[ one ← #1; zero ← #0; two ← #2; cursor ← #0; count ← #0;
+    stop ← limit + one ]) (.loop (.lt cursor stop) body)
+
+/-- The first `k` primality bits, preserving all other memory cells. -/
+def memory (mem : ℕ → ℕ) (base k : ℕ) : ℕ → ℕ :=
+  fun a => if base ≤ a ∧ a < base + k then if (a - base).Prime then 1 else 0 else mem a
+
+private theorem memory_zero (mem : ℕ → ℕ) (base : ℕ) : memory mem base 0 = mem := by
+  funext a
+  simp [memory]
+
+private theorem memory_succ (mem : ℕ → ℕ) (base k : ℕ) :
+    memory mem base (k + 1) = Function.update (memory mem base k)
+      (base + k) (if k.Prime then 1 else 0) := by
+  funext a
+  by_cases ha : a = base + k
+  · subst a
+    simp [memory]
+  · rw [Function.update_of_ne ha]
+    have h : (base ≤ a ∧ a < base + (k + 1)) ↔ (base ≤ a ∧ a < base + k) := by omega
+    simp only [memory, h]
+
+private theorem exec_test {w m : ℕ} (x : Data 15) (htwo : x.regs two = 2)
+    (hbound : x.regs cursor ≤ m) (hfit : Fits w (m + 3)) :
+    ∃ y, ExecWithin w test x y (7 * m + 9) ∧
+      y.regs flag = (if (x.regs cursor).Prime then 1 else 0) ∧ y.mem = x.mem ∧
+      ∀ i, 6 ≤ i.val → y.regs i = x.regs i := by
+  by_cases hp : x.regs cursor < 2
+  · have hnp : ¬(x.regs cursor).Prime := fun h => by have := h.two_le; omega
+    refine ⟨x.set flag 0, ?_, by simp [hnp], rfl, ?_⟩
+    · exact (ExecWithin.branch_pos (by simpa only [Cond.holds_lt, htwo] using hp)
+        (ExecWithin.const (fits_zero w))).mono (by omega)
+    · intro i hi
+      have hne : i ≠ flag := by intro h; subst i; change 6 ≤ 2 at hi; omega
+      simp only [Data.set_regs, Function.update_of_ne hne]
+  · let p := x.set arg (x.regs cursor)
+    have hcopy : ExecWithin (B := .standard) w (ram[ arg ← cursor ]) x p 1 :=
+      ExecWithin.copy
+    obtain ⟨a, ha, haf, haa, ham⟩ := PrimeTrial.exec (p.comap trialRegisters)
+      (by change 2 ≤ x.regs cursor; omega)
+      (hfit.of_le (by change x.regs cursor + 2 ≤ m + 3; omega))
+    let y := a.extend trialRegisters p
+    have hf : y.regs flag = a.regs PrimeTrial.flag := by
+      change Function.extend trialRegisters a.regs p.regs (trialRegisters PrimeTrial.flag) = _
+      exact trial_injective.extend_apply _ _ _
+    refine ⟨y, ?_, hf.trans haf, ham, ?_⟩
+    · exact (ExecWithin.branch_neg (by simpa only [Cond.holds_lt, htwo] using hp)
+        (hcopy.seq (ha.rename trial_injective))).mono (by
+          change 1 + (7 * (x.regs cursor - 2) + 6) + 2 ≤ 7 * m + 9
+          omega)
+    · intro i hi
+      change Function.extend trialRegisters a.regs p.regs i = x.regs i
+      rw [Function.extend_apply']
+      · have hne : i ≠ arg := by intro h; subst i; change 6 ≤ 0 at hi; omega
+        exact Function.update_of_ne hne _ _
+      · rintro ⟨j, hj⟩
+        have := congrArg Fin.val hj
+        dsimp [trialRegisters] at this
+        omega
+
+
+/-- Write every primality bit through `m`, count exactly `Nat.primesLE m`, and preserve
+all cells outside the bitmap. Trial division gives the explicit quadratic bound
+`(m+1)(7m+15)+8`; the final halt is excluded. -/
+theorem exec {w : ℕ} (x : Data 15)
+    (hfit : Fits w (x.regs limit + 3))
+    (haddr : Fits w (x.regs base + (x.regs limit + 1))) :
+    ∃ y, ExecWithin w program x y ((x.regs limit + 1) * (7 * x.regs limit + 15) + 8) ∧
+      y.mem = memory x.mem (x.regs base) (x.regs limit + 1) ∧
+      y.regs count = #(Nat.primesLE (x.regs limit)) ∧
+      y.regs limit = x.regs limit ∧ y.regs base = x.regs base := by
+  let m := x.regs limit
+  let I (y : Data 15) : Prop := y.regs limit = m ∧ y.regs base = x.regs base ∧
+    y.regs stop = m + 1 ∧ y.regs one = 1 ∧ y.regs zero = 0 ∧ y.regs two = 2 ∧
+    y.regs cursor ≤ m + 1 ∧ y.regs count = Nat.primeCounting' (y.regs cursor) ∧
+    y.mem = memory x.mem (x.regs base) (y.regs cursor)
+  have hbody (y : Data 15) (hy : I y) (hc : (Cond.lt cursor stop).Holds y) :
+      ∃ q, ExecWithin w body y q (7 * m + 13) ∧ I q ∧
+        m + 1 - q.regs cursor < m + 1 - y.regs cursor := by
+    obtain ⟨hl, hb, hs, ho, hz, ht, hi, hcount, hm⟩ := hy
+    have hlt : y.regs cursor < m + 1 := by simpa only [Cond.holds_lt, hs] using hc
+    obtain ⟨z, he, hzf, hzm, hzr⟩ := exec_test y ht (show y.regs cursor ≤ m by omega) hfit
+    have hzb : z.regs base = x.regs base := (hzr base (by decide)).trans hb
+    have hzi : z.regs cursor = y.regs cursor := hzr cursor (by decide)
+    have hzc : z.regs count = Nat.primeCounting' (y.regs cursor) :=
+      (hzr count (by decide)).trans hcount
+    have hzo : z.regs one = 1 := (hzr one (by decide)).trans ho
+    have hzf_le : z.regs flag ≤ 1 := by rw [hzf]; split_ifs <;> omega
+    have ha : Fits w (z.regs base + z.regs cursor) := haddr.of_le (by rw [hzb, hzi]; omega)
+    have hn : Fits w (z.regs cursor + 1) := hfit.of_le (by rw [hzi]; change _ ≤ m + 3; omega)
+    have hct : Fits w (z.regs count + z.regs flag) := hfit.of_le (by
+      rw [hzc]
+      have hle : Nat.primeCounting' (y.regs cursor) ≤ y.regs cursor :=
+        Nat.count_le Nat.Prime
+      change _ ≤ m + 3
+      omega)
+    refine ⟨?q, ?run, ?inv, ?decrease⟩
+    case run =>
+      have htail : ExecWithin w (ram[ addr ← base + cursor; *addr ← flag;
+          count ← count + flag; cursor ← cursor + one ]) z ?q 4 := by
+        ram_eval [ha, hn, hct, hzo]
+      exact (he.seq htail).mono (by omega)
+    case inv =>
+      refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩
+      · simpa +decide [Function.update_apply] using (hzr limit (by decide)).trans hl
+      · simpa +decide [Function.update_apply] using hzb
+      · simpa +decide [Function.update_apply] using (hzr stop (by decide)).trans hs
+      · simpa +decide [Function.update_apply] using hzo
+      · simpa +decide [Function.update_apply] using (hzr zero (by decide)).trans hz
+      · simpa +decide [Function.update_apply] using (hzr two (by decide)).trans ht
+      · simpa +decide [Function.update_apply, hzi] using Nat.succ_le_of_lt hlt
+      · change z.regs count + z.regs flag = Nat.primeCounting' (z.regs cursor + 1)
+        rw [hzc, hzf, hzi]
+        simp only [Nat.primeCounting', Nat.count_succ]
+      · change Function.update z.mem (z.regs base + z.regs cursor) (z.regs flag) =
+          memory x.mem (x.regs base) (z.regs cursor + 1)
+        rw [hzm, hm, hzb, hzi, hzf, ← memory_succ]
+    case decrease =>
+      change m + 1 - (z.regs cursor + 1) < m + 1 - y.regs cursor
+      rw [hzi]
+      omega
+  let x₀ := (((((x.set one 1).set zero 0).set two 2).set cursor 0).set count 0).set stop (m + 1)
+  have hx₀ : I x₀ := by
+    simp +decide [I, x₀, m, limit, base, stop, one, zero, two, cursor, count,
+      Nat.primeCounting', Nat.count_zero, memory_zero]
+  obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I
+    (fun y => m + 1 - y.regs cursor) hbody hx₀
+  obtain ⟨hl, hb, hs, ho, hz, ht, hi, hcount, hm⟩ := hy
+  have hend : y.regs cursor = m + 1 := by
+    simp only [Cond.holds_lt, hs] at hc
+    omega
+  have h1 : Fits w 1 := hfit.of_le (by omega)
+  have h2 : Fits w 2 := hfit.of_le (by omega)
+  have hstop : Fits w (x.regs limit + 1) := hfit.of_le (by omega)
+  have hinit : ExecWithin (B := .standard) w
+      (ram[ one ← #1; zero ← #0; two ← #2; cursor ← #0; count ← #0;
+        stop ← limit + one ]) x x₀ 6 := by
+    ram_eval [x₀, m, h1, h2, hstop, limit, base, stop, one, zero, two, cursor, count, Fin.mk.injEq]
+  refine ⟨y, ?_, by simpa only [hend] using hm, ?_, hl, hb⟩
+  · have hi₀ : x₀.regs cursor = 0 := rfl
+    have hcost : 7 * m + 13 + 2 = 7 * m + 15 := by omega
+    exact (hinit.seq he).mono (by rw [hi₀, Nat.sub_zero, hcost]; dsimp only [m]; omega)
+  · rw [hcount, hend, Nat.primesLE_card_eq_primeCounting, Nat.primeCounting]
+
+end Cslib.Algorithms.WordRAM.PrimeTable

Complexity/Algorithms/WordRAM/PrimeTrial.lean

@@ -0,0 +1,142 @@
+/-
+Copyright (c) 2026 Erik Demaine. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Erik Demaine
+-/
+
+module
+-- upstream: candidate
+
+public import Complexity.Computability.WordRAM.Notation
+public import Mathlib.Data.Nat.Prime.Defs
+
+/-!
+# Charged native trial division
+
+This is a finite instruction program, not an invocation of Lean's primality decision procedure.
+The Chan–He prime table will call it only on candidates at least two.
+
+## Main definitions
+
+* `PrimeTrial.program`: trial divisors from two up to the candidate.
+
+## Main results
+
+* `PrimeTrial.exec`: native trial division with an explicit linear instruction bound.
+-/
+
+@[expose] public section
+
+namespace Cslib.Algorithms.WordRAM.PrimeTrial
+
+open Cslib.WordRAM
+
+/-- Candidate, preserved. -/
+abbrev arg : Fin 6 := 0
+/-- Current trial divisor. -/
+abbrev divisor : Fin 6 := 1
+/-- Primality flag. -/
+abbrev flag : Fin 6 := 2
+/-- Remainder scratch. -/
+abbrev remainder : Fin 6 := 3
+/-- Constant one. -/
+abbrev one : Fin 6 := 4
+/-- Constant zero. -/
+abbrev zero : Fin 6 := 5
+
+/-- Test one trial divisor and advance. The conditional also charges its control instructions. -/
+def body : Code .standard 6 := ram[
+  remainder ← arg % divisor;
+  if remainder = zero then [ flag ← zero ] else [ flag ← flag ];
+  divisor ← divisor + one ]
+
+/-- Initialize and test every divisor strictly between one and the candidate. -/
+def program : Code .standard 6 :=
+  .seq (ram[ divisor ← #2; flag ← #1; one ← #1; zero ← #0 ])
+    (.loop (.lt divisor arg) body)
+
+private def Good (p k : ℕ) : Prop := ∀ d, 2 ≤ d → d < k → p % d ≠ 0
+
+private theorem good_succ (p k : ℕ) (hk : 2 ≤ k) :
+    Good p (k + 1) ↔ Good p k ∧ p % k ≠ 0 := by
+  constructor
+  · intro h
+    exact ⟨fun d hd hdk => h d hd (by omega), h k hk (by omega)⟩
+  · rintro ⟨h, hmod⟩ d hd hdk
+    by_cases hdk' : d = k
+    · simpa only [hdk'] using hmod
+    · exact h d hd (by omega)
+
+/-- Test a candidate at least two in at most `7(p-2)+6` native instructions.
+All divisors are actual register values; the native remainder instruction is charged. -/
+theorem exec {w : ℕ} (x : Data 6) (hp : 2 ≤ x.regs arg)
+    (hfit : Fits w (x.regs arg + 2)) :
+    ∃ y, ExecWithin w program x y (7 * (x.regs arg - 2) + 6) ∧
+      y.regs flag = (if (x.regs arg).Prime then 1 else 0) ∧
+      y.regs arg = x.regs arg ∧ y.mem = x.mem := by
+  classical
+  let p := x.regs arg
+  let I (y : Data 6) : Prop := y.regs arg = p ∧ y.regs one = 1 ∧ y.regs zero = 0 ∧
+    2 ≤ y.regs divisor ∧ y.regs divisor ≤ p ∧
+    y.regs flag = (if Good p (y.regs divisor) then 1 else 0) ∧ y.mem = x.mem
+  have hbody (y : Data 6) (hy : I y) (hc : (Cond.lt divisor arg).Holds y) :
+      ∃ z, ExecWithin w body y z 5 ∧ I z ∧ p - z.regs divisor < p - y.regs divisor := by
+    obtain ⟨ha, ho, hz, hd, hdp, hf, hm⟩ := hy
+    have hlt : y.regs divisor < p := by simpa only [Cond.holds_lt, ha] using hc
+    have hmod : Fits w (p % y.regs divisor) := hfit.of_le (by
+      have := Nat.mod_le p (y.regs divisor)
+      change p % y.regs divisor ≤ p + 2
+      omega)
+    have hinc : Fits w (y.regs divisor + 1) := hfit.of_le (by change _ ≤ p + 2; omega)
+    let z := ((y.set remainder (p % y.regs divisor)).set flag
+      (if p % y.regs divisor = 0 then 0 else y.regs flag)).set divisor (y.regs divisor + 1)
+    have he : ExecWithin w body y z 5 := by
+      by_cases hrem : p % y.regs divisor = 0 <;>
+        ram_eval [body, z, ha, ho, hz, hmod, hinc, hrem, arg, divisor, flag,
+          remainder, one, zero, Fin.mk.injEq]
+    have hflag : (if p % y.regs divisor = 0 then 0 else y.regs flag) =
+        if Good p (y.regs divisor + 1) then 1 else 0 := by
+      rw [hf]
+      have hg := good_succ p (y.regs divisor) hd
+      by_cases hrem : p % y.regs divisor = 0
+      · have hg' : ¬Good p (y.regs divisor + 1) := fun h => (hg.mp h).2 hrem
+        simp only [hrem, hg', ↓reduceIte]
+      · have hg' : Good p (y.regs divisor + 1) ↔ Good p (y.regs divisor) := by tauto
+        simp only [hrem, hg', ↓reduceIte]
+    refine ⟨z, he, ?_, ?_⟩
+    · refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_⟩
+      · simpa +decide [z, Function.update_apply, arg, divisor, flag, remainder] using ha
+      · simpa +decide [z, Function.update_apply, one, divisor, flag, remainder] using ho
+      · simpa +decide [z, Function.update_apply, zero, divisor, flag, remainder] using hz
+      · simpa +decide [z, Function.update_apply, divisor] using
+          (show 2 ≤ y.regs divisor + 1 by omega)
+      · simpa +decide [z, Function.update_apply, divisor] using Nat.succ_le_of_lt hlt
+      · simpa +decide [z, Function.update_apply, divisor, flag, remainder] using hflag
+      · exact hm
+    · change p - (y.regs divisor + 1) < p - y.regs divisor
+      omega
+  let x₀ := (((x.set divisor 2).set flag 1).set one 1).set zero 0
+  have hx₀ : I x₀ := by
+    have hgood : Good p 2 := by intro d hd hlt; omega
+    simp +decide [I, x₀, p, arg, divisor, flag, one, zero, hgood, hp]
+  obtain ⟨y, he, hy, hc⟩ := ExecWithin.loop_invariant I (fun y => p - y.regs divisor) hbody hx₀
+  obtain ⟨ha, ho, hz, hd, hdp, hf, hm⟩ := hy
+  have hend : y.regs divisor = p := by simp only [Cond.holds_lt, ha] at hc; omega
+  have hprime : Good p p ↔ p.Prime := by
+    rw [Nat.prime_def_lt']
+    constructor
+    · intro h
+      refine ⟨hp, fun d hd hdp hdiv => ?_⟩
+      exact h d hd hdp (Nat.mod_eq_zero_of_dvd hdiv)
+    · rintro ⟨_, h⟩ d hd hdp hzero
+      exact h d hd hdp (Nat.dvd_of_mod_eq_zero hzero)
+  have h1 : Fits w 1 := hfit.of_le (by omega)
+  have h2 : Fits w 2 := hfit.of_le (by omega)
+  have hinit : ExecWithin (B := .standard) w
+      (ram[ divisor ← #2; flag ← #1; one ← #1; zero ← #0 ]) x x₀ 4 := by
+    ram_eval [x₀, h1, h2]
+  refine ⟨y, ?_, by simpa only [hend, hprime] using hf, ha, hm⟩
+  have hd₀ : x₀.regs divisor = 2 := rfl
+  exact (hinit.seq he).mono (by rw [hd₀]; dsimp only [p]; omega)
+
+end Cslib.Algorithms.WordRAM.PrimeTrial