diff --git a/Complexity.lean b/Complexity.lean index f06264477..9088c1543 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -89,6 +89,9 @@ public import Complexity.Algorithms.ThreeSum.Convolution public import Complexity.Algorithms.ThreeSum.Defs public import Complexity.Algorithms.ThreeSum.ThreeSets public import Complexity.Algorithms.WordRAM.AKS.Congruence +public import Complexity.Algorithms.WordRAM.AKS.CongruenceLoop +public import Complexity.Algorithms.WordRAM.AKS.GcdTotient +public import Complexity.Algorithms.WordRAM.AKS.OrderSearch public import Complexity.Algorithms.WordRAM.AKS.PerfectPower public import Complexity.Algorithms.WordRAM.AKS.PolyArray public import Complexity.Algorithms.WordRAM.AKS.PolyMul @@ -2859,6 +2862,7 @@ public import Complexity.Foundations.MeasureTheory.Constructions.BorelSpace.Orde public import Complexity.Foundations.NumberTheory.AKS.Coefficients public import Complexity.Foundations.NumberTheory.AKS.Correctness public import Complexity.Foundations.NumberTheory.AKS.Criterion +public import Complexity.Foundations.NumberTheory.AKS.IntegerVariant public import Complexity.Foundations.NumberTheory.AKS.Introspective public import Complexity.Foundations.NumberTheory.AKS.PerfectPower public import Complexity.Foundations.NumberTheory.AKS.SmallR diff --git a/Complexity/Algorithms/WordRAM/AKS/CongruenceLoop.lean b/Complexity/Algorithms/WordRAM/AKS/CongruenceLoop.lean new file mode 100644 index 000000000..f1b0a8c86 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/AKS/CongruenceLoop.lean @@ -0,0 +1,202 @@ +/- +Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Samuel Schlesinger +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.AKS.Congruence + +/-! +# Step 5 of the AKS algorithm on the word RAM + +Step 5 of the AKS algorithm [Agrawal, Kayal, Saxena, Section 4][AKS2004] tests the congruence +`(X + a)^n ≡ X^n + a (mod X^r - 1, n)` for `a = 1, …, A`. The program runs the test for one `a` +(`Cslib.Algorithms.WordRAM.AKS.Congruence.program`) for each `a` in turn and clears a flag when +one fails. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.AKS.CongruenceLoop.program`: the loop. +* `Cslib.Algorithms.WordRAM.AKS.CongruenceLoop.time`: its instruction bound. + +## Main results + +* `Cslib.Algorithms.WordRAM.AKS.CongruenceLoop.exec`: the flag is `1` exactly when the congruence + holds for every `1 ≤ a ≤ A`, within `time r n A` instructions. + +## References + +* [M. Agrawal, N. Kayal, N. Saxena, *PRIMES is in P*][AKS2004] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.AKS.CongruenceLoop + +open Cslib.WordRAM Cslib.AKS + +/-- Base address of the result block. -/ +abbrev R : Fin 29 := 0 +/-- Base address of the base block. -/ +abbrev B : Fin 29 := 1 +/-- Base address of the scratch block. -/ +abbrev T : Fin 29 := 2 +/-- The number `r` of coefficients. -/ +abbrev rr : Fin 29 := 3 +/-- The modulus `n`. -/ +abbrev nn : Fin 29 := 4 +/-- The tested value `a`, which is also the loop index. -/ +abbrev av : Fin 29 := 21 +/-- The answer of one test. -/ +abbrev flag : Fin 29 := 24 +/-- The last tested value `A`, preserved. -/ +abbrev A : Fin 29 := 25 +/-- The loop end `A + 1`. -/ +abbrev lim : Fin 29 := 26 +/-- `1` while every test so far has passed. -/ +abbrev ok : Fin 29 := 27 +/-- Constant one. -/ +abbrev c1 : Fin 29 := 28 + +/-- The test for one `a` runs on the first `25` registers. -/ +def congMap : Fin 25 → Fin 29 := Fin.castLE (by norm_num) + +theorem congMap_injective : congMap.Injective := Fin.castLE_injective _ + +/-- The test for one `a`. -/ +def cong : Code .standard 29 := Congruence.program.rename congMap + +/-- Test `a`, record a failure, and advance `a`. -/ +def body : Code .standard 29 := + .seq cong (ram[ if flag = c1 then [ ok ← ok ] else [ ok ← #0 ]; av ← av + c1 ]) + +/-- Test `a = 1, …, A`. -/ +def program : Code .standard 29 := + .seq (ram[ c1 ← #1; lim ← A + c1; av ← c1; ok ← #1 ]) (.loop (.lt av lim) body) + +/-- The instruction bound of `program`. -/ +def time (r n A : ℕ) : ℕ := A * (Congruence.time r n + 6) + 6 + +/-- The registers that the test for one `a` preserves. -/ +def Keep (j : Fin 29) : Prop := + j = R ∨ j = B ∨ j = T ∨ j = rr ∨ j = nn ∨ j = av ∨ j = A ∨ j = lim ∨ j = ok ∨ j = c1 + +instance : DecidablePred Keep := fun _ => by unfold Keep; infer_instance + +private theorem cong_preserves : ∀ j, Keep j → cong.Preserves j := by decide + +private theorem cong_exec {w : ℕ} (x : Data 29) {a b t r n v : ℕ} (ha : x.regs R = a) + (hb : x.regs B = b) (ht : x.regs T = t) (hr : x.regs rr = r) (hn : x.regs nn = n) + (hv : x.regs av = v) (hr0 : 0 < r) (hn0 : 0 < n) (hab : a + r ≤ b ∨ b + r ≤ a) + (hta : t + r ≤ a ∨ a + r ≤ t) (htb : t + r ≤ b ∨ b + r ≤ t) (hfa : Fits w (a + r)) + (hfb : Fits w (b + r)) (hft : Fits w (t + r)) (hf2 : Fits w (2 * r)) + (hfn : Fits w (n * n + n)) (hfv : Fits w v) : + ∃ y, ExecWithin w cong x y (Congruence.time r n) ∧ y.regs flag ≤ 1 ∧ + (y.regs flag = 1 ↔ Congruence n r v) ∧ ∀ j, Keep j → y.regs j = x.regs j := by + obtain ⟨y, hy, hle, hiff, -⟩ := Congruence.exec (x.comap congMap) ha hb ht hr hn hv hr0 hn0 hab + hta htb hfa hfb hft hf2 hfn hfv + have hy' := hy.rename congMap_injective + have hfl : (y.extend congMap x).regs flag = y.regs Congruence.flag := + congMap_injective.extend_apply _ _ Congruence.flag + exact ⟨_, hy', by rw [hfl]; exact hle, by rw [hfl]; exact hiff, + fun j hj => hy'.regs_eq (cong_preserves j hj)⟩ + +/-- **Step 5.** For `r ≥ 1` and `n ≥ 1`, with the three blocks pairwise disjoint and the word size +holding the block ends, `2 r`, `n² + n` and `A + 1`, the program leaves `1` in `ok` exactly when +`(X + a)^n ≡ X^n + a (mod X^r - 1, n)` for every `1 ≤ a ≤ A`, within `time r n A` instructions. -/ +theorem exec {w : ℕ} (x : Data 29) {a b t r n A₀ : ℕ} (ha : x.regs R = a) (hb : x.regs B = b) + (ht : x.regs T = t) (hr : x.regs rr = r) (hn : x.regs nn = n) (hA : x.regs A = A₀) + (hr0 : 0 < r) (hn0 : 0 < n) (hab : a + r ≤ b ∨ b + r ≤ a) (hta : t + r ≤ a ∨ a + r ≤ t) + (htb : t + r ≤ b ∨ b + r ≤ t) (hfa : Fits w (a + r)) (hfb : Fits w (b + r)) + (hft : Fits w (t + r)) (hf2 : Fits w (2 * r)) (hfn : Fits w (n * n + n)) + (hfA : Fits w (A₀ + 1)) : + ∃ y, ExecWithin w program x y (time r n A₀) ∧ y.regs ok ≤ 1 ∧ + (y.regs ok = 1 ↔ ∀ v, 1 ≤ v → v ≤ A₀ → Congruence n r v) ∧ + y.regs R = a ∧ y.regs B = b ∧ y.regs T = t ∧ y.regs rr = r ∧ y.regs nn = n := by + classical + have h1 : Fits w 1 := hfA.of_le (by omega) + let I (y : Data 29) : Prop := y.regs R = a ∧ y.regs B = b ∧ y.regs T = t ∧ y.regs rr = r ∧ + y.regs nn = n ∧ y.regs lim = A₀ + 1 ∧ y.regs c1 = 1 ∧ 1 ≤ y.regs av ∧ + y.regs av ≤ A₀ + 1 ∧ y.regs ok ≤ 1 ∧ + (y.regs ok = 1 ↔ ∀ v, 1 ≤ v → v < y.regs av → Congruence n r v) + have hstep (y : Data 29) (hy : I y) (hc : (Cond.lt av lim).Holds y) : + ∃ z, ExecWithin w body y z (Congruence.time r n + 4) ∧ I z ∧ + A₀ + 1 - z.regs av < A₀ + 1 - y.regs av := by + obtain ⟨hyR, hyB, hyT, hyr, hyn, hyl, hy1, hv1, hvA, hyok, hyiff⟩ := hy + obtain ⟨V, hV⟩ : ∃ V, y.regs av = V := ⟨_, rfl⟩ + have hVA : V ≤ A₀ := by simp only [Cond.holds_lt, hV, hyl] at hc; omega + rw [hV] at hv1 hyiff + obtain ⟨y₁, hrun₁, hle, hiff, hkeep⟩ := cong_exec (w := w) y hyR hyB hyT hyr hyn hV hr0 hn0 + hab hta htb hfa hfb hft hf2 hfn (hfA.of_le (by omega)) + have hk (j : Fin 29) (hj : Keep j) : y₁.regs j = y.regs j := hkeep j hj + let o := if Congruence n r V then y.regs ok else 0 + let z := (y₁.set ok o).set av (V + 1) + have hy₁1 : y₁.regs c1 = 1 := by rw [hk c1 (by simp [Keep]), hy1] + have hy₁ok : y₁.regs ok = y.regs ok := hk ok (by simp [Keep]) + have hy₁v : y₁.regs av = V := by rw [hk av (by simp [Keep]), hV] + have hfV : Fits w (V + 1) := hfA.of_le (by omega) + have hrun₂ : ExecWithin (B := .standard) w + (ram[ if flag = c1 then [ ok ← ok ] else [ ok ← #0 ]; av ← av + c1 ]) y₁ z 4 := by + by_cases hC : Congruence n r V + · have hf1 : y₁.regs flag = 1 := hiff.mpr hC + refine ExecWithin.mono (T := 3) ?_ (by omega) + ram_eval [z, o, hC, hf1, hy₁1, hy₁ok, hy₁v, hfV] + · have hf1 : ¬ y₁.regs flag = 1 := fun h => hC (hiff.mp h) + ram_eval [z, o, hC, hf1, hy₁1, hy₁v, hfV] + have hz (j : Fin 29) (hj : Keep j) (hjo : j ≠ ok) (hjv : j ≠ av) : z.regs j = y.regs j := by + simp only [z, Data.set_regs, Function.update_of_ne hjv, Function.update_of_ne hjo] + exact hk j hj + have hzv : z.regs av = V + 1 := by simp [z] + have hzo : z.regs ok = o := by simp +decide [z] + refine ⟨z, hrun₁.seq hrun₂, ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩, ?_⟩ + · rw [hz R (by simp [Keep]) (by decide) (by decide), hyR] + · rw [hz B (by simp [Keep]) (by decide) (by decide), hyB] + · rw [hz T (by simp [Keep]) (by decide) (by decide), hyT] + · rw [hz rr (by simp [Keep]) (by decide) (by decide), hyr] + · rw [hz nn (by simp [Keep]) (by decide) (by decide), hyn] + · rw [hz lim (by simp [Keep]) (by decide) (by decide), hyl] + · rw [hz c1 (by simp [Keep]) (by decide) (by decide), hy1] + · rw [hzv]; omega + · rw [hzv]; omega + · rw [hzo]; simp only [o]; split_ifs <;> omega + · rw [hzo, hzv] + simp only [o] + constructor + · intro h v hv1' hv + split_ifs at h with hC + rcases Nat.lt_succ_iff_lt_or_eq.mp hv with hv | rfl + · exact hyiff.mp h v hv1' hv + · exact hC + · intro h + have hC : Congruence n r V := h V hv1 (by omega) + simp only [hC, ↓reduceIte] + exact hyiff.mpr fun v hv1' hv => h v hv1' (by omega) + · rw [hzv, hV]; omega + let x₀ := (((x.set c1 1).set lim (A₀ + 1)).set av 1).set ok 1 + have hrun₀ : ExecWithin (B := .standard) w (ram[ c1 ← #1; lim ← A + c1; av ← c1; ok ← #1 ]) + x x₀ 4 := by + ram_eval [x₀, hA, h1, hfA] + have hx₀ : I x₀ := by + refine ⟨by simp +decide [x₀, ha], by simp +decide [x₀, hb], by simp +decide [x₀, ht], + by simp +decide [x₀, hr], by simp +decide [x₀, hn], by simp +decide [x₀], + by simp +decide [x₀], by simp +decide [x₀], by simp +decide [x₀], by simp +decide [x₀], ?_⟩ + have hv : x₀.regs av = 1 := by simp +decide [x₀] + have ho : x₀.regs ok = 1 := by simp +decide [x₀] + rw [hv, ho] + exact ⟨fun _ v hv1 hv => by omega, fun _ => rfl⟩ + obtain ⟨y, hrun, hy, hc⟩ := ExecWithin.loop_invariant I (fun y => A₀ + 1 - y.regs av) hstep hx₀ + obtain ⟨hyR, hyB, hyT, hyr, hyn, hyl, -, -, hvA, hyok, hyiff⟩ := hy + have hend : y.regs av = A₀ + 1 := by simp only [Cond.holds_lt, hyl] at hc; omega + refine ⟨y, (hrun₀.seq hrun).mono ?_, hyok, ?_, hyR, hyB, hyT, hyr, hyn⟩ + · have hv : x₀.regs av = 1 := by simp +decide [x₀] + rw [hv, show A₀ + 1 - 1 = A₀ by omega, show Congruence.time r n + 4 + 2 = + Congruence.time r n + 6 by ring] + unfold time + omega + · rw [hyiff, hend] + exact ⟨fun h v hv1 hv => h v hv1 (by omega), fun h v hv1 hv => h v hv1 (by omega)⟩ + +end Cslib.Algorithms.WordRAM.AKS.CongruenceLoop diff --git a/Complexity/Algorithms/WordRAM/AKS/GcdTotient.lean b/Complexity/Algorithms/WordRAM/AKS/GcdTotient.lean new file mode 100644 index 000000000..055f52590 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/AKS/GcdTotient.lean @@ -0,0 +1,256 @@ +/- +Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Samuel Schlesinger +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.NatGcd +public import Complexity.Computability.WordRAM.RegisterFrame +public import Mathlib.Data.Nat.Totient +import Mathlib.Tactic.FinCases + +/-! +# Step 3 of the AKS algorithm and Euler's totient on the word RAM + +Step 3 of the AKS algorithm [Agrawal, Kayal, Saxena, Section 4][AKS2004] outputs COMPOSITE if +`1 < gcd(a, n) < n` for some `a ≤ r`, and step 5 needs `φ(r)`. One loop over `0 ≤ a ≤ r` does +both, with two runs of Euclid's algorithm (`Cslib.Algorithms.WordRAM.NatGcd.program`) per `a`: it +clears a flag when `gcd(n, a)` is neither `1` nor `n`, and counts the `a` with `gcd(r, a) = 1`. +The index `a = 0` is harmless (`gcd(n, 0) = n`), and for `r ≥ 2` neither `0` nor `r` is coprime +to `r`, so the count is `φ(r)`. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.AKS.GcdTotient.program`: the loop. +* `Cslib.Algorithms.WordRAM.AKS.GcdTotient.time`: its instruction bound. + +## Main results + +* `Cslib.Algorithms.WordRAM.AKS.GcdTotient.exec`: the flag records step 3 and the count is `φ(r)`, + within `time r` instructions. + +## References + +* [M. Agrawal, N. Kayal, N. Saxena, *PRIMES is in P*][AKS2004] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.AKS.GcdTotient + +open Cslib.WordRAM + +/-- The number `n`, preserved. -/ +abbrev n : Fin 11 := 0 +/-- The modulus `r`, preserved. -/ +abbrev r : Fin 11 := 1 +/-- The loop index `a`. -/ +abbrev a : Fin 11 := 2 +/-- The loop end `r + 1`. -/ +abbrev lim : Fin 11 := 3 +/-- `1` while every `gcd(n, a)` seen is `1` or `n`. -/ +abbrev ok : Fin 11 := 4 +/-- The number of `a` seen with `gcd(r, a) = 1`. -/ +abbrev phi : Fin 11 := 5 +/-- Constant one. -/ +abbrev one : Fin 11 := 6 +/-- First argument of the gcd, and its result. -/ +abbrev ga : Fin 11 := 7 +/-- Second argument of the gcd. -/ +abbrev gb : Fin 11 := 8 + +/-- Euclid's algorithm runs on `ga`, `gb` and two scratch registers. -/ +def gcdMap : Fin 4 → Fin 11 := ![ga, gb, 9, 10] + +theorem gcdMap_injective : gcdMap.Injective := by decide + +/-- The gcd on `ga` and `gb`. -/ +def gcd : Code .standard 11 := NatGcd.program.rename gcdMap + +/-- Test `gcd(n, a)`, count `gcd(r, a) = 1`, and advance `a`. -/ +def body : Code .standard 11 := + .seq (ram[ ga ← n; gb ← a ]) (.seq gcd + (.seq (ram[ if ga = one then [ ok ← ok ] else [ if ga = n then [ ok ← ok ] else [ ok ← #0 ] ]; + ga ← r; gb ← a ]) (.seq gcd + (ram[ if ga = one then [ phi ← phi + one ] else [ phi ← phi ]; a ← a + one ])))) + +/-- Run the body for `a = 0, 1, …, r`. -/ +def program : Code .standard 11 := + .seq (ram[ one ← #1; lim ← r + one; a ← #0; ok ← #1; phi ← #0 ]) (.loop (.lt a lim) body) + +/-- The instruction bound of `program`. -/ +def time (r : ℕ) : ℕ := (r + 1) * (10 * r + 21) + 7 + +/-- The gcd changes only `ga`, `gb` and its scratch registers. -/ +private theorem gcd_regs {w T : ℕ} {x y : Data 11} (h : ExecWithin w gcd x y T) (j : Fin 11) + (hj : j.val < 7) : y.regs j = x.regs j := by + refine h.regs_eq (Code.preserves_rename_of_forall_ne _ _ _ ?_) + intro b hb + subst hb + fin_cases b <;> simp [gcdMap] at hj + +private theorem gcd_exec {w : ℕ} (x : Data 11) {u v : ℕ} (hu : x.regs ga = u) (hv : x.regs gb = v) + (hfu : Fits w u) (hfv : Fits w v) : + ∃ y, ExecWithin w gcd x y (5 * v + 3) ∧ y.regs ga = Nat.gcd u v ∧ + ∀ j : Fin 11, j.val < 7 → y.regs j = x.regs j := by + obtain ⟨y, hy, hg, -⟩ := NatGcd.exec_rename gcdMap gcdMap_injective x + (by rw [show gcdMap NatGcd.a = ga from rfl, hu]; exact hfu) + (by rw [show gcdMap NatGcd.b = gb from rfl, hv]; exact hfv) + refine ⟨y, ?_, ?_, fun j hj => gcd_regs hy j hj⟩ + · rw [show gcdMap NatGcd.b = gb from rfl, hv] at hy; exact hy + · rw [show ga = gcdMap NatGcd.a from rfl, hg, show gcdMap NatGcd.a = ga from rfl, + show gcdMap NatGcd.b = gb from rfl, hu, hv] + +/-- Counting the residues coprime to `r ≥ 2` up to `r` inclusive gives `φ(r)`. -/ +private theorem card_filter_coprime_range_succ {R : ℕ} (hR : 2 ≤ R) : + ((Finset.range (R + 1)).filter R.Coprime).card = Nat.totient R := by + rw [Finset.range_add_one, Finset.filter_insert] + simp [Nat.coprime_self, show R ≠ 1 by omega, Nat.totient_eq_card_coprime] + +/-- The loop invariant. -/ +private def Inv (N R : ℕ) (y : Data 11) : Prop := + y.regs n = N ∧ y.regs r = R ∧ y.regs one = 1 ∧ y.regs lim = R + 1 ∧ y.regs a ≤ R + 1 ∧ + y.regs ok ≤ 1 ∧ (y.regs ok = 1 ↔ ∀ b < y.regs a, Nat.gcd N b = 1 ∨ Nat.gcd N b = N) ∧ + y.regs phi = ((Finset.range (y.regs a)).filter R.Coprime).card + +private theorem body_step {w N R : ℕ} (hfN : Fits w N) (hfR : Fits w (R + 1)) (y : Data 11) + (hy : Inv N R y) (hc : (Cond.lt a lim).Holds y) : + ∃ z, ExecWithin w body y z (10 * R + 19) ∧ Inv N R z ∧ + R + 1 - z.regs a < R + 1 - y.regs a := by + classical + obtain ⟨hyn, hyr, hy1, hyl, hya, hyok, hyiff, hyphi⟩ := hy + obtain ⟨A, hA⟩ : ∃ A, y.regs a = A := ⟨_, rfl⟩ + have hAR : A ≤ R := by simp only [Cond.holds_lt, hA, hyl] at hc; omega + rw [hA] at hyiff hyphi + have hfA : Fits w A := hfR.of_le (by omega) + have h1f : Fits w 1 := hfR.of_le (by omega) + -- `gcd(n, a)`. + let y₁ := (y.set ga N).set gb A + have hrun₁ : ExecWithin (B := .standard) w (ram[ ga ← n; gb ← a ]) y y₁ 2 := by + ram_eval [y₁, hyn, hA] + have hy₁ (j : Fin 11) (hj : j.val < 7) : y₁.regs j = y.regs j := by + have h9 : j ≠ ga ∧ j ≠ gb := by constructor <;> (intro h; subst h; simp at hj) + simp [y₁, h9.1, h9.2] + obtain ⟨y₂, hrun₂, hg₂, hy₂⟩ := gcd_exec (w := w) y₁ (u := N) (v := A) (by simp +decide [y₁]) + (by simp +decide [y₁]) hfN hfA + have hy₂' (j : Fin 11) (hj : j.val < 7) : y₂.regs j = y.regs j := by rw [hy₂ j hj, hy₁ j hj] + -- Record step 3. + let o := if Nat.gcd N A = 1 ∨ Nat.gcd N A = N then y.regs ok else 0 + let y₃ := ((y₂.set ok o).set ga R).set gb A + have hy₂ok : y₂.regs ok = y.regs ok := hy₂' ok (by simp) + have hy₂n : y₂.regs n = N := by rw [hy₂' n (by simp), hyn] + have hy₂1 : y₂.regs one = 1 := by rw [hy₂' one (by simp), hy1] + have hy₂r : y₂.regs r = R := by rw [hy₂' r (by simp), hyr] + have hy₂a : y₂.regs a = A := by rw [hy₂' a (by simp), hA] + have hrun₃ : ExecWithin (B := .standard) w + (ram[ if ga = one then [ ok ← ok ] else [ if ga = n then [ ok ← ok ] else [ ok ← #0 ] ]; + ga ← r; gb ← a ]) y₂ y₃ 7 := by + by_cases hg1 : Nat.gcd N A = 1 + · refine ExecWithin.mono (T := 4) ?_ (by omega) + ram_eval [y₃, o, hg₂, hg1, hy₂1, hy₂ok, hy₂r, hy₂a] + · by_cases hgN : Nat.gcd N A = N + · have hN1 : ¬ N = 1 := fun h => hg1 (hgN.trans h) + refine ExecWithin.mono (T := 6) ?_ (by omega) + ram_eval [y₃, o, hg₂, hg1, hgN, hN1, hy₂1, hy₂n, hy₂ok, hy₂r, hy₂a] + · ram_eval [y₃, o, hg₂, hg1, hgN, hy₂1, hy₂n, hy₂r, hy₂a] + -- `gcd(r, a)`. + obtain ⟨y₄, hrun₄, hg₄, hy₄⟩ := gcd_exec (w := w) y₃ (u := R) (v := A) (by simp +decide [y₃]) + (by simp +decide [y₃]) (hfR.of_le (by omega)) hfA + have hy₄' (j : Fin 11) (hj : j.val < 7) (hjo : j ≠ ok) : y₄.regs j = y.regs j := by + rw [hy₄ j hj] + have h9 : j ≠ ga ∧ j ≠ gb := by constructor <;> (intro h; subst h; simp at hj) + simp only [y₃, Data.set_regs, Function.update_of_ne h9.1, Function.update_of_ne h9.2, + Function.update_of_ne hjo] + exact hy₂' j hj + have hy₄ok : y₄.regs ok = o := by rw [hy₄ ok (by simp)]; simp +decide [y₃] + have hy₄1 : y₄.regs one = 1 := by rw [hy₄' one (by simp) (by decide), hy1] + have hy₄phi : y₄.regs phi = y.regs phi := hy₄' phi (by simp) (by decide) + have hy₄a : y₄.regs a = A := by rw [hy₄' a (by simp) (by decide), hA] + -- Count and advance. + have hphiA : ((Finset.range A).filter R.Coprime).card ≤ A := + (Finset.card_filter_le _ _).trans (Finset.card_range A).le + let ph := if Nat.gcd R A = 1 then y.regs phi + 1 else y.regs phi + let z := (y₄.set phi ph).set a (A + 1) + have hfph : Fits w (y.regs phi + 1) := hfR.of_le (by omega) + have hfA1 : Fits w (A + 1) := hfR.of_le (by omega) + have hrun₅ : ExecWithin (B := .standard) w + (ram[ if ga = one then [ phi ← phi + one ] else [ phi ← phi ]; a ← a + one ]) y₄ z 4 := by + by_cases hg1 : Nat.gcd R A = 1 + · refine ExecWithin.mono (T := 3) ?_ (by omega) + ram_eval [z, ph, hg₄, hg1, hy₄1, hy₄phi, hy₄a, hfph, hfA1] + · ram_eval [z, ph, hg₄, hg1, hy₄1, hy₄phi, hy₄a, hfA1] + have hza : z.regs a = A + 1 := by simp [z] + refine ⟨z, (hrun₁.seq (hrun₂.seq (hrun₃.seq (hrun₄.seq hrun₅)))).mono (by omega), ?_, ?_⟩ + · refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩ + · simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₄' n (by simp) (by decide), hyn] + · simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₄' r (by simp) (by decide), hyr] + · simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte]; exact hy₄1 + · simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₄' lim (by simp) (by decide), hyl] + · rw [hza]; omega + · simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₄ok]; simp only [o]; split_ifs <;> omega + · have hzo : z.regs ok = o := by simp +decide [z, hy₄ok] + rw [hzo, hza] + simp only [o] + constructor + · intro h b hb + split_ifs at h with hc + · rcases Nat.lt_succ_iff_lt_or_eq.mp hb with hb | rfl + · exact hyiff.mp h b hb + · exact hc + · intro h + have hc : Nat.gcd N A = 1 ∨ Nat.gcd N A = N := h A (by omega) + simp only [hc, ↓reduceIte] + exact hyiff.mpr fun b hb => h b (by omega) + · have hzp : z.regs phi = ph := by simp +decide [z] + rw [hzp, hza, Finset.range_add_one, Finset.filter_insert] + simp only [ph, hyphi, Nat.Coprime] + split_ifs with hg + · rw [Finset.card_insert_of_notMem (by simp)] + · rfl + · rw [hza, hA]; omega + +/-- **Step 3 and the totient.** For `r ≥ 2`, with the word size holding `n` and `r + 1`, within +`time r` instructions the flag `ok` is `1` exactly when every `1 ≤ a ≤ r` has `gcd(a, n) = 1` or +`gcd(a, n) = n`, and `phi` is `φ(r)`. For `n ≥ 1` the flag condition says that no `1 ≤ a ≤ r` has +`1 < gcd(a, n) < n`; for `n = 0` it is stronger, since `gcd(a, 0) = a`. -/ +theorem exec {w : ℕ} (x : Data 11) {N R : ℕ} (hn : x.regs n = N) (hr : x.regs r = R) + (hR : 2 ≤ R) (hfN : Fits w N) (hfR : Fits w (R + 1)) : + ∃ y, ExecWithin w program x y (time R) ∧ y.regs ok ≤ 1 ∧ + (y.regs ok = 1 ↔ ∀ b, 1 ≤ b → b ≤ R → b.gcd N = 1 ∨ b.gcd N = N) ∧ + y.regs phi = Nat.totient R ∧ y.regs n = N ∧ y.regs r = R := by + have h1f : Fits w 1 := hfR.of_le (by omega) + let x₀ := ((((x.set one 1).set lim (R + 1)).set a 0).set ok 1).set phi 0 + have hrun₀ : ExecWithin (B := .standard) w + (ram[ one ← #1; lim ← r + one; a ← #0; ok ← #1; phi ← #0 ]) x x₀ 5 := by + ram_eval [x₀, hr, h1f, hfR] + have hx₀ : Inv N R x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩ <;> simp +decide [x₀, hn, hr] + obtain ⟨y, hrun, hy, hc⟩ := ExecWithin.loop_invariant (Inv N R) (fun y => R + 1 - y.regs a) + (body_step hfN hfR) hx₀ + obtain ⟨hyn, hyr, -, hyl, hya, hyok, hyiff, hyphi⟩ := hy + have hend : y.regs a = R + 1 := by simp only [Cond.holds_lt, hyl] at hc; omega + rw [hend] at hyiff hyphi + refine ⟨y, (hrun₀.seq hrun).mono ?_, hyok, ?_, ?_, hyn, hyr⟩ + · have : x₀.regs a = 0 := by simp +decide [x₀] + rw [this, Nat.sub_zero, show 10 * R + 19 + 2 = 10 * R + 21 by ring] + unfold time + omega + · rw [hyiff] + constructor + · intro h b hb1 hbR + rw [Nat.gcd_comm] + exact h b (by omega) + · intro h b hb + rcases Nat.eq_zero_or_pos b with rfl | hb0 + · right; exact Nat.gcd_zero_right N + · rw [Nat.gcd_comm]; exact h b hb0 (by omega) + · rw [hyphi, card_filter_coprime_range_succ hR] + +end Cslib.Algorithms.WordRAM.AKS.GcdTotient diff --git a/Complexity/Algorithms/WordRAM/AKS/OrderSearch.lean b/Complexity/Algorithms/WordRAM/AKS/OrderSearch.lean new file mode 100644 index 000000000..1e6c7f8aa --- /dev/null +++ b/Complexity/Algorithms/WordRAM/AKS/OrderSearch.lean @@ -0,0 +1,293 @@ +/- +Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Samuel Schlesinger +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.NatGcd +public import Complexity.Computability.WordRAM.RegisterFrame +import Mathlib.Tactic.Linarith +import Mathlib.Tactic.FinCases + +/-! +# The search for `r` in the AKS algorithm on the word RAM + +Step 2 of the AKS algorithm [Agrawal, Kayal, Saxena, Section 4][AKS2004] finds the least `r` such +that the multiplicative order of `n` modulo `r` exceeds a bound `L`. As in the proof of their +Theorem 5.1, the program tries `r = 2, 3, …` in turn. It accepts `r` when `gcd(n, r) = 1` and +`n^k mod r ≠ 1` for every `1 ≤ k ≤ L`; for `r ≥ 2` this is the order condition +(`Cslib.AKS.lt_orderOf_iff`). The powers are computed by one multiplication modulo `r` each, and +the gcd by Euclid's algorithm (`Cslib.Algorithms.WordRAM.NatGcd.program`). Given any `r₀ ≥ 2` +satisfying the condition, the search stops at some `2 ≤ r ≤ r₀` satisfying it, within a number of +instructions polynomial in `r₀` and `L`. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.AKS.OrderSearch.Good`: the condition tested for one `r`. +* `Cslib.Algorithms.WordRAM.AKS.OrderSearch.program`: the search. +* `Cslib.Algorithms.WordRAM.AKS.OrderSearch.time`: its instruction bound. + +## Main results + +* `Cslib.Algorithms.WordRAM.AKS.OrderSearch.exec`: the search stops at an `r ≤ r₀` satisfying the + condition, within `time L r₀` instructions. + +## References + +* [M. Agrawal, N. Kayal, N. Saxena, *PRIMES is in P*][AKS2004] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.AKS.OrderSearch + +open Cslib.WordRAM + +/-- The condition tested for `r`: `n` is coprime to `r` and no `n^k` with `1 ≤ k ≤ L` is `1` +modulo `r`. -/ +def Good (n L r : ℕ) : Prop := n.Coprime r ∧ ∀ k, 1 ≤ k → k ≤ L → n ^ k % r ≠ 1 + +/-- The number `n`, preserved. -/ +abbrev n : Fin 13 := 0 +/-- The bound `L` on the order, preserved. -/ +abbrev L : Fin 13 := 1 +/-- The candidate `r`, and the output. -/ +abbrev r : Fin 13 := 2 +/-- `1` once a candidate is accepted. -/ +abbrev found : Fin 13 := 3 +/-- Constant one. -/ +abbrev one : Fin 13 := 4 +/-- `n mod r`. -/ +abbrev nr : Fin 13 := 5 +/-- The power `n^k mod r`. -/ +abbrev p : Fin 13 := 6 +/-- The exponent `k`. -/ +abbrev k : Fin 13 := 7 +/-- `1` while no power `n^j` with `1 ≤ j ≤ k` is `1` modulo `r`. -/ +abbrev ok : Fin 13 := 8 +/-- First argument of the gcd, and its result. -/ +abbrev ga : Fin 13 := 9 +/-- Second argument of the gcd. -/ +abbrev gb : Fin 13 := 10 + +/-- Euclid's algorithm runs on `ga`, `gb` and two scratch registers. -/ +def gcdMap : Fin 4 → Fin 13 := ![ga, gb, 11, 12] + +theorem gcdMap_injective : gcdMap.Injective := by decide + +/-- One more power: `p := p · (n mod r) mod r`, and record whether it is `1`. -/ +def powBody : Code .standard 13 := ram[ + p ← p * nr; p ← p % r; k ← k + one; if p = one then [ ok ← #0 ] else [ ok ← ok ] ] + +/-- Test the powers `n^k mod r` for `1 ≤ k ≤ L`. -/ +def powers : Code .standard 13 := + .seq (ram[ nr ← n % r; p ← one; k ← #0; ok ← #1 ]) (.loop (.lt k L) powBody) + +/-- Accept `r` if `gcd(n, r) = 1` and no small power is `1`; otherwise try `r + 1`. -/ +def accept : Code .standard 13 := ram[ + if ga = one then [ if ok = one then [ found ← #1 ] else [ r ← r + one ] ] else [ r ← r + one ] ] + +/-- Test one candidate. -/ +def body : Code .standard 13 := + .seq powers (.seq (ram[ ga ← n; gb ← r ]) (.seq (NatGcd.program.rename gcdMap) accept)) + +/-- Try `r = 2, 3, …` until one is accepted. -/ +def program : Code .standard 13 := + .seq (ram[ one ← #1; r ← #2; found ← #0 ]) (.loop (.lt found one) body) + +/-- The instruction bound of `program`. -/ +def time (L r₀ : ℕ) : ℕ := (2 * r₀ + 1) * (8 * L + 5 * r₀ + 20) + 5 + +/-- The registers that the power test preserves. -/ +private theorem powers_preserves : + ∀ j ∈ [n, L, r, found, one], powers.Preserves j := by decide + +private theorem powers_exec {w : ℕ} (x : Data 13) {N L₀ R : ℕ} (hn : x.regs n = N) + (hL : x.regs L = L₀) (hr : x.regs r = R) (h1 : x.regs one = 1) (hR : 2 ≤ R) + (hfN : Fits w N) (hfR : Fits w (R * R)) (hfL : Fits w (L₀ + 1)) : + ∃ y, ExecWithin w powers x y (8 * L₀ + 6) ∧ y.regs ok ≤ 1 ∧ + (y.regs ok = 1 ↔ ∀ j, 1 ≤ j → j ≤ L₀ → N ^ j % R ≠ 1) ∧ + ∀ j ∈ [n, L, r, found, one], y.regs j = x.regs j := by + have hR0 : 0 < R := by omega + have h1f : Fits w 1 := hfL.of_le (by omega) + let J (y : Data 13) : Prop := y.regs L = L₀ ∧ y.regs r = R ∧ y.regs one = 1 ∧ + y.regs nr = N % R ∧ y.regs k ≤ L₀ ∧ y.regs p = N ^ y.regs k % R ∧ y.regs ok ≤ 1 ∧ + (y.regs ok = 1 ↔ ∀ j, 1 ≤ j → j ≤ y.regs k → N ^ j % R ≠ 1) + have hstep (y : Data 13) (hy : J y) (hc : (Cond.lt k L).Holds y) : + ∃ z, ExecWithin w powBody y z 6 ∧ J z ∧ L₀ - z.regs k < L₀ - y.regs k := by + obtain ⟨hyL, hyr, hy1, hynr, hyk, hyp, hyok, hyiff⟩ := hy + obtain ⟨K, hK⟩ : ∃ K, y.regs k = K := ⟨_, rfl⟩ + have hKL : K < L₀ := by simpa only [Cond.holds_lt, hK, hyL] using hc + rw [hK] at hyp hyiff + have hPR : N ^ K % R < R := Nat.mod_lt _ hR0 + have hNR : N % R < R := Nat.mod_lt _ hR0 + have hmul : N ^ K % R * (N % R) < R * R := Nat.mul_lt_mul'' hPR hNR + have hmod : N ^ K % R * (N % R) % R = N ^ (K + 1) % R := by + rw [← Nat.mul_mod, pow_succ] + have hf₁ : Fits w (N ^ K % R * (N % R)) := hfR.of_le hmul.le + have hf₂ : Fits w (N ^ (K + 1) % R) := hfR.of_le + ((Nat.mod_lt _ hR0).le.trans (Nat.le_mul_self R)) + have hf₃ : Fits w (K + 1) := hfL.of_le (by omega) + have hiff (o : ℕ) : ((∀ j, 1 ≤ j → j ≤ K + 1 → N ^ j % R ≠ 1) ↔ + (∀ j, 1 ≤ j → j ≤ K → N ^ j % R ≠ 1) ∧ N ^ (K + 1) % R ≠ 1) := by + constructor + · intro h; exact ⟨fun j hj hjK => h j hj (by omega), h (K + 1) (by omega) le_rfl⟩ + · rintro ⟨h, hK1⟩ j hj hjK + rcases Nat.lt_or_ge j (K + 1) with hlt | hge + · exact h j hj (by omega) + · rw [show j = K + 1 by omega]; exact hK1 + let o := if N ^ (K + 1) % R = 1 then 0 else y.regs ok + let z := (((y.set p (N ^ K % R * (N % R))).set p (N ^ (K + 1) % R)).set k (K + 1)).set ok o + have hrun : ExecWithin (B := .standard) w powBody y z 6 := by + by_cases hP : N ^ (K + 1) % R = 1 + · refine ExecWithin.mono (T := 5) ?_ (by omega) + ram_eval [powBody, z, o, hyp, hynr, hyr, hK, hy1, hmod, hP, hf₁, hf₂, hf₃, h1f] + · ram_eval [powBody, z, o, hyp, hynr, hyr, hK, hy1, hmod, hP, hf₁, hf₂, hf₃, h1f] + have hzk : z.regs k = K + 1 := by simp +decide [z] + refine ⟨z, hrun, ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩, ?_⟩ + · simp +decide [z, hyL] + · simp +decide [z, hyr] + · simp +decide [z, hy1] + · simp +decide [z, hynr] + · rw [hzk]; omega + · rw [hzk]; simp +decide [z] + · simp only [z, o, Data.set_regs, Function.update_self] + split_ifs <;> omega + · rw [hzk, hiff 0] + simp only [z, o, Data.set_regs, Function.update_self] + split_ifs with hP + · simp [hP] + · rw [hyiff]; simp [hP] + · rw [hzk, hK]; omega + let x₀ := (((x.set nr (N % R)).set p 1).set k 0).set ok 1 + have hrun₀ : ExecWithin (B := .standard) w (ram[ nr ← n % r; p ← one; k ← #0; ok ← #1 ]) + x x₀ 4 := by + ram_eval [x₀, hn, hr, h1, h1f, hfN.of_le (Nat.mod_le N R)] + have hx₀ : J x₀ := by + have hk : x₀.regs k = 0 := by simp +decide [x₀] + have hp : x₀.regs p = 1 := by simp +decide [x₀] + have ho : x₀.regs ok = 1 := by simp +decide [x₀] + refine ⟨by simp +decide [x₀, hL], by simp +decide [x₀, hr], by simp +decide [x₀, h1], + by simp +decide [x₀], by rw [hk]; omega, ?_, by rw [ho], ?_⟩ + · rw [hk, hp, pow_zero, Nat.mod_eq_of_lt (by omega : 1 < R)] + · rw [ho, hk] + exact ⟨fun _ j hj hj0 => by omega, fun _ => rfl⟩ + obtain ⟨y, hrun, hy, hc⟩ := ExecWithin.loop_invariant J (fun y => L₀ - y.regs k) hstep hx₀ + obtain ⟨hyL, -, -, -, hyk, -, hyok, hyiff⟩ := hy + have hkL : y.regs k = L₀ := by simp only [Cond.holds_lt, hyL] at hc; omega + have hrun' := hrun₀.seq hrun + refine ⟨y, hrun'.mono ?_, hyok, by rw [hyiff, hkL], + fun j hj => hrun'.regs_eq (powers_preserves j hj)⟩ + have : x₀.regs k = 0 := by simp +decide [x₀] + rw [this] + omega + +/-- **The search for `r`.** If some `r₀ ≥ 2` satisfies `Good n L r₀`, and the word size holds +`n`, `r₀²` and `L + 1`, the search stops within `time L r₀` instructions at some `2 ≤ r ≤ r₀` +satisfying `Good n L r`, preserving `n` and `L`. -/ +theorem exec {w : ℕ} (x : Data 13) {N L₀ r₀ : ℕ} (hn : x.regs n = N) (hL : x.regs L = L₀) + (hr₀ : 2 ≤ r₀) (hgood : Good N L₀ r₀) (hfN : Fits w N) (hfR : Fits w (r₀ * r₀)) + (hfL : Fits w (L₀ + 1)) : + ∃ y, ExecWithin w program x y (time L₀ r₀) ∧ 2 ≤ y.regs r ∧ y.regs r ≤ r₀ ∧ + Good N L₀ (y.regs r) ∧ y.regs n = N ∧ y.regs L = L₀ := by + have h1f : Fits w 1 := hfL.of_le (by omega) + have hfr₁ : Fits w (r₀ + 1) := hfR.of_le (by nlinarith) + let I (y : Data 13) : Prop := y.regs n = N ∧ y.regs L = L₀ ∧ y.regs one = 1 ∧ 2 ≤ y.regs r ∧ + y.regs r ≤ r₀ ∧ y.regs found ≤ 1 ∧ (y.regs found = 1 → Good N L₀ (y.regs r)) + have hstep (y : Data 13) (hy : I y) (hc : (Cond.lt found one).Holds y) : + ∃ z, ExecWithin w body y z (8 * L₀ + 5 * r₀ + 16) ∧ I z ∧ + 2 * (r₀ - z.regs r) + (1 - z.regs found) < 2 * (r₀ - y.regs r) + (1 - y.regs found) := by + obtain ⟨hyn, hyL, hy1, hr2, hrr₀, hyf, -⟩ := hy + obtain ⟨R, hR⟩ : ∃ R, y.regs r = R := ⟨_, rfl⟩ + rw [hR] at hr2 hrr₀ + have hf0 : y.regs found = 0 := by simp only [Cond.holds_lt, hy1] at hc; omega + -- The powers. + obtain ⟨y₁, hrun₁, hok1, hokiff, hy₁⟩ := powers_exec (w := w) y hyn hyL hR hy1 hr2 hfN + (hfR.of_le (Nat.mul_le_mul hrr₀ hrr₀)) hfL + have hy₁n : y₁.regs n = N := by rw [hy₁ n (by simp), hyn] + have hy₁r : y₁.regs r = R := by rw [hy₁ r (by simp), hR] + -- The gcd. + let y₂ := (y₁.set ga N).set gb R + have hrun₂ : ExecWithin (B := .standard) w (ram[ ga ← n; gb ← r ]) y₁ y₂ 2 := by + ram_eval [y₂, hy₁n, hy₁r] + obtain ⟨y₃, hrun₃, hg, -⟩ := NatGcd.exec_rename gcdMap gcdMap_injective y₂ + (by simpa +decide [y₂, gcdMap] using hfN) + (by simpa +decide [y₂, gcdMap] using hfR.of_le (hrr₀.trans (Nat.le_mul_self r₀))) + have hgcd : y₃.regs ga = Nat.gcd N R := by + rw [show ga = gcdMap NatGcd.a from rfl, hg]; simp +decide [y₂, gcdMap] + have hy₃ (j : Fin 13) (hj : j.val < 9) : y₃.regs j = y₁.regs j := by + rw [hrun₃.regs_eq (Code.preserves_rename_of_forall_ne _ _ _ ?_)] + · have h9 : j ≠ ga ∧ j ≠ gb := by + constructor <;> (intro h; subst h; simp at hj) + simp [y₂, h9.1, h9.2] + · intro a h; subst h; fin_cases a <;> simp [gcdMap] at hj + have hy₃1 : y₃.regs one = 1 := by rw [hy₃ one (by simp), hy₁ one (by simp), hy1] + have hy₃r : y₃.regs r = R := by rw [hy₃ r (by simp), hy₁r] + have hy₃ok : y₃.regs ok = y₁.regs ok := hy₃ ok (by simp) + -- Accept or advance. + have hfR1 : Fits w (R + 1) := hfr₁.of_le (by omega) + have hgood' : Good N L₀ R ↔ Nat.gcd N R = 1 ∧ y₁.regs ok = 1 := by + rw [hokiff]; rfl + by_cases hG : Good N L₀ R + · obtain ⟨hg1, ho1⟩ := hgood'.mp hG + let z := y₃.set found 1 + have hrun₄ : ExecWithin (B := .standard) w accept y₃ z 5 := by + refine ExecWithin.mono (T := 3) ?_ (by omega) + ram_eval [accept, z, hgcd, hg1, hy₃1, hy₃ok, ho1, h1f] + refine ⟨z, (hrun₁.seq (hrun₂.seq (hrun₃.seq hrun₄))).mono ?_, ?_, ?_⟩ + · rw [show y₂.regs (gcdMap NatGcd.b) = R by simp +decide [y₂, gcdMap]] + omega + · have hz (j : Fin 13) (hj : j ≠ found) : z.regs j = y₃.regs j := by simp [z, hj] + have hzr : z.regs r = R := by rw [hz r (by decide), hy₃r] + refine ⟨?_, ?_, ?_, by rw [hzr]; exact hr2, by rw [hzr]; exact hrr₀, + by simp +decide [z], fun _ => by rw [hzr]; exact hG⟩ + · rw [hz n (by decide), hy₃ n (by simp), hy₁n] + · rw [hz L (by decide), hy₃ L (by simp), hy₁ L (by simp), hyL] + · rw [hz one (by decide), hy₃1] + · simp +decide [z, hy₃r, hR, hf0] + · have hne : R ≠ r₀ := fun h => hG (h ▸ hgood) + let z := y₃.set r (R + 1) + have hrun₄ : ExecWithin (B := .standard) w accept y₃ z 5 := by + by_cases hg1 : Nat.gcd N R = 1 + · have ho1 : ¬ y₁.regs ok = 1 := fun h => hG (hgood'.mpr ⟨hg1, h⟩) + ram_eval [accept, z, hgcd, hg1, hy₃1, hy₃ok, ho1, hy₃r, hfR1] + · refine ExecWithin.mono (T := 3) ?_ (by omega) + ram_eval [accept, z, hgcd, hg1, hy₃1, hy₃r, hfR1] + refine ⟨z, (hrun₁.seq (hrun₂.seq (hrun₃.seq hrun₄))).mono ?_, ?_, ?_⟩ + · rw [show y₂.regs (gcdMap NatGcd.b) = R by simp +decide [y₂, gcdMap]] + omega + · have hf : z.regs found = 0 := by + simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₃ found (by simp), hy₁ found (by simp), hf0] + have hz (j : Fin 13) (hj : j ≠ r) : z.regs j = y₃.regs j := by simp [z, hj] + have hzr : z.regs r = R + 1 := by simp [z] + refine ⟨?_, ?_, ?_, by rw [hzr]; omega, by rw [hzr]; omega, by omega, by omega⟩ + · rw [hz n (by decide), hy₃ n (by simp), hy₁n] + · rw [hz L (by decide), hy₃ L (by simp), hy₁ L (by simp), hyL] + · rw [hz one (by decide), hy₃1] + · have hf : z.regs found = 0 := by + simp +decide only [z, Data.set_regs, Function.update_apply, ↓reduceIte] + rw [hy₃ found (by simp), hy₁ found (by simp), hf0] + rw [hf, hf0, hR, show z.regs r = R + 1 by simp [z]] + omega + let x₀ := ((x.set one 1).set r 2).set found 0 + have hrun₀ : ExecWithin (B := .standard) w (ram[ one ← #1; r ← #2; found ← #0 ]) x x₀ 3 := by + ram_eval [x₀, h1f, hfr₁.of_le (by omega : 2 ≤ r₀ + 1)] + have hx₀ : I x₀ := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_⟩ <;> simp +decide [x₀, hn, hL, hr₀] + obtain ⟨y, hrun, hy, hc⟩ := ExecWithin.loop_invariant I + (fun y => 2 * (r₀ - y.regs r) + (1 - y.regs found)) hstep hx₀ + obtain ⟨hyn, hyL, hy1, hr2, hrr₀, hyf1, hgoody⟩ := hy + have hf1 : y.regs found = 1 := by simp only [Cond.holds_lt, hy1] at hc; omega + refine ⟨y, (hrun₀.seq hrun).mono ?_, hr2, hrr₀, hgoody hf1, hyn, hyL⟩ + have hx₀r : x₀.regs r = 2 := by simp +decide [x₀] + have hx₀f : x₀.regs found = 0 := by simp +decide [x₀] + have hμ : 2 * (r₀ - x₀.regs r) + (1 - x₀.regs found) ≤ 2 * r₀ + 1 := by omega + have := Nat.mul_le_mul hμ (le_refl (8 * L₀ + 5 * r₀ + 16 + 2)) + unfold time + nlinarith + +end Cslib.Algorithms.WordRAM.AKS.OrderSearch diff --git a/Complexity/Foundations/NumberTheory/AKS/Criterion.lean b/Complexity/Foundations/NumberTheory/AKS/Criterion.lean index 845c948c5..1a2382748 100644 --- a/Complexity/Foundations/NumberTheory/AKS/Criterion.lean +++ b/Complexity/Foundations/NumberTheory/AKS/Criterion.lean @@ -29,8 +29,10 @@ Checking the identity directly takes `n` coefficients. The AKS algorithm checks `X^r - 1`, for a small `r`, and for a few values of `a`. Given `r`, the algorithm outputs PRIME exactly when `n` is not a perfect power, no `a ≤ r` has `1 < gcd(a, n) < n`, and either `n ≤ r` or `(X + a)^n ≡ X^n + a (mod X^r - 1, n)` for every `1 ≤ a ≤ ⌊√φ(r) log₂ n⌋` -(`Cslib.AKS.Passes`). Every prime passes, for every `r` (`Cslib.AKS.passes_of_prime`): this is -the easy direction of [Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004]. The converse, for `r` with +(`Cslib.AKS.Passes`). `Cslib.AKS.PassesUpTo n r A` states the same tests with the congruence +checked for `1 ≤ a ≤ A`, and `Passes n r` is its case `A = ⌊√φ(r) log₂ n⌋`. Every prime passes, +for every `r` and `A` (`Cslib.AKS.passesUpTo_of_prime`, `Cslib.AKS.passes_of_prime`): this is the +easy direction of [Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004]. The converse, for `r` with multiplicative order of `n` modulo `r` above `log₂² n`, is the substance of their theorem. ## Main definitions @@ -38,6 +40,7 @@ multiplicative order of `n` modulo `r` above `log₂² n`, is the substance of t * `Cslib.AKS.IsPerfectPower`: `n = a ^ b` for some `b > 1`. * `Cslib.AKS.numTests`: the number `⌊√φ(r) log₂ n⌋` of tested values of `a`. * `Cslib.AKS.Congruence`: `(X + a)^n ≡ X^n + a (mod X^r - 1, n)`. +* `Cslib.AKS.PassesUpTo`: the AKS tests, with the congruence checked for `1 ≤ a ≤ A`. * `Cslib.AKS.Passes`: the conditions under which the AKS algorithm, given `r`, outputs PRIME. ## Main results @@ -45,7 +48,7 @@ multiplicative order of `n` modulo `r` above `log₂² n`, is the substance of t * `Cslib.AKS.not_dvd_choose_of_prime_dvd`: a composite `n` does not divide `C(n, q)` for a prime factor `q < n`. * `Cslib.AKS.prime_iff_add_pow_eq`: the AKS identity characterizes primes. -* `Cslib.AKS.passes_of_prime`: every prime passes the AKS tests. +* `Cslib.AKS.passesUpTo_of_prime`, `Cslib.AKS.passes_of_prime`: every prime passes the AKS tests. ## References @@ -141,14 +144,19 @@ def Congruence (n r a : ℕ) : Prop := AdjoinRoot.mk (X ^ r - 1 : (ZMod n)[X]) ((X + C (a : ZMod n)) ^ n) = AdjoinRoot.mk (X ^ r - 1 : (ZMod n)[X]) (X ^ n + C (a : ZMod n)) -/-- **The conditions under which the AKS algorithm outputs PRIME**, given `r` -[Agrawal, Kayal, Saxena, Section 4][AKS2004]: `n` is not a perfect power (step 1), no `a ≤ r` has -`1 < gcd(a, n) < n` (step 3), and either `n ≤ r` (step 4) or the congruence holds for every -`1 ≤ a ≤ ⌊√φ(r) log₂ n⌋` (step 5). The algorithm chooses `r` as the least number for which the -multiplicative order of `n` modulo `r` exceeds `log₂² n` (step 2). -/ -def Passes (n r : ℕ) : Prop := +/-- **The AKS tests with the congruence checked up to `A`**: `n` is not a perfect power (step 1), +every `1 ≤ a ≤ r` has `gcd(a, n) = 1` or `gcd(a, n) = n` (step 3; for `n ≥ 1`, no `a ≤ r` has +`1 < gcd(a, n) < n`), and either `n ≤ r` (step 4) or the congruence holds for every `1 ≤ a ≤ A` +(step 5). -/ +def PassesUpTo (n r A : ℕ) : Prop := ¬ IsPerfectPower n ∧ (∀ a, 1 ≤ a → a ≤ r → a.gcd n = 1 ∨ a.gcd n = n) ∧ - (n ≤ r ∨ ∀ a, 1 ≤ a → a ≤ numTests n r → Congruence n r a) + (n ≤ r ∨ ∀ a, 1 ≤ a → a ≤ A → Congruence n r a) + +/-- **The conditions under which the AKS algorithm outputs PRIME**, given `r` +[Agrawal, Kayal, Saxena, Section 4][AKS2004]: the tests of `PassesUpTo`, with the congruence +checked for every `1 ≤ a ≤ ⌊√φ(r) log₂ n⌋` in step 5. The algorithm chooses `r` as the least +number for which the multiplicative order of `n` modulo `r` exceeds `log₂² n` (step 2). -/ +def Passes (n r : ℕ) : Prop := PassesUpTo n r (numTests n r) /-- A prime is not a perfect power. -/ theorem not_isPerfectPower_of_prime {p : ℕ} (hp : p.Prime) : ¬ IsPerfectPower p := by @@ -159,13 +167,18 @@ theorem not_isPerfectPower_of_prime {p : ℕ} (hp : p.Prime) : ¬ IsPerfectPower theorem congruence_of_prime {p : ℕ} (hp : p.Prime) (r a : ℕ) : Congruence p r a := by rw [Congruence, add_pow_eq_of_prime hp] -/-- **Every prime passes the AKS tests**, for every `r`: the easy direction of -[Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004]. -/ -theorem passes_of_prime {p : ℕ} (hp : p.Prime) (r : ℕ) : Passes p r := by +/-- **Every prime passes the AKS tests**, for every `r` and every bound `A` on the tested `a`: +the easy direction of [Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004]. -/ +theorem passesUpTo_of_prime {p : ℕ} (hp : p.Prime) (r A : ℕ) : PassesUpTo p r A := by refine ⟨not_isPerfectPower_of_prime hp, fun a _ _ => ?_, Or.inr fun a _ _ => congruence_of_prime hp r a⟩ rcases (Nat.coprime_or_dvd_of_prime hp a) with h | h · exact Or.inl (Nat.Coprime.symm h) · exact Or.inr (Nat.gcd_eq_right h) +/-- **Every prime passes the AKS tests**, for every `r`: the easy direction of +[Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004]. -/ +theorem passes_of_prime {p : ℕ} (hp : p.Prime) (r : ℕ) : Passes p r := + passesUpTo_of_prime hp r (numTests p r) + end Cslib.AKS diff --git a/Complexity/Foundations/NumberTheory/AKS/IntegerVariant.lean b/Complexity/Foundations/NumberTheory/AKS/IntegerVariant.lean new file mode 100644 index 000000000..6281427c9 --- /dev/null +++ b/Complexity/Foundations/NumberTheory/AKS/IntegerVariant.lean @@ -0,0 +1,123 @@ +/- +Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Samuel Schlesinger +-/ + +module +-- upstream: candidate + +public import Complexity.Foundations.NumberTheory.AKS.Correctness +import Mathlib.Data.Nat.Log + +/-! +# The AKS tests with integer thresholds + +The AKS algorithm [Agrawal, Kayal, Saxena, Section 4][AKS2004] compares the multiplicative order +of `n` modulo `r` with `log₂² n` (step 2) and tests the congruence for `a ≤ ⌊√φ(r) log₂ n⌋` +(step 5). A program decides these comparisons more easily with the integer threshold +`c = ⌈log₂ n⌉` in place of `log₂ n`: step 2 asks for `ord_r(n) > c²`, and step 5 tests +`a ≤ ⌊√φ(r) · c⌋ = ⌊√(φ(r) c²)⌋`. This is a standard integer variant. It is still correct, by +Theorem 4.1 of the paper (`Cslib.AKS.prime_iff_passes`): `c² ≥ log₂² n`, so the order condition +is stronger than the paper's, and every prime passes the congruence for every `a`, so testing more +values of `a` loses nothing. + +`Cslib.AKS.PassesUpTo n r A` states the tests with the congruence checked for `1 ≤ a ≤ A`; the +paper's `Cslib.AKS.Passes n r` is its case `A = ⌊√φ(r) log₂ n⌋`. For `n ≥ 2`, `ord_r(n) > c²` and +`A ≥ ⌊√φ(r) log₂ n⌋`, it holds exactly when `n` is prime +(`Cslib.AKS.prime_iff_passesUpTo`). The order condition is decided by small powers: for `r ≥ 2`, +`ord_r(n) > L` exactly when `n` is coprime to `r` and no `n^k` with `1 ≤ k ≤ L` is `1` modulo `r` +(`Cslib.AKS.lt_orderOf_iff`). + +## Main definitions + +* `Cslib.AKS.numTestsNat`: the integer bound `⌊√(φ(r) ⌈log₂ n⌉²)⌋` of step 5. + +## Main results + +* `Cslib.AKS.logb_le_clog`: `log₂ n ≤ ⌈log₂ n⌉`. +* `Cslib.AKS.numTests_le_numTestsNat`: the integer bound tests at least the paper's values. +* `Cslib.AKS.lt_orderOf_iff`: the order condition by small powers. +* `Cslib.AKS.prime_iff_passesUpTo`: the tests with integer thresholds characterize the primes. + +## References + +* [M. Agrawal, N. Kayal, N. Saxena, *PRIMES is in P*][AKS2004] +-/ + +@[expose] public section + +namespace Cslib.AKS + +/-- The integer bound `⌊√(φ(r) c²)⌋ = ⌊√φ(r) · c⌋` of step 5, with `c = ⌈log₂ n⌉`. -/ +def numTestsNat (n r : ℕ) : ℕ := Nat.sqrt (Nat.totient r * Nat.clog 2 n ^ 2) + +/-- `log₂ n ≤ ⌈log₂ n⌉`. -/ +theorem logb_le_clog (n : ℕ) : Real.logb 2 n ≤ Nat.clog 2 n := by + rcases Nat.eq_zero_or_pos n with rfl | hn + · simp + have hn0 : (0 : ℝ) < n := by exact_mod_cast hn + rw [Real.logb_le_iff_le_rpow one_lt_two hn0, Real.rpow_natCast] + exact_mod_cast Nat.le_pow_clog one_lt_two n + +/-- The integer bound tests at least the paper's values: `⌊√φ(r) log₂ n⌋ ≤ ⌊√(φ(r) c²)⌋`. -/ +theorem numTests_le_numTestsNat (n r : ℕ) : numTests n r ≤ numTestsNat n r := by + rw [numTests, numTestsNat, ← Real.nat_floor_real_sqrt_eq_nat_sqrt] + refine Nat.floor_le_floor ?_ + rw [Nat.cast_mul, Nat.cast_pow, Real.sqrt_mul (Nat.cast_nonneg _), + Real.sqrt_sq (Nat.cast_nonneg _)] + rcases Nat.eq_zero_or_pos n with rfl | hn + · simp + have hlog : 0 ≤ Real.logb 2 n := + Real.logb_nonneg one_lt_two (by exact_mod_cast hn) + exact mul_le_mul_of_nonneg_left (logb_le_clog n) (Real.sqrt_nonneg _) + +/-- **The order condition by small powers**: for `r ≥ 2`, the multiplicative order of `n` modulo +`r` exceeds `L` exactly when `n` is coprime to `r` and no `n^k` with `1 ≤ k ≤ L` is `1` modulo +`r`. -/ +theorem lt_orderOf_iff {n r L : ℕ} (hr : 2 ≤ r) : + L < orderOf (n : ZMod r) ↔ n.Coprime r ∧ ∀ k, 1 ≤ k → k ≤ L → n ^ k % r ≠ 1 := by + have hpow (k : ℕ) : n ^ k % r = 1 ↔ (n : ZMod r) ^ k = 1 := by + rw [← Nat.cast_pow, ← Nat.cast_one (R := ZMod r), ZMod.natCast_eq_natCast_iff', + Nat.one_mod_eq_one.mpr (by omega)] + constructor + · intro h + refine ⟨coprime_of_orderOf_pos (by omega), fun k hk hkL hk1 => ?_⟩ + have := Nat.le_of_dvd (by omega) (orderOf_dvd_of_pow_eq_one ((hpow k).mp hk1)) + omega + · rintro ⟨hco, hsmall⟩ + have : NeZero r := ⟨by omega⟩ + have hpos : 0 < orderOf (n : ZMod r) := by + rw [← ZMod.coe_unitOfCoprime n hco, orderOf_units] + exact orderOf_pos _ + by_contra hle + exact hsmall _ hpos (by omega) ((hpow _).mpr (pow_orderOf_eq_one _)) + +/-- **The tests with integer thresholds characterize the primes**: for `n ≥ 2`, `r` such that the +order of `n` modulo `r` exceeds `⌈log₂ n⌉²`, and `A ≥ ⌊√φ(r) log₂ n⌋`, `n` passes the tests up to +`A` exactly when it is prime. This is [Agrawal, Kayal, Saxena, Theorem 4.1][AKS2004] with integer +thresholds. -/ +theorem prime_iff_passesUpTo {n r A : ℕ} (hn : 2 ≤ n) + (hr : Nat.clog 2 n ^ 2 < orderOf (n : ZMod r)) (hA : numTests n r ≤ A) : + n.Prime ↔ PassesUpTo n r A := by + refine ⟨fun hp => passesUpTo_of_prime hp r A, fun ⟨h1, h3, h5⟩ => ?_⟩ + have hlog : Real.logb 2 n ^ 2 < orderOf (n : ZMod r) := by + have h0 : 0 ≤ Real.logb 2 n := + Real.logb_nonneg one_lt_two (by exact_mod_cast (by omega : 1 ≤ n)) + have h1 : Real.logb 2 n ^ 2 ≤ (Nat.clog 2 n : ℝ) ^ 2 := + pow_le_pow_left₀ h0 (logb_le_clog n) 2 + have h2 : ((Nat.clog 2 n ^ 2 : ℕ) : ℝ) < orderOf (n : ZMod r) := by exact_mod_cast hr + push_cast at h2 + linarith + refine prime_of_passes hn hlog ⟨h1, h3, ?_⟩ + rcases h5 with h | h + · exact Or.inl h + · exact Or.inr fun a ha1 ha => h a ha1 (ha.trans hA) + +/-- The tests with the integer bound `⌊√(φ(r) ⌈log₂ n⌉²)⌋` characterize the primes. -/ +theorem prime_iff_passesUpTo_numTestsNat {n r : ℕ} (hn : 2 ≤ n) + (hr : Nat.clog 2 n ^ 2 < orderOf (n : ZMod r)) : + n.Prime ↔ PassesUpTo n r (numTestsNat n r) := + prime_iff_passesUpTo hn hr (numTests_le_numTestsNat n r) + +end Cslib.AKS diff --git a/ComplexityTests/AKSCriterion.lean b/ComplexityTests/AKSCriterion.lean index f9ebb874a..4bb4107a7 100644 --- a/ComplexityTests/AKSCriterion.lean +++ b/ComplexityTests/AKSCriterion.lean @@ -6,8 +6,8 @@ Authors: Samuel Schlesinger import Complexity.Foundations.NumberTheory.AKS.Criterion -/-! The AKS identity on a prime and a composite, the binomial divisibility, and the AKS tests on -a prime. -/ +/-! The AKS identity on a prime and a composite, the binomial divisibility, the AKS tests on a +prime, and `Passes` as `PassesUpTo` at the paper's bound. -/ namespace ComplexityTests.AKSCriterion @@ -27,4 +27,9 @@ example (r : ℕ) : passes_of_prime (by decide) r, fun h => h.1 ⟨2, 3, by norm_num, by norm_num⟩⟩ +-- `Passes` is `PassesUpTo` at the paper's bound `⌊√φ(r) log₂ n⌋`, and the prime `7` passes the +-- tests up to every bound. +example (n r A : ℕ) : (Passes n r ↔ PassesUpTo n r (numTests n r)) ∧ PassesUpTo 7 r A := + ⟨Iff.rfl, passesUpTo_of_prime (by decide) r A⟩ + end ComplexityTests.AKSCriterion diff --git a/ComplexityTests/WordRAMAKSSearch.lean b/ComplexityTests/WordRAMAKSSearch.lean new file mode 100644 index 000000000..2622be1ad --- /dev/null +++ b/ComplexityTests/WordRAMAKSSearch.lean @@ -0,0 +1,103 @@ +/- +Copyright (c) 2026 Samuel Schlesinger. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Samuel Schlesinger +-/ + +import Complexity.Algorithms.WordRAM.AKS.CongruenceLoop +import Complexity.Algorithms.WordRAM.AKS.GcdTotient +import Complexity.Algorithms.WordRAM.AKS.OrderSearch +import Complexity.Foundations.NumberTheory.AKS.IntegerVariant + +/-! +# The order search, step 3 and the integer thresholds of the AKS test + +Concrete runs of the search for `r`, of the step-3 and totient loop and of the step-5 loop, +instances of their correctness theorems, and checked instances of the integer-threshold +conditions. +-/ + +namespace ComplexityTests.WordRAMAKSSearch + +open Cslib.WordRAM Cslib.Algorithms.WordRAM + +/-! ## The search for `r` -/ + +private def searchInput (n L : ℕ) : Data 13 := + (Data.init []).set AKS.OrderSearch.n n |>.set AKS.OrderSearch.L L + +/-- For `n = 5` and `L = 2`: `r = 2` fails (`5 ≡ 1`), `r = 3` fails (`5² ≡ 1`), `r = 4` fails +(`5 ≡ 1`), and `r = 5` fails (`gcd = 5`); `r = 6` fails (`5² ≡ 1`), and `r = 7` is accepted. -/ +example : (Code.eval 16 16 AKS.OrderSearch.program (searchInput 5 2)).map + (·.1.regs AKS.OrderSearch.r) = some 7 := by + decide +kernel + +/-- `Good 5 2 7`: `5` is coprime to `7`, and `5 ≡ 5`, `5² ≡ 4` modulo `7`. -/ +theorem good_5_2_7 : AKS.OrderSearch.Good 5 2 7 := by + refine ⟨by decide, fun k hk1 hk2 => ?_⟩ + interval_cases k <;> decide + +/-- The correctness theorem with `r₀ = 7`: the search stops at a good `2 ≤ r ≤ 7`. -/ +example : ∃ y, ExecWithin 16 AKS.OrderSearch.program (searchInput 5 2) y + (AKS.OrderSearch.time 2 7) ∧ AKS.OrderSearch.Good 5 2 (y.regs AKS.OrderSearch.r) := by + obtain ⟨y, hy, -, -, hgood, -⟩ := AKS.OrderSearch.exec (w := 16) (searchInput 5 2) + (by decide) (by decide) (by decide) good_5_2_7 (by decide) (by decide) (by decide) + exact ⟨y, hy, hgood⟩ + +/-- `Good` is the order condition: `ord₇(5) = 6 > 2`. -/ +example : 2 < orderOf ((5 : ℕ) : ZMod 7) := + (Cslib.AKS.lt_orderOf_iff (by norm_num)).mpr good_5_2_7 + +/-! ## Step 3 and the totient -/ + +private def gcdInput (n r : ℕ) : Data 11 := + (Data.init []).set AKS.GcdTotient.n n |>.set AKS.GcdTotient.r r + +/-- For `n = 15` and `r = 4`, `gcd(3, 15) = 3` fails step 3; for `n = 7` it passes; `φ(4) = 2`. -/ +example : (Code.eval 16 8 AKS.GcdTotient.program (gcdInput 15 4)).map + (fun yn => (yn.1.regs AKS.GcdTotient.ok, yn.1.regs AKS.GcdTotient.phi)) = some (0, 2) ∧ + (Code.eval 16 8 AKS.GcdTotient.program (gcdInput 7 4)).map + (fun yn => (yn.1.regs AKS.GcdTotient.ok, yn.1.regs AKS.GcdTotient.phi)) = some (1, 2) := by + decide +kernel + +/-- The correctness theorem computes `φ(4) = 2`. -/ +example : ∃ y, ExecWithin 16 AKS.GcdTotient.program (gcdInput 7 4) y (AKS.GcdTotient.time 4) ∧ + y.regs AKS.GcdTotient.phi = 2 := by + obtain ⟨y, hy, -, -, hphi, -⟩ := AKS.GcdTotient.exec (w := 16) (gcdInput 7 4) (N := 7) (R := 4) + (by decide) + (by decide) (by decide) (by decide) (by decide) + exact ⟨y, hy, hphi.trans (by decide)⟩ + +/-! ## Step 5 -/ + +/-- Blocks at `0`, `r` and `2r`, with `r`, `n` and the last tested value `A` in their registers. -/ +private def loopInput (r n A : ℕ) : Data 29 := + ⟨fun j => if j = AKS.CongruenceLoop.B then r else if j = AKS.CongruenceLoop.T then 2 * r else + if j = AKS.CongruenceLoop.rr then r else if j = AKS.CongruenceLoop.nn then n else + if j = AKS.CongruenceLoop.A then A else 0, + fun _ => 0⟩ + +/-- Modulo `X - 1`, `4` fails at `a = 1` (`2⁴ ≡ 0 ≢ 2`), so the flag is cleared. -/ +example : (Code.eval 16 4 AKS.CongruenceLoop.program (loopInput 1 4 1)).map + (·.1.regs AKS.CongruenceLoop.ok) = some 0 := by + decide +kernel + +/-- The correctness theorem: for the prime `5`, every `a ≤ 3` passes modulo `X² - 1`. -/ +example : ∃ y, ExecWithin 16 AKS.CongruenceLoop.program (loopInput 2 5 3) y + (AKS.CongruenceLoop.time 2 5 3) ∧ y.regs AKS.CongruenceLoop.ok = 1 := by + obtain ⟨y, hy, -, hiff, -⟩ := AKS.CongruenceLoop.exec (w := 16) (loopInput 2 5 3) (a := 0) + (b := 2) (t := 4) (r := 2) (n := 5) (A₀ := 3) rfl rfl rfl rfl rfl rfl (by decide) + (by decide) (by decide) (by decide) (by decide) (by decide) (by decide) (by decide) + (by decide) (by decide) (by decide) + exact ⟨y, hy, hiff.mpr fun v _ _ => Cslib.AKS.congruence_of_prime (by decide) 2 v⟩ + +/-! ## Integer thresholds -/ + +/-- Step 3 rejects `15` with `r = 3`, since `gcd(3, 15) = 3`. -/ +example (A : ℕ) : ¬ Cslib.AKS.PassesUpTo 15 3 A := by + rintro ⟨-, h3, -⟩ + have := h3 3 (by norm_num) le_rfl + revert this + decide + +end ComplexityTests.WordRAMAKSSearch