Skip to content
  1. CLRS-Lean
  2. Chapters 1–35 · CLRS Fourth Edition
  3. Chapter 31. Number-Theoretic Algorithms
Browse chapters

Chapter 31 — Number-Theoretic Algorithms

CLRS, fourth edition · Lean 4 formalization

The proofs below use the models and assumptions described in the scope and implementation notes.

Imports
import Mathlib

31.1. Elementary Number-Theoretic Notions

CLRS §31.1: divisibility, the division theorem, the greatest common divisor, coprime, and prime numbers.

Main results:

  • Lemma 31.1 (divides_refl, divides_zero, divides_trans, divides_mul_right, divides_add, divides_sub): the basic divisibility facts used throughout the chapter.

  • Theorem 31.1 (division_theorem): for b > 0, a = q·b + r with 0 ≤ r < b has a unique quotient q and remainder r.

  • IsGCD / nat_gcd_isGCD / IsGCD.eq_gcd: the greatest-common-divisor property predicate and its agreement with Mathlib's Nat.gcd.

  • coprime_iff_gcd_eq_one / coprime_iff_no_common_divisor: characterizations of Nat.Coprime.

  • prime_def_gt_one, prime_two, and exists_prime_ge (Euclid's theorem: there are infinitely many primes).

  • lcm_dvd_left / lcm_dvd_right / lcm_dvd_of_dvd / lcm_dvd_iff / lcm_comm / lcm_assoc / lcm_eq_zero_iff: the basic least-common-multiple facts.

  • gcd_mul_lcm_eq: the identity gcd(a, b) · lcm(a, b) = a · b.

  • lcm_eq_mul_of_coprime: for coprime a b, lcm(a, b) = a · b.

Notation:

  • a ∣ b : a divides b.

  • Nat.gcd a b : the greatest common divisor.

  • Nat.lcm a b : the least common multiple.

  • Nat.Coprime a b : gcd a b = 1.

  • Nat.Prime p : p is prime.

Deferred: none.

namespace CLRSnamespace Chapter31

g is a greatest common divisor of a and b: a common divisor that is divisible by every common divisor (the universal property). For g > 0 this is equivalent to CLRS's wording "the largest d with d ∣ a and d ∣ b"; the universal form also behaves at g = 0.

def IsGCD (g a b : ℕ) : Prop := g ∣ a ∧ g ∣ b ∧ ∀ d : ℕ, d ∣ a → d ∣ b → d ∣ g

Nat.gcd a b is a greatest common divisor of a and b.

theorem nat_gcd_isGCD (a b : ℕ) : IsGCD (Nat.gcd a b) a b := ⟨Nat.gcd_dvd_left a b, Nat.gcd_dvd_right a b, fun Variable name `d` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`d hda hdb => Nat.dvd_gcd hda hdb⟩

The greatest-common-divisor property determines the value: if g is a greatest common divisor of a and b, then g = Nat.gcd a b.

theorem IsGCD.eq_gcd {g a b : ℕ} (hg : IsGCD g a b) : g = Nat.gcd a b := Nat.dvd_antisymm (Nat.dvd_gcd hg.1 hg.2.1) (hg.2.2 (Nat.gcd a b) (Nat.gcd_dvd_left a b) (Nat.gcd_dvd_right a b))

For a positive greatest common divisor, the universal property gives the CLRS "largest common divisor" form: every common divisor is ≤ g.

theorem IsGCD.greatest {g a b : ℕ} (hgpos : 0 < g) (hg : IsGCD g a b) : ∀ d : ℕ, d ∣ a → d ∣ b → d ≤ g := by intro d hda hdb exact Nat.le_of_dvd hgpos (hg.2.2 d hda hdb)

Lemma 31.1: every number divides itself.

theorem divides_refl (a : ℕ) : a ∣ a := Nat.dvd_refl a

Lemma 31.1: every number divides zero.

theorem divides_zero (a : ℕ) : a ∣ 0 := Nat.dvd_zero a

Lemma 31.1: divisibility is transitive.

theorem divides_trans {a b c : ℕ} (hab : a ∣ b) (hbc : b ∣ c) : a ∣ c := Nat.dvd_trans hab hbc

Lemma 31.1: if a ∣ b then a ∣ b·c.

theorem divides_mul_right {a b c : ℕ} (hab : a ∣ b) : a ∣ b * c := Nat.dvd_mul_right_of_dvd hab c

Lemma 31.1: if a ∣ b and a ∣ c then a ∣ b + c.

theorem divides_add {a b c : ℕ} (hab : a ∣ b) (hac : a ∣ c) : a ∣ b + c := Nat.dvd_add hab hac

Lemma 31.1: if a ∣ b and a ∣ c then a ∣ b - c.

theorem divides_sub {a b c : ℕ} (hab : a ∣ b) (hac : a ∣ c) : a ∣ b - c := Nat.dvd_sub hab hac

If a divides both b and c, then a divides every integer linear combination x·b + y·c of them.

theorem divides_linear_combination {a b c x y : ℕ} (hab : a ∣ b) (hac : a ∣ c) : a ∣ x * b + y * c := Nat.dvd_add (by simpa [Nat.mul_comm] using Nat.dvd_mul_right_of_dvd hab x) (by simpa [Nat.mul_comm] using Nat.dvd_mul_right_of_dvd hac y)

Two representations of a as q·b + r with 0 ≤ r < b agree: the division theorem's (q, r) is unique.

lemma division_unique (a b q₁ r₁ q₂ r₂ : ℕ) (hb : 0 < b) (h₁ : a = q₁ * b + r₁) (hr₁ : r₁ < b) (h₂ : a = q₂ * b + r₂) (hr₂ : r₂ < b) : q₁ = q₂ ∧ r₁ = r₂ := by have h : q₁ * b + r₁ = q₂ * b + r₂ := h₁.symm.trans h₂ have hmod₁ : (q₁ * b + r₁) % b = r₁ := by rw [Nat.add_mod, Nat.mul_mod] simp [Nat.mod_eq_of_lt hr₁] have hmod₂ : (q₂ * b + r₂) % b = r₂ := by rw [Nat.add_mod, Nat.mul_mod] simp [Nat.mod_eq_of_lt hr₂] have hr : r₁ = r₂ := by calc r₁ = (q₁ * b + r₁) % b := hmod₁.symm _ = (q₂ * b + r₂) % b := congrArg (fun x => x % b) h _ = r₂ := hmod₂ have hb_mul : q₁ * b = q₂ * b := by exact Nat.add_right_cancel (by rwa [← hr] at h) have hq : q₁ = q₂ := Nat.mul_right_cancel hb hb_mul exact ⟨hq, hr⟩

Theorem 31.1 (Division theorem). For a and b > 0, there is a unique pair (q, r) with a = q·b + r and 0 ≤ r < b (over ℕ, 0 ≤ r is automatic). q is the quotient and r the remainder.

theorem division_theorem (a b : ℕ) (hb : 0 < b) : ∃! q : ℕ, ∃! r : ℕ, a = q * b + r ∧ r < b := by have hda : a = a / b * b + a % b := by rw [Nat.mul_comm] exact (Nat.div_add_mod a b).symm refine ⟨a / b, ?_⟩ constructor · refine ⟨a % b, ?_, ?_⟩ · exact ⟨hda, Nat.mod_lt a hb⟩ · intro r' hr' exact (division_unique a b (a / b) (a % b) (a / b) r' hb hda (Nat.mod_lt a hb) hr'.1 hr'.2).2.symm · intro q hq rcases hq with ⟨r, hr, hrq⟩ exact (division_unique a b q r (a / b) (a % b) hb hr.1 hr.2 hda (Nat.mod_lt a hb)).1

a and b are coprime exactly when their gcd is one.

theorem coprime_iff_gcd_eq_one (a b : ℕ) : Nat.Coprime a b ↔ Nat.gcd a b = 1 := Nat.coprime_iff_gcd_eq_one

a and b are coprime exactly when they have no common divisor other than one.

theorem coprime_iff_no_common_divisor (a b : ℕ) : Nat.Coprime a b ↔ ∀ d : ℕ, d ∣ a → d ∣ b → d = 1 := by rw [coprime_iff_gcd_eq_one, Nat.gcd_eq_one_iff]

CLRS's definition of prime: p > 1 and the only divisors of p are 1 and p itself.

theorem prime_def_gt_one (p : ℕ) : Nat.Prime p ↔ 1 < p ∧ ∀ m : ℕ, m ∣ p → m = 1 ∨ m = p := by rw [Nat.prime_def] constructor · intro hp exact ⟨Nat.lt_of_lt_of_le (by norm_num) hp.1, hp.2⟩ · intro hp exact ⟨by omega, hp.2⟩

2 is prime.

theorem prime_two : Nat.Prime 2 := by decide

Euclid's theorem: for every n there is a prime ≥ n, so there are infinitely many primes.

theorem exists_prime_ge (n : ℕ) : ∃ p : ℕ, Nat.Prime p ∧ n ≤ p := by rcases Nat.exists_infinite_primes n with ⟨p, hp, hprime⟩ exact ⟨p, hprime, hp⟩

Least common multiple

The least common multiple is a multiple of a.

theorem lcm_dvd_left (a b : ℕ) : a ∣ Nat.lcm a b := dvd_lcm_left a b

The least common multiple is a multiple of b.

theorem lcm_dvd_right (a b : ℕ) : b ∣ Nat.lcm a b := dvd_lcm_right a b

Any common multiple of a and b is a multiple of their least common multiple.

theorem lcm_dvd_of_dvd {a b c : ℕ} (ha : a ∣ c) (hb : b ∣ c) : Nat.lcm a b ∣ c := lcm_dvd ha hb

Nat.lcm a b divides c exactly when both a and b divide c.

theorem lcm_dvd_iff {a b c : ℕ} : Nat.lcm a b ∣ c ↔ a ∣ c ∧ b ∣ c := _root_.lcm_dvd_iff

The least common multiple is commutative.

theorem lcm_comm (a b : ℕ) : Nat.lcm a b = Nat.lcm b a := _root_.lcm_comm a b

The least common multiple is associative.

theorem lcm_assoc (a b c : ℕ) : Nat.lcm (Nat.lcm a b) c = Nat.lcm a (Nat.lcm b c) := _root_.lcm_assoc a b c

The least common multiple is zero exactly when one of the arguments is zero.

theorem lcm_eq_zero_iff (a b : ℕ) : Nat.lcm a b = 0 ↔ a = 0 ∨ b = 0 := _root_.lcm_eq_zero_iff a b

CLRS identity: gcd(a, b) · lcm(a, b) = a · b.

theorem gcd_mul_lcm_eq (a b : ℕ) : Nat.gcd a b * Nat.lcm a b = a * b := associated_iff_eq.mp (_root_.gcd_mul_lcm a b)

For coprime a and b, the least common multiple is the product a · b.

theorem lcm_eq_mul_of_coprime {a b : ℕ} (h : Nat.Coprime a b) : Nat.lcm a b = a * b := h.lcm_eq_mul
end Chapter31end CLRS
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_1_Elementary_Number_Theory

31.2. Greatest Common Divisor

CLRS §31.2: the Euclid recursion, Bezout's lemma, the smallest-positive-linear-combination characterization of the gcd, and the EUCLID and EXTENDED-EUCLID algorithms.

Main results:

  • Lemma 31.2 (euclid_recursion, gcd_zero_left, gcd_zero_right): gcd(a, b) = gcd(b mod a, a) and the base cases gcd(0, b) = b, gcd(a, 0) = a.

  • euclid + euclid_eq_gcd + euclid_terminates: the EUCLID algorithm is a total function and always returns Nat.gcd a b.

  • Running time (Lamé / Fibonacci): euclidDivisions counts the recursive calls of EUCLID. Lemma 31.10 (fib_le_of_euclidDivisions) gives the Fibonacci lower bounds a ≥ F_{k+2} and b ≥ F_{k+1} for k calls; Theorem 31.11, Lamé's theorem (euclidDivisions_lt), is the running-time bound b < F_{k+1} ⇒ fewer than k calls; and Corollary 31.12 (euclidDivisions_le_two_log) records the O(log b) bound. The helper lemmas fib_two_step_ge_pow_two and pow_two_le_fib prove the exponential Fibonacci growth 2^(n/2) ≤ F_{n+2} used by Corollary 31.12.

  • Lemma 31.3 (gcd_is_linear_combination): Bezout's identity — gcd a b is an integer linear combination of a and b.

  • Theorem 31.2 (gcd_is_smallest_positive_linear_combination): gcd a b is the smallest positive linear combination of a and b (when a ≠ 0 ∨ b ≠ 0).

  • Corollary 31.3 (gcd_dvd_linear_combination): gcd a b divides every linear combination of a and b.

  • Corollary 31.4 (gcd_eq_one_iff_coprime, coprime_iff_one_linear_combination, gcd_div_gcd_coprime): coprime characterizations, including coprime (a/g) (b/g) for g = gcd a b.

  • extendedEuclid + extendedEuclid_spec: EXTENDED-EUCLID returns (d, x, y) with d = gcd a b = a·x + b·y.

Notation:

  • Nat.gcd a b : the greatest common divisor.

  • Nat.gcdA a b / Nat.gcdB a b : the Bezout coefficients.

  • x·a + y·b : integer linear combinations (coefficients in ℤ).

This section is complete. Modular arithmetic, primality testing, and RSA are covered in their own sections (§31.3–31.9).

open Natnamespace CLRSnamespace Chapter31

Lemma 31.2 (Euclid's recursion): gcd(a, b) = gcd(b mod a, a).

theorem euclid_recursion (a b : ℕ) : Nat.gcd a b = Nat.gcd (b % a) a := Nat.gcd_rec a b

Lemma 31.2: gcd(0, b) = b.

theorem gcd_zero_left (b : ℕ) : Nat.gcd 0 b = b := Nat.gcd_zero_left b

Lemma 31.2: gcd(a, 0) = a.

theorem gcd_zero_right (a : ℕ) : Nat.gcd a 0 = a := Nat.gcd_zero_right a

EUCLID (CLRS §31.2). Compute gcd a b by the recursion gcd(a, b) = gcd(b mod a, a) for a > 0 and gcd(0, b) = b. Total by well-founded recursion on the first argument (it strictly decreases to b mod (a+1) < a+1), so termination is part of the definition.

def euclid : ℕ → ℕ → ℕ | 0, b => b | a' + 1, b => euclid (b % (a' + 1)) (a' + 1) termination_by a Variable name `b` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`b => a decreasing_by exact Nat.mod_lt b (Nat.succ_pos a')

EUCLID is correct: euclid a b returns Nat.gcd a b.

theorem euclid_eq_gcd (a b : ℕ) : euclid a b = Nat.gcd a b := by induction a, b using Nat.gcd.induction with | H0 b => simp [euclid] | H1 a b ha ih => cases a with | zero => exact (Nat.ne_of_gt ha rfl).elim | succ a' => simp [euclid] rw [ih] rw [← Nat.gcd_succ]

EUCLID terminates: it is total, so every call returns a value.

theorem euclid_terminates (a b : ℕ) : ∃ r : ℕ, euclid a b = r := ⟨euclid a b, rfl⟩

EUCLID division count (CLRS §31.2). The number of recursive calls that EUCLID(a, b) makes under CLRS's recursion EUCLID(a, b) = EUCLID(b, a mod b) for b > 0, with EUCLID(a, 0) = a. Total by well-founded recursion on the second argument (each call drops to a mod b < b).

def euclidDivisions : ℕ → ℕ → ℕ | _, 0 => 0 | a, b + 1 => 1 + euclidDivisions (b + 1) (a % (b + 1)) termination_by _ b => b decreasing_by exact Nat.mod_lt a (Nat.succ_pos b)

Lemma 31.10 (Lamé, core direction). If EUCLID(a, b) with a > b ≥ 1 invokes k recursive calls, then the Fibonacci bounds b ≥ F_{k+1} and a ≥ F_{k+2} hold.

theorem fib_le_of_euclidDivisions (a b : ℕ) (hb0 : 0 < b) (hba : b < a) : fib (euclidDivisions a b + 1) ≤ b ∧ fib (euclidDivisions a b + 2) ≤ a := by revert a hb0 hba induction b using Nat.strong_induction_on with | h b ih => intro a hb0 hba cases b with | zero => omega | succ b' => let r := a % (b' + 1) have hmain : euclidDivisions a (b' + 1) = 1 + euclidDivisions (b' + 1) r := by simp [euclidDivisions, r] by_cases hr0 : r = 0 · have hk1 : euclidDivisions a (b' + 1) = 1 := by rw [hmain] simp [euclidDivisions, r, hr0] constructor · rw [hk1] simp [fib_two] · rw [hk1] simp [fib_add_two] have hdvd : b' + 1 ∣ a := by exact Nat.dvd_of_mod_eq_zero (by simpa [r] using hr0) rcases hdvd with ⟨q, ha⟩ have hqb : b' + 1 < (b' + 1) * q := by simpa [ha] using hba have hq : 1 < q := by exact (Nat.mul_lt_mul_left (by omega : 0 < b' + 1)).mp (by simpa using hqb) have hle : 2 ≤ (b' + 1) * q := by simpa using (Nat.mul_le_mul (by omega : 1 ≤ b' + 1) (Nat.succ_le_of_lt hq)) rw [ha] exact hle · have hrpos : 0 < r := Nat.pos_of_ne_zero hr0 have hrlt : r < b' + 1 := by simpa [r] using (Nat.mod_lt a (Nat.succ_pos b')) have ih' := ih r hrlt (b' + 1) hrpos hrlt let k' := euclidDivisions (b' + 1) r have hk' : euclidDivisions a (b' + 1) = 1 + k' := by simpa [k'] using hmain constructor · rw [hk'] have harg : (1 + k') + 1 = k' + 2 := by omega rw [harg] simpa [k'] using ih'.2 · rw [hk'] have hfib3 : fib (k' + 3) = fib (k' + 1) + fib (k' + 2) := by simpa [Nat.add_assoc] using (fib_add_two (n := k' + 1)) have hfib_sum : fib (k' + 1) + fib (k' + 2) ≤ r + (b' + 1) := by exact Nat.add_le_add (by simpa [k'] using ih'.1) (by simpa [k'] using ih'.2) have hra : r + (b' + 1) ≤ a := by have hmod : a = a / (b' + 1) * (b' + 1) + r := by change a = a / (b' + 1) * (b' + 1) + a % (b' + 1) rw [mul_comm] exact (Nat.div_add_mod a (b' + 1)).symm have hq1 : 1 ≤ a / (b' + 1) := by rw [Nat.le_div_iff_mul_le (Nat.succ_pos b')] simpa using (Nat.le_of_lt hba) rw [hmod] have hle : b' + 1 ≤ a / (b' + 1) * (b' + 1) := by simpa using (Nat.mul_le_mul_right (b' + 1) hq1) omega have harg : (1 + k') + 2 = k' + 3 := by omega calc fib ((1 + k') + 2) = fib (k' + 3) := by rw [harg] _ = fib (k' + 1) + fib (k' + 2) := hfib3 _ ≤ r + (b' + 1) := hfib_sum _ ≤ a := hra

Theorem 31.11 (Lamé's theorem). For k ≥ 1, if a > b ≥ 1 and b < F_{k+1}, then EUCLID(a, b) makes fewer than k recursive calls — equivalently, EUCLID needs O(log b) calls for inputs a > b.

theorem euclidDivisions_lt {a b k : ℕ} (_hk : 1 ≤ k) (hb0 : 0 < b) (hba : b < a) (hbf : b < fib (k + 1)) : euclidDivisions a b < k := by have hcore := fib_le_of_euclidDivisions a b hb0 hba by_contra hnot have hge : k ≤ euclidDivisions a b := Nat.le_of_not_gt hnot have hmono : fib (k + 1) ≤ fib (euclidDivisions a b + 1) := by exact fib_mono (by omega) have : fib (k + 1) ≤ b := le_trans hmono hcore.1 exact (not_lt_of_ge this) hbf

Exponential Fibonacci growth. For every t, fib(2t+2) ≥ 2^t and fib(2t+3) ≥ 2^t; hence the Fibonacci sequence grows at least like 2^(n/2). This is the exponential bound that turns Lamé's theorem into the O(log b) running-time bound.

theorem fib_two_step_ge_pow_two (t : ℕ) : (2 ^ t ≤ fib (2 * t + 2)) ∧ (2 ^ t ≤ fib (2 * t + 3)) := by induction t with | zero => constructor <;> norm_num [fib_two, fib_add_two, fib_one] | succ t ih => have hA0 : 2 ^ t ≤ fib (2 * t + 2) := ih.1 have hB0 : 2 ^ t ≤ fib (2 * t + 3) := ih.2 have hA1 : 2 ^ (t + 1) ≤ fib (2 * (t + 1) + 2) := by have hfib : fib (2 * (t + 1) + 2) = fib (2 * t + 2) + fib (2 * t + 3) := by rw [show 2 * (t + 1) + 2 = 2 * t + 4 by ring] have h := fib_add_two (n := 2 * t + 2) rw [show (2 * t + 2) + 2 = 2 * t + 4 by omega, show (2 * t + 2) + 1 = 2 * t + 3 by omega] at h exact h rw [hfib] rw [show 2 ^ (t + 1) = 2 ^ t + 2 ^ t by rw [pow_succ]; ring] exact Nat.add_le_add hA0 hB0 have hB1 : 2 ^ (t + 1) ≤ fib (2 * (t + 1) + 3) := by have hfib : fib (2 * (t + 1) + 3) = fib (2 * t + 3) + fib (2 * t + 4) := by rw [show 2 * (t + 1) + 3 = 2 * t + 5 by ring] have h := fib_add_two (n := 2 * t + 3) rw [show (2 * t + 3) + 2 = 2 * t + 5 by omega, show (2 * t + 3) + 1 = 2 * t + 4 by omega] at h exact h rw [hfib] have hA1' : 2 ^ (t + 1) ≤ fib (2 * t + 4) := by simpa [show 2 * (t + 1) + 2 = 2 * t + 4 by ring] using hA1 exact le_trans (Nat.le_add_left _ _) (by simpa [Nat.add_comm] using (Nat.add_le_add hB0 hA1')) exact ⟨hA1, hB1⟩

The Fibonacci sequence grows exponentially: 2^(n/2) ≤ fib (n+2) for all n.

theorem pow_two_le_fib (n : ℕ) : 2 ^ (n / 2) ≤ fib (n + 2) := by rcases Nat.even_or_odd n with ⟨t, rfl⟩ | ⟨t, rfl⟩ · have h2 : t + t = 2 * t := by omega simpa [h2] using (fib_two_step_ge_pow_two t).1 · have hdiv : (2 * t + 1) / 2 = t := by simpa [show 1 / 2 = 0 by norm_num] using (Nat.mul_add_div (by decide : 2 > 0) t 1) rw [hdiv] have harg : (2 * t + 1) + 2 = 2 * t + 3 := by omega rw [harg] exact (fib_two_step_ge_pow_two t).2

Corollary 31.12. For a > b ≥ 1, EUCLID(a, b) makes at most 2·log₂ b + 2 recursive calls — i.e. O(log b). Combining Lemma 31.10 with b ≥ F_{k+1} ≥ 2^{(k−1)/2} bounds the division count logarithmically.

theorem euclidDivisions_le_two_log (a b : ℕ) (hb0 : 0 < b) (hba : b < a) : euclidDivisions a b ≤ 2 * Nat.log 2 b + 2 := by have hkge1 : 1 ≤ euclidDivisions a b := by rcases b with _ | b' · omega · simp [euclidDivisions] have hcore := fib_le_of_euclidDivisions a b hb0 hba have hpow_le_fib : 2 ^ ((euclidDivisions a b - 1) / 2) ≤ fib (euclidDivisions a b + 1) := by have h := pow_two_le_fib (euclidDivisions a b - 1) rwa [show (euclidDivisions a b - 1) + 2 = euclidDivisions a b + 1 by omega] at h have hpow_le_b : 2 ^ ((euclidDivisions a b - 1) / 2) ≤ b := le_trans hpow_le_fib hcore.1 have hlog : (euclidDivisions a b - 1) / 2 ≤ Nat.log 2 b := Nat.le_log_of_pow_le (by decide : 1 < 2) hpow_le_b have hk1 : euclidDivisions a b - 1 ≤ 2 * Nat.log 2 b + 1 := by rw [Nat.div_le_iff_le_mul_add_pred (by decide : 0 < 2)] at hlog exact hlog omega

Lemma 31.3 (Bezout's identity): gcd a b is an integer linear combination of a and b.

theorem gcd_is_linear_combination (a b : ℕ) : (Nat.gcd a b : ℤ) = (a : ℤ) * Nat.gcdA a b + (b : ℤ) * Nat.gcdB a b := Nat.gcd_eq_gcd_ab a b

Corollary 31.3: gcd a b divides every integer linear combination of a and b.

theorem gcd_dvd_linear_combination (a b : ℕ) (x y : ℤ) : (Nat.gcd a b : ℤ) ∣ x * (a : ℤ) + y * (b : ℤ) := by have hga : (Nat.gcd a b : ℤ) ∣ (a : ℤ) := by exact_mod_cast Nat.gcd_dvd_left a b have hgb : (Nat.gcd a b : ℤ) ∣ (b : ℤ) := by exact_mod_cast Nat.gcd_dvd_right a b exact dvd_add (dvd_mul_of_dvd_right hga x) (dvd_mul_of_dvd_right hgb y)

Theorem 31.2, lower-bound part: gcd a b is no larger than any positive linear combination of a and b.

theorem gcd_le_positive_linear_combination (a b : ℕ) {z : ℤ} (hzpos : 0 < z) (hz : ∃ x y : ℤ, z = x * (a : ℤ) + y * (b : ℤ)) : (Nat.gcd a b : ℤ) ≤ z := by rcases hz with ⟨x, y, rfl⟩ exact Int.le_of_dvd hzpos (gcd_dvd_linear_combination a b x y)

Theorem 31.2. gcd a b is the smallest positive linear combination of a and b (for a ≠ 0 ∨ b ≠ 0). Bezout's identity shows it is itself a positive linear combination; Corollary 31.3 shows it divides (hence is ≤) every positive one.

theorem gcd_is_smallest_positive_linear_combination (a b : ℕ) (hab : a ≠ 0 ∨ b ≠ 0) : IsLeast {z : ℤ | 0 < z ∧ ∃ x y : ℤ, z = x * (a : ℤ) + y * (b : ℤ)} (Nat.gcd a b : ℤ) := by constructor · constructor · rcases hab with ha | hb · exact_mod_cast Nat.gcd_pos_of_pos_left b (Nat.pos_of_ne_zero ha) · exact_mod_cast Nat.gcd_pos_of_pos_right a (Nat.pos_of_ne_zero hb) · refine ⟨Nat.gcdA a b, Nat.gcdB a b, ?_⟩ rw [← mul_comm (a : ℤ), ← mul_comm (b : ℤ)] exact gcd_is_linear_combination a b · intro z hz exact gcd_le_positive_linear_combination a b hz.1 hz.2

Corollary 31.4: a and b are coprime exactly when gcd a b = 1.

theorem gcd_eq_one_iff_coprime (a b : ℕ) : Nat.gcd a b = 1 ↔ Nat.Coprime a b := Nat.coprime_iff_gcd_eq_one.symm

Corollary 31.4: a and b are coprime exactly when 1 is an integer linear combination of them.

theorem coprime_iff_one_linear_combination (a b : ℕ) : Nat.Coprime a b ↔ ∃ x y : ℤ, 1 = x * (a : ℤ) + y * (b : ℤ) := by constructor · intro hcop have h := Nat.gcd_eq_gcd_ab a b rw [hcop.gcd_eq_one] at h refine ⟨Nat.gcdA a b, Nat.gcdB a b, ?_⟩ rw [← mul_comm (a : ℤ), ← mul_comm (b : ℤ)] exact h · rintro ⟨x, y, h⟩ have hdiv : (Nat.gcd a b : ℤ) ∣ 1 := by rw [h] exact gcd_dvd_linear_combination a b x y have hgcd : Nat.gcd a b = 1 := by have hdiv' : (Nat.gcd a b : ℤ).natAbs ∣ (1 : ℤ).natAbs := Int.natAbs_dvd_natAbs.mpr hdiv have hg1 : (Nat.gcd a b : ℤ).natAbs = 1 := Nat.dvd_one.mp (by simpa using hdiv') exact (by simpa using hg1) exact hgcd

Corollary 31.4: dividing a and b by their gcd yields coprime numbers.

theorem gcd_div_gcd_coprime (a b : ℕ) (h : 0 < Nat.gcd a b) : Nat.Coprime (a / Nat.gcd a b) (b / Nat.gcd a b) := Nat.coprime_div_gcd_div_gcd h

EXTENDED-EUCLID (CLRS §31.2). Return (d, x, y) with d = gcd a b and d = a·x + b·y. Mathlib's Nat.gcdA/Nat.gcdB supply the Bezout coefficients.

def extendedEuclid (a b : ℕ) : ℕ × ℤ × ℤ := (Nat.gcd a b, Nat.gcdA a b, Nat.gcdB a b)

EXTENDED-EUCLID is correct: it returns d = gcd a b together with coefficients x, y satisfying d = a·x + b·y.

theorem extendedEuclid_spec (a b : ℕ) : let (d, x, y) := extendedEuclid a b d = Nat.gcd a b ∧ (d : ℤ) = (a : ℤ) * x + (b : ℤ) * y := by simp [extendedEuclid, Nat.gcd_eq_gcd_ab]
end Chapter31end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_31.Section_31_2_Greatest_Common_Divisor.Execution

Euclid follows the public first-argument recursion. Its returned counter charges one remainder/division step at each nonterminal call. The value refines the public function and the count equals the legacy division recurrence with swapped arguments; no bit-cost model is claimed.

namespace CLRS.Chapter31def euclidWithCount : Nat → Nat → Nat × Nat | 0, b => (b, 0) | a + 1, b => let next := euclidWithCount (b % (a + 1)) (a + 1) (next.1, next.2 + 1) termination_by a _ => a decreasing_by exact Nat.mod_lt _ (Nat.succ_pos _) theorem euclidWithCount_spec (a b : Nat) : euclidWithCount a b = (euclid a b, euclidDivisions b a) := by induction a using Nat.strong_induction_on generalizing b with | h a ih => cases a with | zero => simp [euclidWithCount, euclid, euclidDivisions] | succ a => rw [euclidWithCount, ih _ (Nat.mod_lt _ (Nat.succ_pos _)), euclid, euclidDivisions] simp [Nat.add_comm] @[simp] theorem euclidWithCount_value (a b : Nat) : (euclidWithCount a b).1 = euclid a b := by rw [euclidWithCount_spec] @[simp] theorem euclidWithCount_count (a b : Nat) : (euclidWithCount a b).2 = euclidDivisions b a := by rw [euclidWithCount_spec]end CLRS.Chapter31
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_2_Greatest_Common_Divisor

31.3. Modular Arithmetic

CLRS §31.3: the arithmetic of the ring Z_n of residues modulo n — the well-definedness of addition and multiplication, additive and multiplicative inverses, cancellation, and the solvability criterion for linear congruences.

Main results:

  • Theorem 31.5 (mod_add, mod_mul): (a + b) mod n and (a · b) mod n are well-defined on residues, so + and · descend to Z_n.

  • Theorem 31.6 (exists_mul_inverse_mod): if gcd(a, n) = 1, then a has a multiplicative inverse modulo n.

  • Theorem 31.9 (mul_left_cancel_mod): gcd(c, n) = 1 and a·c ≡ b·c (mod n) imply a ≡ b (mod n) (cancellation in Z_n).

  • Theorem 31.11 (modular_linear_solvable): the congruence a·x ≡ b (mod n) has a solution exactly when gcd(a, n) ∣ b.

Notation:

  • a ≡ b [MOD n] : Nat.ModEq — a and b leave the same remainder modulo n.

  • ZMod n : the ring of residues modulo n.

  • Nat.gcd a n : the greatest common divisor.

Deferred: none (the enumeration of all solutions of a linear congruence is proved in §31.4 via CLRS.Chapter31.modularLinearEquationSolver).

namespace CLRSnamespace Chapter31

Theorem 31.5: addition modulo n is well-defined on residues.

theorem mod_add (a b n : ℕ) : (a + b) % n = (a % n + b % n) % n := Nat.add_mod a b n

Theorem 31.5: multiplication modulo n is well-defined on residues.

theorem mod_mul (a b n : ℕ) : (a * b) % n = ((a % n) * (b % n)) % n := Nat.mul_mod a b n

Congruence preserves addition: a ≡ b and c ≡ d imply a + c ≡ b + d.

theorem modEq_add {a b c d n : ℕ} (hab : a ≡ b [MOD n]) (hcd : c ≡ d [MOD n]) : a + c ≡ b + d [MOD n] := Nat.ModEq.add hab hcd

Congruence preserves multiplication: a ≡ b and c ≡ d imply a·c ≡ b·d.

theorem modEq_mul {a b c d n : ℕ} (hab : a ≡ b [MOD n]) (hcd : c ≡ d [MOD n]) : a * c ≡ b * d [MOD n] := Nat.ModEq.mul hab hcd

Congruence preserves powers: a ≡ b implies aᵏ ≡ bᵏ.

theorem modEq_pow {a b n : ℕ} (hab : a ≡ b [MOD n]) (k : ℕ) : a ^ k ≡ b ^ k [MOD n] := Nat.ModEq.pow k hab

Congruence is an equivalence relation.

theorem modEq_refl (a n : ℕ) : a ≡ a [MOD n] := Nat.ModEq.refl a
theorem modEq_symm {a b n : ℕ} (h : a ≡ b [MOD n]) : b ≡ a [MOD n] := h.symmtheorem modEq_trans {a b c n : ℕ} (hab : a ≡ b [MOD n]) (hbc : b ≡ c [MOD n]) : a ≡ c [MOD n] := hab.trans hbc

Theorem 31.6. If gcd(a, n) = 1, then a has a multiplicative inverse modulo n: there is x with a·x ≡ 1 (mod n).

try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false` theorem exists_mul_inverse_mod {a n : ℕ} [NeZero n] (hcop : Nat.Coprime a n) : ∃ x : ℕ, a * x ≡ 1 [MOD n] := by let x : ℕ := ((a : ZMod n)⁻¹).val refine ⟨x, ?_⟩ rw [← ZMod.natCast_eq_natCast_iff] have hx : (x : ZMod n) = (a : ZMod n)⁻¹ := by dsimp [x] try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa using (ZMod.natCast_val (R := ZMod n) ((a : ZMod n)⁻¹)) rw [Nat.cast_mul, hx] simpa using (ZMod.coe_mul_inv_eq_one a hcop)

Theorem 31.9 (cancellation). If gcd(c, n) = 1 and a·c ≡ b·c (mod n), then a ≡ b (mod n).

theorem mul_left_cancel_mod {a b c n : ℕ} (hcop : Nat.Coprime c n) (h : a * c ≡ b * c [MOD n]) : a ≡ b [MOD n] := by have h' : c * a ≡ c * b [MOD n] := by simpa [Nat.mul_comm] using h exact Nat.ModEq.cancel_left_of_coprime (Nat.gcd_comm n c ▸ hcop) h'

Theorem 31.11. The linear congruence a·x ≡ b (mod n) has a solution exactly when gcd(a, n) ∣ b.

theorem modular_linear_solvable (a b n : ℕ) [NeZero n] : (∃ x : ℕ, a * x ≡ b [MOD n]) ↔ Nat.gcd a n ∣ b := by constructor · rintro ⟨x, h⟩ have hzn : (a * x : ZMod n) = (b : ZMod n) := by simpa using ((ZMod.natCast_eq_natCast_iff (a * x) b n).mpr h) let f : ZMod n →+* ZMod (Nat.gcd a n) := ZMod.castHom (Nat.gcd_dvd_right a n) (ZMod (Nat.gcd a n)) have hf : f (a * x : ZMod n) = f (b : ZMod n) := congrArg f hzn have ha0 : (a : ZMod (Nat.gcd a n)) = 0 := (ZMod.natCast_eq_zero_iff a (Nat.gcd a n)).mpr (Nat.gcd_dvd_left a n) have hb0 : (b : ZMod (Nat.gcd a n)) = 0 := by have hfax : f (a * x : ZMod n) = 0 := by rw [map_mul] simp [ha0] simpa using (hf.symm.trans hfax) exact (ZMod.natCast_eq_zero_iff b (Nat.gcd a n)).mp hb0 · intro h rcases h with ⟨k, hk⟩ have hbez := Nat.gcd_eq_gcd_ab a n let x0 : ℤ := Nat.gcdA a n * (k : ℤ) let y0 : ℤ := Nat.gcdB a n * (k : ℤ) have hb : (b : ℤ) = (a : ℤ) * x0 + (n : ℤ) * y0 := by calc (b : ℤ) = (Nat.gcd a n : ℤ) * (k : ℤ) := by rw [hk]; norm_num _ = ((a : ℤ) * Nat.gcdA a n + (n : ℤ) * Nat.gcdB a n) * (k : ℤ) := by rw [hbez] _ = (a : ℤ) * (Nat.gcdA a n * (k : ℤ)) + (n : ℤ) * (Nat.gcdB a n * (k : ℤ)) := by ring _ = (a : ℤ) * x0 + (n : ℤ) * y0 := by rfl have hcong : (a : ℤ) * x0 ≡ (b : ℤ) [ZMOD n] := by rw [hb] exact (Int.modEq_iff_dvd.2 ⟨y0, by ring⟩) let x : ℕ := Int.toNat (x0 % n) have hxrep : (x : ℤ) ≡ x0 [ZMOD n] := by dsimp [x] rw [Int.ModEq] have hnn : (n : ℤ) ≠ 0 := by exact_mod_cast (NeZero.ne n) have hnonneg : 0 ≤ x0 % n := Int.emod_nonneg x0 hnn rw [Int.toNat_of_nonneg hnonneg] rw [Int.emod_emod] have hcongx : (a : ℤ) * (x : ℤ) ≡ (b : ℤ) [ZMOD n] := (Int.ModEq.mul_left (a : ℤ) hxrep).trans hcong exact ⟨x, by exact_mod_cast hcongx⟩
end Chapter31end CLRS
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_3_Modular_Arithmetic

31.4. Solving Modular Linear Equations

CLRS §31.4: the structure of the solutions to the linear congruence a·x ≡ b (mod n). Once one solution x₀ is known, every solution is x₀ + k·(n/d) for d = gcd(a, n), so the d distinct solutions modulo n are x₀, x₀ + n/d, …, x₀ + (d−1)·(n/d).

Main results:

  • linear_congruence_shift: if x solves a·x ≡ b (mod n), then x + k·(n/d) solves it too — shifting by n/d preserves solutions.

  • linear_congruence_all_solutions: if x₀ and x both solve a·x ≡ b (mod n), then x ≡ x₀ (mod n/d) — every solution differs from x₀ by a multiple of n/d.

  • Theorem linear_congruence_solutions (Theorem 31.10): the solutions are exactly the residue class x₀ mod (n/d).

  • Theorem linear_congruence_distinct (Theorem 31.10): the d values k·(n/d) for 0 ≤ k < d are pairwise incongruent, so the congruence has exactly d distinct solutions.

Notation:

  • a ≡ b [MOD n] : Nat.ModEq.

  • Nat.gcd a n : the greatest common divisor.

Deferred: none (the executable enumerator modularLinearEquationSolver and its d-solution length bound are proved).

namespace CLRSnamespace Chapter31

a·(n/d) = (a/d)·n for d = gcd a n: the cross-term used in the shift.

lemma mul_nat_div_eq (a n : ℕ) [NeZero n] : a * (n / Nat.gcd a n) = (a / Nat.gcd a n) * n := by let d := Nat.gcd a n have hd : 0 < d := Nat.gcd_pos_of_pos_right a (NeZero.pos n) apply Nat.mul_right_cancel hd calc (a * (n / d)) * d = a * ((n / d) * d) := by ring _ = a * n := by rw [Nat.div_mul_cancel (Nat.gcd_dvd_right a n)] _ = ((a / d) * d) * n := by rw [Nat.div_mul_cancel (Nat.gcd_dvd_left a n)] _ = (a / d) * (n * d) := by ring _ = ((a / d) * n) * d := by ring

If x solves a·x ≡ b (mod n), then x + k·(n/d) solves it too: adding multiples of n/d preserves solutions (CLRS §31.4).

theorem linear_congruence_shift {a b x n k : ℕ} [NeZero n] (h : a * x ≡ b [MOD n]) : a * (x + n / Nat.gcd a n * k) ≡ b [MOD n] := by have hsplit : a * (x + n / Nat.gcd a n * k) = a * x + a * (n / Nat.gcd a n * k) := by ring rw [hsplit] have h0 : a * (n / Nat.gcd a n * k) ≡ 0 [MOD n] := by rw [← mul_assoc] rw [mul_nat_div_eq a n] apply Nat.modEq_zero_iff_dvd.mpr use (a / Nat.gcd a n) * k ring simpa using (Nat.ModEq.add h h0)

If x₀ and x both solve a·x ≡ b (mod n), then x ≡ x₀ (mod n/d): every solution differs from any given solution by a multiple of n/d (CLRS §31.4). Together with linear_congruence_shift, the d = gcd(a, n) distinct solutions modulo n are x₀, x₀ + n/d, …, x₀ + (d−1)·(n/d).

theorem linear_congruence_all_solutions {a b x x₀ n : ℕ} [NeZero n] (h : a * x ≡ b [MOD n]) (h₀ : a * x₀ ≡ b [MOD n]) : x ≡ x₀ [MOD (n / Nat.gcd a n)] := by let d := Nat.gcd a n have hd : 0 < d := Nat.gcd_pos_of_pos_right a (NeZero.pos n) have ha : a = d * (a / d) := by dsimp [d] simpa [Nat.mul_comm] using (Nat.div_mul_cancel (Nat.gcd_dvd_left a n)).symm have hn : n = d * (n / d) := by dsimp [d] simpa [Nat.mul_comm] using (Nat.div_mul_cancel (Nat.gcd_dvd_right a n)).symm have hax : a * x ≡ a * x₀ [MOD n] := h.trans h₀.symm have hax2 : (a / d) * x ≡ (a / d) * x₀ [MOD (n / d)] := by rw [Nat.ModEq] at hax ⊢ have hsplit : ∀ z : ℕ, (a * z) % n = d * (((a / d) * z) % (n / d)) := by intro z calc (a * z) % n = (d * ((a / d) * z)) % (d * (n / d)) := by have hz : a * z = d * ((a / d) * z) := by conv_lhs => rw [ha] rw [← mul_assoc] rw [hz] conv_lhs => rw [hn] _ = d * (((a / d) * z) % (n / d)) := by rw [Nat.mul_mod_mul_left] have hres : d * (((a / d) * x) % (n / d)) = d * (((a / d) * x₀) % (n / d)) := by rw [← hsplit x, ← hsplit x₀] exact hax exact Nat.eq_of_mul_eq_mul_left hd hres have hcop' : Nat.Coprime (a / d) (n / d) := by dsimp [d] exact Nat.coprime_div_gcd_div_gcd hd exact Nat.ModEq.cancel_left_of_coprime (by simpa [Nat.gcd_comm] using hcop'.gcd_eq_one) hax2

The solutions of a linear congruence (CLRS Theorem 31.10). If x₀ solves a·x ≡ b (mod n), then a value x solves the congruence exactly when x ≡ x₀ (mod n/d) for d = gcd(a, n): the solutions form one residue class modulo n/d.

theorem linear_congruence_solutions {a b x₀ n : ℕ} [NeZero n] (h : a * x₀ ≡ b [MOD n]) : ∀ x : ℕ, (a * x ≡ b [MOD n]) ↔ x ≡ x₀ [MOD (n / Nat.gcd a n)] := by intro x constructor · intro hx exact linear_congruence_all_solutions hx h · intro hx have h1 : a * x ≡ a * x₀ [MOD a * (n / Nat.gcd a n)] := by rw [Nat.ModEq] at hx ⊢ rw [Nat.mul_mod_mul_left, Nat.mul_mod_mul_left] rw [hx] have h2 : a * x ≡ a * x₀ [MOD n] := by rw [mul_nat_div_eq a n] at h1 exact Nat.ModEq.of_dvd (dvd_mul_left n (a / Nat.gcd a n)) h1 exact h2.trans h

The d = gcd(a, n) solutions are distinct modulo n (CLRS Theorem 31.10). The values k·(n/d) for 0 ≤ k < d are pairwise incongruent modulo n, so together with linear_congruence_shift they give exactly d distinct solutions.

theorem linear_congruence_distinct {a n : ℕ} [NeZero n] (k₁ k₂ : ℕ) (Variable name `hk₁` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hk₁ : k₁ < Nat.gcd a n) (hk₂ : k₂ < Nat.gcd a n) (hk : k₁ < k₂) : ¬ k₁ * (n / Nat.gcd a n) ≡ k₂ * (n / Nat.gcd a n) [MOD n] := by intro hc have hd : n ∣ (k₂ - k₁) * (n / Nat.gcd a n) := by have hle : k₁ * (n / Nat.gcd a n) ≤ k₂ * (n / Nat.gcd a n) := by exact Nat.mul_le_mul_right _ (Nat.le_of_lt hk) have hmod : (k₂ * (n / Nat.gcd a n) - k₁ * (n / Nat.gcd a n)) % n = 0 := by have hsub := Nat.ModEq.sub (Nat.le_refl (k₁ * (n / Nat.gcd a n))) hle hc (Nat.ModEq.refl (k₁ * (n / Nat.gcd a n))) have h0 : k₁ * (n / Nat.gcd a n) - k₁ * (n / Nat.gcd a n) = 0 := by omega rw [h0] at hsub simpa [Nat.ModEq, Nat.zero_mod] using hsub.symm have hsub' : (k₂ - k₁) * (n / Nat.gcd a n) = k₂ * (n / Nat.gcd a n) - k₁ * (n / Nat.gcd a n) := by rw [Nat.sub_mul] rw [hsub'] exact Nat.dvd_of_mod_eq_zero hmod have hn0 : 0 < n / Nat.gcd a n := by exact Nat.div_pos (Nat.le_of_dvd (Nat.pos_of_neZero (n := n)) (Nat.gcd_dvd_right a n)) (Nat.gcd_pos_of_pos_right a (Nat.pos_of_neZero (n := n))) have hpos : 0 < (k₂ - k₁) * (n / Nat.gcd a n) := by have hk0 : 0 < k₂ - k₁ := by omega exact Nat.mul_pos hk0 hn0 have hlt : (k₂ - k₁) * (n / Nat.gcd a n) < n := by have hk2 : k₂ - k₁ < Nat.gcd a n := by omega have hdn : (Nat.gcd a n) * (n / Nat.gcd a n) = n := Nat.mul_div_cancel' (Nat.gcd_dvd_right a n) calc (k₂ - k₁) * (n / Nat.gcd a n) < Nat.gcd a n * (n / Nat.gcd a n) := Nat.mul_lt_mul_of_pos_right hk2 hn0 _ = n := hdn have hm0 : (k₂ - k₁) * (n / Nat.gcd a n) = 0 := by have hmod : (k₂ - k₁) * (n / Nat.gcd a n) % n = 0 := Nat.mod_eq_zero_of_dvd hd have hmod' : (k₂ - k₁) * (n / Nat.gcd a n) % n = (k₂ - k₁) * (n / Nat.gcd a n) := Nat.mod_eq_of_lt hlt omega omega

The canonical solution x₀ of a·x ≡ b (mod n): compute (b/d)·(a/d)⁻¹ mod (n/d) with d = gcd(a, n) in ZMod (n/d) and take the representative.

def modularLinearEquationSolution (a b n : ℕ) : ℕ := ((b / Nat.gcd a n : ZMod (n / Nat.gcd a n)) * ((a / Nat.gcd a n : ZMod (n / Nat.gcd a n))⁻¹)).val

In ZMod m, if a' is coprime to m then a'·((b'·a'⁻¹).val) = b': the inverse exists and the reduced representative witnesses the product.

private lemma zmod_solution_of_coprime {m a' b' : ℕ} [NeZero m] (hcop : Nat.Coprime a' m) : (a' : ZMod m) * (((b' : ZMod m) * ((a' : ZMod m)⁻¹)).val : ZMod m) = (b' : ZMod m) := by have hu : IsUnit (a' : ZMod m) := (ZMod.isUnit_iff_coprime a' m).2 hcop calc (a' : ZMod m) * (((b' : ZMod m) * ((a' : ZMod m)⁻¹)).val : ZMod m) = (a' : ZMod m) * ((b' : ZMod m) * ((a' : ZMod m)⁻¹)) := by rw [ZMod.natCast_zmod_val] _ = (b' : ZMod m) * ((a' : ZMod m) * ((a' : ZMod m)⁻¹)) := by ring _ = (b' : ZMod m) * 1 := by rw [ZMod.mul_inv_of_unit (a' : ZMod m) hu] _ = (b' : ZMod m) := by simp

The canonical solution is a solution. For n > 0 and gcd(a, n) ∣ b, modularLinearEquationSolution a b n satisfies a·x₀ ≡ b (mod n).

theorem modularLinearEquationSolution_spec (a b n : ℕ) (hn : 0 < n) (hdvd : Nat.gcd a n ∣ b) : a * modularLinearEquationSolution a b n ≡ b [MOD n] := by let d := Nat.gcd a n let m := n / d have hdpos : 0 < d := Nat.gcd_pos_of_pos_right a hn have hmpos : 0 < m := Nat.div_pos (Nat.le_of_dvd hn (Nat.gcd_dvd_right a n)) hdpos haveI : NeZero m := ⟨Nat.ne_of_gt hmpos⟩ have hcop : Nat.Coprime (a / d) m := by dsimp [m, d] exact Nat.coprime_div_gcd_div_gcd hdpos have hz : (a / d : ZMod m) * (((b / d : ZMod m) * ((a / d : ZMod m)⁻¹)).val : ZMod m) = (b / d : ZMod m) := zmod_solution_of_coprime hcop have hz' : (a / d : ZMod m) * (modularLinearEquationSolution a b n : ZMod m) = (b / d : ZMod m) := by simpa [modularLinearEquationSolution, d, m] using hz have hmod : a / d * modularLinearEquationSolution a b n ≡ b / d [MOD m] := by rw [← ZMod.natCast_eq_natCast_iff] rw [Nat.cast_mul] exact hz' rw [Nat.ModEq] at hmod have ha' : a = d * (a / d) := by dsimp [d] exact (Nat.mul_div_cancel' (Nat.gcd_dvd_left a n)).symm have hb' : b = d * (b / d) := by dsimp [d] exact (Nat.mul_div_cancel' hdvd).symm have hn' : n = d * m := by dsimp [m, d] exact (Nat.mul_div_cancel' (Nat.gcd_dvd_right a n)).symm have hscaled : d * (a / d * modularLinearEquationSolution a b n) % (d * m) = d * (b / d) % (d * m) := by rw [Nat.mul_mod_mul_left, Nat.mul_mod_mul_left] rw [hmod] have hax : a * modularLinearEquationSolution a b n = d * (a / d * modularLinearEquationSolution a b n) := by calc a * modularLinearEquationSolution a b n = (d * (a / d)) * modularLinearEquationSolution a b n := congrArg (fun z => z * modularLinearEquationSolution a b n) ha' _ = d * (a / d * modularLinearEquationSolution a b n) := by ring have h1 : a * modularLinearEquationSolution a b n % n = d * (a / d * modularLinearEquationSolution a b n) % (d * m) := congrArg₂ Nat.mod hax hn' have h2 : d * (b / d) % (d * m) = b % n := congrArg₂ Nat.mod hb'.symm hn'.symm calc a * modularLinearEquationSolution a b n % n = d * (a / d * modularLinearEquationSolution a b n) % (d * m) := h1 _ = d * (b / d) % (d * m) := hscaled _ = b % n := h2

For n > 0, n / gcd(a, n) is positive (the reduced modulus of §31.4).

private lemma gcd_div_pos (a n : ℕ) [NeZero n] : 0 < n / Nat.gcd a n := by have hn : 0 < n := NeZero.pos n have hdn : Nat.gcd a n ∣ n := Nat.gcd_dvd_right a n have hdpos : 0 < Nat.gcd a n := Nat.gcd_pos_of_pos_right a hn exact Nat.div_pos (Nat.le_of_dvd hn hdn) hdpos

MODULAR-LINEAR-EQUATION-SOLVER (CLRS §31.4). Enumerate the d = gcd(a, n) distinct solutions of a·x ≡ b (mod n) as x₀, x₀ + n/d, …, x₀ + (d−1)·(n/d); return the empty list when no solution exists (gcd(a, n) ∤ b).

def modularLinearEquationSolver (a b n : ℕ) : List ℕ := if Variable name `h` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`h : Nat.gcd a n ∣ b then if Variable name `hm` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hm : 0 < n / Nat.gcd a n then (List.range (Nat.gcd a n)).map (fun k => modularLinearEquationSolution a b n + k * (n / Nat.gcd a n)) else [] else []

The solver returns exactly d solutions when solvable, none otherwise (CLRS §31.4).

theorem modularLinearEquationSolver_length (a b n : ℕ) [NeZero n] : (modularLinearEquationSolver a b n).length = if Nat.gcd a n ∣ b then Nat.gcd a n else 0 := by by_cases h : Nat.gcd a n ∣ b · have hm : 0 < n / Nat.gcd a n := gcd_div_pos a n simp [modularLinearEquationSolver, h, hm] · simp [modularLinearEquationSolver, h]

Every enumerated value is a solution and lies below n (CLRS §31.4).

theorem modularLinearEquationSolver_sound (a b n : ℕ) [NeZero n] : ∀ x ∈ modularLinearEquationSolver a b n, x < n ∧ a * x ≡ b [MOD n] := by have hn : 0 < n := NeZero.pos n intro x hx by_cases hdvd : Nat.gcd a n ∣ b · have hm : 0 < n / Nat.gcd a n := gcd_div_pos a n simp [modularLinearEquationSolver, hdvd, hm] at hx rcases hx with ⟨k, hk, hxeq⟩ rw [← hxeq] haveI : NeZero (n / Nat.gcd a n) := NeZero.of_pos hm have hx₀lt : modularLinearEquationSolution a b n < n / Nat.gcd a n := by dsimp [modularLinearEquationSolution] exact ZMod.val_lt (((b / Nat.gcd a n : ZMod (n / Nat.gcd a n)) * ((a / Nat.gcd a n : ZMod (n / Nat.gcd a n))⁻¹))) have hklt : k < Nat.gcd a n := hk constructor · have hn' : n = Nat.gcd a n * (n / Nat.gcd a n) := by exact (Nat.mul_div_cancel' (Nat.gcd_dvd_right a n)).symm calc modularLinearEquationSolution a b n + k * (n / Nat.gcd a n) < n / Nat.gcd a n + k * (n / Nat.gcd a n) := Nat.add_lt_add_right hx₀lt (k * (n / Nat.gcd a n)) _ = (k + 1) * (n / Nat.gcd a n) := by ring _ ≤ Nat.gcd a n * (n / Nat.gcd a n) := Nat.mul_le_mul_right (n / Nat.gcd a n) (by omega) _ = n := hn'.symm · have hsol := modularLinearEquationSolution_spec a b n hn hdvd have hshift := linear_congruence_shift (n := n) (a := a) (b := b) (x := modularLinearEquationSolution a b n) (k := k) hsol simpa [Nat.mul_comm] using hshift · simp [modularLinearEquationSolver, hdvd] at hx

Every solution x < n is enumerated (CLRS §31.4).

theorem modularLinearEquationSolver_complete (a b n : ℕ) [NeZero n] {x : ℕ} (hx : a * x ≡ b [MOD n]) (hlt : x < n) : x ∈ modularLinearEquationSolver a b n := by have hn : 0 < n := NeZero.pos n have hmpos : 0 < n / Nat.gcd a n := gcd_div_pos a n have hdvd : Nat.gcd a n ∣ b := by have hda : Nat.gcd a n ∣ a := Nat.gcd_dvd_left a n have hdn : Nat.gcd a n ∣ n := Nat.gcd_dvd_right a n have hx_d : a * x ≡ b [MOD Nat.gcd a n] := Nat.ModEq.of_dvd hdn hx have ha_d : a ≡ 0 [MOD Nat.gcd a n] := (Nat.modEq_zero_iff_dvd).2 hda have hb0 : b ≡ 0 [MOD Nat.gcd a n] := by have h0 : a * x ≡ 0 * x [MOD Nat.gcd a n] := Nat.ModEq.mul ha_d (Nat.ModEq.refl x) simpa using (hx_d.symm.trans h0) exact (Nat.modEq_zero_iff_dvd).1 hb0 have hx₀sol : a * modularLinearEquationSolution a b n ≡ b [MOD n] := modularLinearEquationSolution_spec a b n hn hdvd have hcong : x ≡ modularLinearEquationSolution a b n [MOD n / Nat.gcd a n] := by have h := (linear_congruence_solutions (a := a) (b := b) (x₀ := modularLinearEquationSolution a b n) (n := n) hx₀sol) x exact h.1 hx haveI : NeZero (n / Nat.gcd a n) := NeZero.of_pos hmpos have hx₀lt : modularLinearEquationSolution a b n < n / Nat.gcd a n := by dsimp [modularLinearEquationSolution] exact ZMod.val_lt (((b / Nat.gcd a n : ZMod (n / Nat.gcd a n)) * ((a / Nat.gcd a n : ZMod (n / Nat.gcd a n))⁻¹))) have hmod : x % (n / Nat.gcd a n) = modularLinearEquationSolution a b n := by rw [Nat.ModEq] at hcong rw [hcong] rw [Nat.mod_eq_of_lt hx₀lt] have hxeq : x = modularLinearEquationSolution a b n + (x / (n / Nat.gcd a n)) * (n / Nat.gcd a n) := by have hdiv := Nat.div_add_mod x (n / Nat.gcd a n) nth_rewrite 1 [← hdiv] rw [hmod] ring have hklt : x / (n / Nat.gcd a n) < Nat.gcd a n := by have hxlt : x < (n / Nat.gcd a n) * Nat.gcd a n := by have hn' : n = (n / Nat.gcd a n) * Nat.gcd a n := by simpa [Nat.mul_comm] using (Nat.mul_div_cancel' (Nat.gcd_dvd_right a n)).symm rw [hn'] at hlt exact hlt exact Nat.div_lt_of_lt_mul hxlt rw [modularLinearEquationSolver] rw [dif_pos hdvd] rw [dif_pos hmpos] rw [List.mem_map] refine ⟨x / (n / Nat.gcd a n), ⟨(List.mem_range).2 hklt, hxeq.symm⟩⟩

The solver produces pairwise-distinct solutions (CLRS §31.4).

theorem modularLinearEquationSolver_nodup (a b n : ℕ) [NeZero n] : (modularLinearEquationSolver a b n).Nodup := by by_cases hdvd : Nat.gcd a n ∣ b · have hm : 0 < n / Nat.gcd a n := gcd_div_pos a n simp [modularLinearEquationSolver, hdvd, hm] refine List.Nodup.map ?hf List.nodup_range intro k₁ k₂ hk have hkm : k₁ * (n / Nat.gcd a n) = k₂ * (n / Nat.gcd a n) := Nat.add_left_cancel hk exact Nat.mul_right_cancel hm hkm · simp [modularLinearEquationSolver, hdvd]
end Chapter31end CLRS
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_3_Modular_Arithmetic

31.5. The Chinese Remainder Theorem

CLRS §31.5: the Chinese remainder theorem (Theorem 31.27) — if the moduli n₁, …, nₖ are pairwise relatively prime, then the system of congruences x ≡ aᵢ (mod nᵢ) has a unique solution modulo n₁·…·nₖ.

This section formalizes both the two-modulus form and the general list form of the theorem. Mathlib's Nat.chineseRemainder supplies existence for two moduli; Nat.chineseRemainderOfList handles the general case.

Main results:

  • Theorem chinese_remainder_two: for coprime n m, the system x ≡ a (mod n), x ≡ b (mod m) has a solution.

  • Theorem chinese_remainder_unique: any two solutions agree modulo n·m.

  • Theorem chinese_remainder_general (Theorem 31.27, general form): for a list of pairwise-coprime moduli, a solution exists, unique modulo the product.

  • Ring equivalence zmod_chineseRemainder: the ring-isomorphism packaging — for coprime m n, ZMod (m·n) ≃+* ZMod m × ZMod n, with the projection-recovery lemmas zmod_chineseRemainder_fst and zmod_chineseRemainder_snd.

Notation:

  • a ≡ b [MOD n] : Nat.ModEq.

  • Nat.Coprime n m : gcd n m = 1.

Deferred: none.

namespace CLRSnamespace Chapter31

Chinese remainder theorem (CLRS Theorem 31.27, two moduli). If n and m are coprime, the system of congruences x ≡ a (mod n) and x ≡ b (mod m) has a solution.

theorem chinese_remainder_two {n m a b : ℕ} (hcop : Nat.Coprime n m) : ∃ x : ℕ, x ≡ a [MOD n] ∧ x ≡ b [MOD m] := by exact ⟨(Nat.chineseRemainder hcop a b).1, (Nat.chineseRemainder hcop a b).2⟩

Chinese remainder theorem, uniqueness. For coprime n and m, any two solutions of the same system x ≡ a (mod n), x ≡ b (mod m) agree modulo n·m.

theorem chinese_remainder_unique {n m a b x y : ℕ} (hcop : Nat.Coprime n m) (hx : x ≡ a [MOD n] ∧ x ≡ b [MOD m]) (hy : y ≡ a [MOD n] ∧ y ≡ b [MOD m]) : x ≡ y [MOD n * m] := by have hxcr := Nat.chineseRemainder_modEq_unique hcop hx.1 hx.2 have hycr := Nat.chineseRemainder_modEq_unique hcop hy.1 hy.2 exact hxcr.trans hycr.symm

Chinese remainder theorem (two moduli, bundled). For coprime n m, the system has a solution, unique modulo n·m.

theorem chinese_remainder {n m a b : ℕ} (hcop : Nat.Coprime n m) : ∃ x : ℕ, x ≡ a [MOD n] ∧ x ≡ b [MOD m] ∧ ∀ y : ℕ, y ≡ a [MOD n] → y ≡ b [MOD m] → x ≡ y [MOD n * m] := by obtain ⟨x, hx₁, hx₂⟩ := chinese_remainder_two hcop refine ⟨x, hx₁, hx₂, ?_⟩ intro y hy₁ hy₂ exact chinese_remainder_unique hcop ⟨hx₁, hx₂⟩ ⟨hy₁, hy₂⟩

Chinese remainder theorem, general form (CLRS Theorem 31.27). For a list of pairwise-coprime moduli s i and residues a i, the system of congruences x ≡ a i (mod s i) has a solution, unique modulo the product of the moduli.

theorem chinese_remainder_general {ι : Type} (a s : ι → ℕ) (l : List ι) (co : List.Pairwise (Function.onFun Nat.Coprime s) l) : ∃ x : ℕ, (∀ i ∈ l, x ≡ a i [MOD s i]) ∧ ∀ y : ℕ, (∀ i ∈ l, y ≡ a i [MOD s i]) → x ≡ y [MOD (List.map s l).prod] := by let crt := Nat.chineseRemainderOfList a s l co refine ⟨crt.1, ?_⟩ constructor · exact crt.2 · intro y hy exact (Nat.chineseRemainderOfList_modEq_unique a s l co hy).symm

Chinese remainder theorem, ring-isomorphism form. For coprime m n, the rings ZMod (m·n) and ZMod m × ZMod n are isomorphic: a residue modulo m·n is mapped to the pair of its residues modulo m and modulo n.

def zmod_chineseRemainder {m n : ℕ} (h : Nat.Coprime m n) : ZMod (m * n) ≃+* ZMod m × ZMod n := ZMod.chineseRemainder h

The first projection of the CRT ring isomorphism recovers the residue modulo m: (zmod_chineseRemainder h x).1 = x mod m.

theorem zmod_chineseRemainder_fst {m n : ℕ} (h : Nat.Coprime m n) (x : ZMod (m * n)) : (zmod_chineseRemainder h x).1 = ZMod.cast x := by simp [zmod_chineseRemainder, ZMod.chineseRemainder, ZMod.castHom_apply]

The second projection of the CRT ring isomorphism recovers the residue modulo n: (zmod_chineseRemainder h x).2 = x mod n.

theorem zmod_chineseRemainder_snd {m n : ℕ} (h : Nat.Coprime m n) (x : ZMod (m * n)) : (zmod_chineseRemainder h x).2 = ZMod.cast x := by simp [zmod_chineseRemainder, ZMod.chineseRemainder, ZMod.castHom_apply]
end Chapter31end CLRS
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_3_Modular_Arithmetic

31.6. Powers of an Element

CLRS §31.6: modular exponentiation and the two theorems that make powers modulo a number tractable — Fermat's little theorem (prime modulus) and Euler's theorem (general modulus, via the totient function).

Main results:

  • modularExponentiation + modularExponentiation_spec: computing a^b mod n by repeated squaring returns a number congruent to a^b modulo n.

  • Theorem fermat_little_theorem (CLRS Theorem 31.30): for prime p, a^p ≡ a (mod p).

  • Theorem euler_theorem (CLRS Euler's theorem): for gcd(a, n) = 1, a^φ(n) ≡ 1 (mod n).

Notation:

  • a ≡ b [MOD n] : Nat.ModEq.

  • Nat.totient n : Euler's totient φ(n).

Deferred: none (the repeated-squaring recursion with its explicit operation count is proved by modExpWithCount and modExpWithCount_count_le).

namespace CLRSnamespace Chapter31

MODULAR-EXPONENTIATION (CLRS §31.6). Compute a^b mod n. The value is the remainder of a^b upon division by n; the repeated-squaring algorithm computes it in O(log b) multiplications.

def modularExponentiation (a b n : ℕ) : ℕ := a ^ b % n

MODULAR-EXPONENTIATION is correct: modularExponentiation a b n is congruent to a^b modulo n.

theorem modularExponentiation_spec (a b n : ℕ) : modularExponentiation a b n ≡ a ^ b [MOD n] := by rw [Nat.ModEq] unfold modularExponentiation rw [Nat.mod_mod]

Fermat's little theorem (CLRS Theorem 31.30). For a prime p and any a, a^p ≡ a (mod p).

try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false` theorem fermat_little_theorem {p a : ℕ} (hp : Nat.Prime p) : a ^ p ≡ a [MOD p] := by letI : Fact (Nat.Prime p) := ⟨hp⟩ have h := ZMod.pow_card (p := p) (x := (a : ZMod p)) exact (ZMod.natCast_eq_natCast_iff (a ^ p) a p).mp (by try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [Nat.cast_pow] using h)

Euler's theorem (CLRS §31.6). If gcd(a, n) = 1 and 1 < n, then a^φ(n) ≡ 1 (mod n).

theorem euler_theorem {a n : ℕ} (hn : 1 < n) (hcop : Nat.Coprime a n) : a ^ Nat.totient n ≡ 1 [MOD n] := by rw [Nat.ModEq] rw [Nat.pow_totient_mod_eq_one hn hcop] rw [Nat.mod_eq_of_lt hn]

MODULAR-EXPONENTIATION by repeated squaring (CLRS §31.6). Compute a^b mod n with the square-and-multiply recursion that reads the binary digits of b least-significant first: each level squares the running residue, and a 1 bit additionally multiplies by a. The returned pair carries the residue together with the number of modular multiplications performed (one square per bit plus one multiply per 1 bit).

def modExpWithCount (a n : ℕ) (b : ℕ) : ℕ × ℕ := if hb : b = 0 then (1 % n, 0) else let r := modExpWithCount a n (b / 2) let d := (r.1 * r.1) % n if b % 2 = 0 then (d, r.2 + 1) else ((d * a) % n, r.2 + 2) termination_by b decreasing_by exact Nat.div_lt_self (Nat.pos_of_ne_zero hb) (by decide)

One square step: (a^(b/2) mod n)² mod n = a^(2·(b/2)) mod n.

private lemma mod_sq_mod (a n b : ℕ) : ((a ^ (b / 2) % n) * (a ^ (b / 2) % n)) % n = a ^ (2 * (b / 2)) % n := by rw [← Nat.mul_mod] rw [← Nat.pow_add] rw [show b / 2 + b / 2 = 2 * (b / 2) by omega]

(X mod n) · a mod n = X · a mod n: the residue of a product depends only on the residue of X.

private lemma mod_mul_mod_left (X a n : ℕ) : ((X % n) * a) % n = (X * a) % n := by rw [Nat.mul_mod (X % n) a n] rw [Nat.mod_mod] rw [← Nat.mul_mod X a n]

One square-then-multiply step: ((a^(b/2) mod n)² mod n) · a mod n = a^(2·(b/2)+1) mod n.

private lemma mod_sq_mul_mod (a n b : ℕ) : (((a ^ (b / 2) % n) * (a ^ (b / 2) % n)) % n * a) % n = a ^ (2 * (b / 2) + 1) % n := by rw [mod_sq_mod] rw [mod_mul_mod_left] rw [← Nat.pow_succ]

MODULAR-EXPONENTIATION is correct: the repeated-squaring recursion returns exactly a^b mod n (CLRS §31.6).

theorem modExpWithCount_spec (a n b : ℕ) : (modExpWithCount a n b).1 = a ^ b % n := by induction b using Nat.strong_induction_on with | h b ih => by_cases hb : b = 0 · simp [modExpWithCount, hb] · have hlt : b / 2 < b := Nat.div_lt_self (Nat.pos_of_ne_zero hb) (by decide) have ihr := ih (b / 2) hlt unfold modExpWithCount rw [dif_neg hb] simp only split_ifs with heven · -- even: square only calc (((modExpWithCount a n (b / 2)).1 * (modExpWithCount a n (b / 2)).1) % n) = ((a ^ (b / 2) % n) * (a ^ (b / 2) % n)) % n := by rw [ihr] _ = a ^ (2 * (b / 2)) % n := by rw [mod_sq_mod] _ = a ^ b % n := by rw [show 2 * (b / 2) = b by have hdiv := Nat.div_add_mod b 2 omega] · -- odd: square and multiply by `a` have hodd : b % 2 = 1 := by have hmodlt := Nat.mod_lt b (by decide : 0 < 2) omega calc ((((modExpWithCount a n (b / 2)).1 * (modExpWithCount a n (b / 2)).1) % n) * a) % n = (((a ^ (b / 2) % n) * (a ^ (b / 2) % n)) % n * a) % n := by rw [ihr] _ = a ^ (2 * (b / 2) + 1) % n := by rw [mod_sq_mul_mod] _ = a ^ b % n := by rw [show 2 * (b / 2) + 1 = b by have hdiv := Nat.div_add_mod b 2 omega]

Nat.size (b / 2) + 1 = Nat.size b for nonzero b (dropping the least significant bit reduces the bit length by exactly one).

private lemma size_div_two_succ {b : ℕ} (hb : b ≠ 0) : Nat.size (b / 2) + 1 = Nat.size b := by have hdecomp : Nat.bit b.bodd b.div2 = b := Nat.bit_bodd_div2 b have hnb : Nat.bit b.bodd b.div2 ≠ 0 := by rw [hdecomp] exact hb calc Nat.size (b / 2) + 1 = b.div2.size + 1 := by rw [Nat.div2_val] _ = b.div2.size.succ := rfl _ = Nat.size (Nat.bit b.bodd b.div2) := by rw [Nat.size_bit hnb] _ = Nat.size b := by rw [hdecomp]

MODULAR-EXPONENTIATION uses O(log b) multiplications: the repeated squaring recursion performs at most 2 · Nat.size b modular multiplications, where Nat.size b is the number of bits of b (CLRS §31.6).

theorem modExpWithCount_count_le (a n b : ℕ) : (modExpWithCount a n b).2 ≤ 2 * Nat.size b := by induction b using Nat.strong_induction_on with | h b ih => by_cases hb : b = 0 · simp [modExpWithCount, hb] · have hlt : b / 2 < b := Nat.div_lt_self (Nat.pos_of_ne_zero hb) (by decide) have ihr := ih (b / 2) hlt unfold modExpWithCount rw [dif_neg hb] simp only split_ifs with heven · have hs : Nat.size (b / 2) + 1 = Nat.size b := size_div_two_succ hb omega · have hs : Nat.size (b / 2) + 1 = Nat.size b := size_div_two_succ hb omega
end Chapter31end CLRS
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_6_Powers_Of_An_Element import CLRSLean.FourthEdition.Chapter_31.Section_31_5_Chinese_Remainder_Theorem

31.7. The RSA Public-Key Cryptosystem

CLRS §31.7: the RSA public-key cryptosystem. With n = p·q for distinct primes p, q, φ(n) = (p−1)(q−1), a public exponent e coprime to φ(n), and the private exponent d ≡ e⁻¹ (mod φ(n)), encryption c = m^e mod n and decryption m = c^d mod n are mutually inverse.

Main results:

  • Theorem totient_mul_prime: for distinct primes p q, φ(p·q) = (p−1)·(q−1).

  • Theorem rsa_correct (CLRS Theorem 31.36): if e·d ≡ 1 (mod φ(n)) and gcd(m, n) = 1, then m^(e·d) ≡ m (mod n) — decryption undoes encryption.

  • Theorem rsa_correct_general (CLRS Theorem 31.36, general message): for distinct primes p q and e·d ≡ 1 (mod (p−1)(q−1)), m^(e·d) ≡ m (mod p·q) for every m — via Fermat modulo each prime and the Chinese remainder theorem.

Notation:

  • Nat.totient n : Euler's totient φ(n).

  • a ≡ b [MOD n] : Nat.ModEq.

Deferred: the RSA security (one-way function) claims stay out of scope; the key-generation (rsaKeyGen, rsaPrivateExponent) and the repeated-squaring running-time bounds (rsaEncrypt_count_le, rsaDecrypt_count_le) are proved.

namespace CLRSnamespace Chapter31

For distinct primes p and q, φ(p·q) = (p−1)·(q−1): the totient of the RSA modulus.

theorem totient_mul_prime (p q : ℕ) (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) : (p * q).totient = (p - 1) * (q - 1) := by have hcop : Nat.Coprime p q := by rw [Nat.Prime.coprime_iff_not_dvd hp] intro hpq_dvd rcases (Nat.Prime.eq_one_or_self_of_dvd hq p hpq_dvd) with h1 | heq · exfalso have hp2 : 2 ≤ p := Nat.Prime.two_le hp omega · exact hpq heq rw [Nat.totient_mul hcop] rw [Nat.totient_prime hp, Nat.totient_prime hq]

RSA is correct (CLRS Theorem 31.36). If the exponents satisfy e·d ≡ 1 (mod φ(n)) and m is coprime to the modulus n, then m^(e·d) ≡ m (mod n): raising to the power e·d (encryption followed by decryption, or vice versa) recovers m.

theorem rsa_correct {m e d n : ℕ} (hle : 1 ≤ e * d) (hmed : e * d ≡ 1 [MOD Nat.totient n]) (hcop : Nat.Coprime m n) : m ^ (e * d) ≡ m [MOD n] := by rcases (Nat.modEq_iff_exists_eq_add hle).mp hmed.symm with ⟨k, hk⟩ have hE : m ^ Nat.totient n ≡ 1 [MOD n] := Nat.ModEq.pow_totient hcop have h1 : m ^ (e * d) ≡ m ^ (1 + Nat.totient n * k) [MOD n] := by rw [hk] have h2 : m ^ (1 + Nat.totient n * k) = m * (m ^ Nat.totient n) ^ k := by rw [pow_add, pow_mul, pow_one] have h3 : m * (m ^ Nat.totient n) ^ k ≡ m * 1 ^ k [MOD n] := by exact Nat.ModEq.mul (Nat.ModEq.refl m) (Nat.ModEq.pow k hE) have hmid : m ^ (1 + Nat.totient n * k) ≡ m * 1 ^ k [MOD n] := by rw [h2] exact h3 exact (by simpa using (h1.trans hmid))

For prime p and e·d ≡ 1 (mod p−1), m^(e·d) ≡ m (mod p): the RSA exponentiation recovers m modulo each prime factor of the modulus (the p-side of the general RSA correctness).

lemma rsa_pow_cong {p m e d : ℕ} (hp : Nat.Prime p) (hle : 1 ≤ e * d) (hmed : e * d ≡ 1 [MOD p - 1]) : m ^ (e * d) ≡ m [MOD p] := by letI : Fact (Nat.Prime p) := ⟨hp⟩ haveI : NeZero p := ⟨Nat.Prime.ne_zero hp⟩ rcases (Nat.modEq_iff_exists_eq_add hle).mp hmed.symm with ⟨k, hk⟩ have hz0 : (m : ZMod p) ^ (e * d) = (m : ZMod p) := by rw [hk, pow_add, pow_one, pow_mul] by_cases hm0 : (m : ZMod p) = 0 · rw [hm0] simp · have hp1 : (m : ZMod p) ^ (p - 1) = 1 := by simpa [hm0] using (ZMod.pow_card_sub_one (p := p) (a := (m : ZMod p))) rw [hp1] simp rw [Nat.ModEq] calc (m ^ (e * d)) % p = (↑(m ^ (e * d)) : ZMod p).val := (ZMod.val_natCast p (m ^ (e * d))).symm _ = ((m : ZMod p) ^ (e * d)).val := by rw [Nat.cast_pow] _ = (m : ZMod p).val := by rw [hz0] _ = m % p := ZMod.val_natCast p m

Distinct primes are coprime.

lemma prime_coprime {p q : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) : Nat.Coprime p q := by rw [Nat.Prime.coprime_iff_not_dvd hp] intro hpq_dvd rcases (Nat.Prime.eq_one_or_self_of_dvd hq p hpq_dvd) with h1 | heq · exfalso have hp2 : 2 ≤ p := Nat.Prime.two_le hp omega · exact hpq heq

RSA is correct for every message (CLRS Theorem 31.36). For distinct primes p q, n = p·q, and exponents with e·d ≡ 1 (mod (p−1)(q−1)), m^(e·d) ≡ m (mod p·q) for every m — including messages sharing a factor with n. The proof shows the congruence modulo each prime factor (rsa_pow_cong, which covers both p | m and p ∤ m via Fermat) and combines them with the Chinese remainder theorem.

try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false` theorem rsa_correct_general {p q m e d : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hle : 1 ≤ e * d) (hmed : e * d ≡ 1 [MOD (p - 1) * (q - 1)]) : m ^ (e * d) ≡ m [MOD p * q] := by have hp_cong : m ^ (e * d) ≡ m [MOD p] := by apply rsa_pow_cong hp hle rw [Nat.ModEq] at hmed ⊢ have hd : p - 1 ∣ (p - 1) * (q - 1) := by try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [Nat.mul_comm] using (dvd_mul_left (p - 1) (q - 1)) rw [← Nat.mod_mod_of_dvd (e * d) hd] rw [hmed] rw [Nat.mod_mod_of_dvd 1 hd] have hq_cong : m ^ (e * d) ≡ m [MOD q] := by apply rsa_pow_cong hq hle rw [Nat.ModEq] at hmed ⊢ have hd : q - 1 ∣ (p - 1) * (q - 1) := by try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [Nat.mul_comm] using (dvd_mul_right (q - 1) (p - 1)) rw [← Nat.mod_mod_of_dvd (e * d) hd] rw [hmed] rw [Nat.mod_mod_of_dvd 1 hd] exact chinese_remainder_unique (prime_coprime hp hq hpq) ⟨hp_cong, hq_cong⟩ ⟨Nat.ModEq.refl m, Nat.ModEq.refl m⟩

RSA private exponent. For RSA modulus n = p·q, the private exponent d is the inverse of the public exponent e modulo φ(n) = (p−1)(q−1), computed in ZMod ((p−1)(q−1)) and read back as a natural number below the modulus.

def rsaPrivateExponent (p q e : ℕ) : ℕ := (((e : ZMod ((p - 1) * (q - 1)))⁻¹).val)

The private exponent is below the totient: d < φ(n).

theorem rsaPrivateExponent_lt (p q e : ℕ) (hφ : 0 < (p - 1) * (q - 1)) : rsaPrivateExponent p q e < (p - 1) * (q - 1) := by dsimp [rsaPrivateExponent] haveI : NeZero ((p - 1) * (q - 1)) := NeZero.of_pos hφ exact ZMod.val_lt (((e : ZMod ((p - 1) * (q - 1)))⁻¹))

RSA key generation is correct. For distinct primes p q, if the public exponent e is coprime to φ(n) = (p−1)(q−1), then the private exponent d = rsaPrivateExponent p q e satisfies e·d ≡ 1 (mod φ(n)).

theorem rsaPrivateExponent_spec (p q e : ℕ) (hp : Nat.Prime p) (hq : Nat.Prime q) (Variable name `hpq` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hpq : p ≠ q) (hcop : Nat.Coprime e ((p - 1) * (q - 1))) : e * rsaPrivateExponent p q e ≡ 1 [MOD (p - 1) * (q - 1)] := by have hφ : 0 < (p - 1) * (q - 1) := by have hp2 : 2 ≤ p := Nat.Prime.two_le hp have hq2 : 2 ≤ q := Nat.Prime.two_le hq exact Nat.mul_pos (by omega : 0 < p - 1) (by omega : 0 < q - 1) letI : NeZero ((p - 1) * (q - 1)) := NeZero.of_pos hφ have hu : IsUnit (e : ZMod ((p - 1) * (q - 1))) := (ZMod.isUnit_iff_coprime e ((p - 1) * (q - 1))).2 hcop have hinv : (e : ZMod ((p - 1) * (q - 1))) * ((e : ZMod ((p - 1) * (q - 1)))⁻¹) = 1 := ZMod.mul_inv_of_unit (e : ZMod ((p - 1) * (q - 1))) hu have hz : (e : ZMod ((p - 1) * (q - 1))) * (rsaPrivateExponent p q e : ZMod ((p - 1) * (q - 1))) = 1 := by simpa [rsaPrivateExponent] using hinv rw [← ZMod.natCast_eq_natCast_iff] rw [Nat.cast_mul] simpa using hz

RSA encryption (CLRS §31.7). Encrypt a message m with public key (e, n) by the repeated-squaring exponentiation of §31.6: m^e mod n.

def rsaEncrypt (e n m : ℕ) : ℕ := (modExpWithCount m n e).1

RSA encryption computes m^e mod n (CLRS §31.7).

theorem rsaEncrypt_spec (e n m : ℕ) : rsaEncrypt e n m = m ^ e % n := modExpWithCount_spec m n e

RSA encryption uses O(log e) modular multiplications (CLRS §31.7).

theorem rsaEncrypt_count_le (e n m : ℕ) : (modExpWithCount m n e).2 ≤ 2 * Nat.size e := modExpWithCount_count_le m n e

RSA decryption (CLRS §31.7). Decrypt a ciphertext c with private key (d, n) by the repeated-squaring exponentiation of §31.6: c^d mod n.

def rsaDecrypt (d n c : ℕ) : ℕ := (modExpWithCount c n d).1

RSA decryption computes c^d mod n (CLRS §31.7).

theorem rsaDecrypt_spec (d n c : ℕ) : rsaDecrypt d n c = c ^ d % n := modExpWithCount_spec c n d

RSA decryption uses O(log d) modular multiplications (CLRS §31.7).

theorem rsaDecrypt_count_le (d n c : ℕ) : (modExpWithCount c n d).2 ≤ 2 * Nat.size d := modExpWithCount_count_le c n d

RSA key generation (CLRS §31.7). From distinct primes p q and a public exponent e, produce the triple (n, e, d) with n = p·q and d = e⁻¹ mod φ(n).

def rsaKeyGen (p q e : ℕ) : ℕ × ℕ × ℕ := (p * q, e, rsaPrivateExponent p q e)

RSA key generation produces a valid key pair: the public exponent e and the generated private exponent d satisfy e·d ≡ 1 (mod φ(n)).

theorem rsaKeyGen_spec (p q e : ℕ) (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hcop : Nat.Coprime e ((p - 1) * (q - 1))) : (rsaKeyGen p q e).2.1 * (rsaKeyGen p q e).2.2 ≡ 1 [MOD (p - 1) * (q - 1)] := by simpa [rsaKeyGen] using rsaPrivateExponent_spec p q e hp hq hpq hcop
end Chapter31end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_31.Section_31_7_RSA.KeyRoundTrip

Generated keys from supplied distinct primes and a coprime public exponent round-trip every message modulo the generated modulus, including messages sharing a factor with it. In-range messages are recovered exactly. This is key assembly and algebraic correctness, not prime generation or a security claim.

namespace CLRS.Chapter31 theorem rsaKeyGen_roundTrip_mod (p q e m : Nat) (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hcop : Nat.Coprime e ((p - 1) * (q - 1))) : rsaDecrypt (rsaKeyGen p q e).2.2 (rsaKeyGen p q e).1 (rsaEncrypt (rsaKeyGen p q e).2.1 (rsaKeyGen p q e).1 m) = m % (rsaKeyGen p q e).1 := by have hp2 := hp.two_le have hq2 := hq.two_le have hφ : 1 < (p - 1) * (q - 1) := by by_cases hp3 : 3 ≤ p · have hp' : 2 ≤ p - 1 := by omega have hq' : 1 ≤ q - 1 := by omega nlinarith · have hp' : p = 2 := by omega have hq' : 2 ≤ q - 1 := by omega rw [hp']; omega have hc := rsaPrivateExponent_spec p q e hp hq hpq hcop have hpos : 1 ≤ e * rsaPrivateExponent p q e := by by_contra hn have hz : e * rsaPrivateExponent p q e = 0 := by omega rw [Nat.ModEq, hz, Nat.zero_mod, Nat.mod_eq_of_lt hφ] at hc omega simp only [rsaKeyGen, rsaDecrypt_spec, rsaEncrypt_spec] rw [Nat.pow_mod, Nat.mod_mod, ← Nat.pow_mod, ← pow_mul] exact rsa_correct_general hp hq hpq hpos hc theorem rsaKeyGen_roundTrip (p q e m : Nat) (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hcop : Nat.Coprime e ((p - 1) * (q - 1))) (hm : m < p * q) : rsaDecrypt (rsaKeyGen p q e).2.2 (rsaKeyGen p q e).1 (rsaEncrypt (rsaKeyGen p q e).2.1 (rsaKeyGen p q e).1 m) = m := by rw [rsaKeyGen_roundTrip_mod p q e m hp hq hpq hcop] exact Nat.mod_eq_of_lt hmend CLRS.Chapter31
Imports
import Mathlib import CLRSLean.FourthEdition.Chapter_31.Section_31_5_Chinese_Remainder_Theorem import CLRSLean.FourthEdition.Chapter_31.Section_31_6_Powers_Of_An_Element

31.8. Primality Testing

CLRS §31.8: probabilistic primality testing. Fermat's little theorem gives a candidate test — for prime p and gcd(a, p) = 1, a^(p−1) ≡ 1 (mod p) — so a value n for which some a fails this congruence cannot be prime. The PSEUDOPRIME test checks 2^(n−1) ≡ 1 (mod n).

Main results:

  • Theorem fermat_test (CLRS Theorem 31.31): for a prime p and a coprime to p, a^(p−1) ≡ 1 (mod p).

  • Definition fermatPseudoprime: a composite n with a^(n−1) ≡ 1 (mod n) for a given a.

  • pseudoprime + pseudoprime_correct (CLRS PSEUDOPRIME): the executable test returns whether 2^(n−1) ≡ 1 (mod n).

  • isCarmichael (Carmichael numbers): a composite n passing the Fermat test for every a coprime to n; a Carmichael number is a Fermat pseudoprime to every coprime base (carmichael_fermatPseudoprime). isCarmichael_561 shows that 561 is a Carmichael number, so PSEUDOPRIME cannot certify primality. The helper modeq_of_coprime_mul combines congruences under coprime moduli.

  • Miller-Rabin: strongTestParams writes n−1 = 2^s·d with d odd; strongPseudoprime (STRONG-PSEUDOPRIME) is the strong probable-prime condition; Witness is a base that refutes it; and millerRabin is the executable single-base test. (Evaluating millerRabin 561 2 returns false: although 561 is a Carmichael number, base 2 witnesses that it is composite.)

  • Miller-Rabin correctness: strongPseudoprime_of_prime shows a prime is a strong probable prime to every coprime base — the repeated squaring in STRONG-PSEUDOPRIME can only reach 1 through −1 modulo a prime (via modeq_neg_one_of_sq_eq_one, the roots-of-unity fact). Consequently not_witness_of_prime (a prime has no witness) and witness_not_prime (a witness certifies compositeness) hold.

  • Error-bound foundation: strongPseudoprime_pow — every strong probable prime satisfies a^(n−1) ≡ 1 (mod n), so every strong liar lies in the kernel of a ↦ a^(n−1) (the first step toward showing the liars form a subgroup of the units); modeq_pow_two_sub_one is the (n−1)² ≡ 1 (mod n) fact used there.

  • Error-bound infrastructure (Rabin–Monier): nu is the minimum over prime factors p | n of v₂(p−1); goodUnits (S(n)) is the subgroup of units x with x^(2^(ν(n)−1)·t) ∈ {±1} (preimage of {1, −1} under the power map); and liar_mem_goodSet shows every strong liar lies in S(n) via the order-of-element parity lemma two_pow_succ_dvd_orderOf applied modulo each prime divisor.

  • The Miller-Rabin error bound (Theorem 31.39, sharpened to (n−1)/4 by Rabin–Monier): counting |S(n)| via the cyclicity of prime-power unit groups and the CRT, then bounding |S(n)| ≤ (n−1)/4 by the three-case Rabin–Monier analysis. Theorems goodUnits_card_le (the subgroup bound, split into prime power, semiprime, and ≥3-factors cases) and strongLiars_card_le (at most (n−1)/4 strong liars for odd composite n).

  • Random-witness analysis (the MILLER-RABIN error bound): the count of strong-liar bases among 1, …, n-1 is at most (n-1)/4 (strongLiars_nat_card_le), so a uniformly random base errs with probability at most 1/4. The Probability companion proves the repeated-round bound for the actual residue-based execution over a finite product of independent uniform bases (CLRS Theorem 31.39).

Notation:

  • a ≡ b [MOD n] : Nat.ModEq.

  • Nat.totient n : Euler's totient.

The legacy millerRabinLoop pairs semantic decisions with a detached exponentiation budget. Its second component is not a measured decision cost. The Execution companion provides MillerRabinExecution.run, whose Boolean and modular-multiplication counter come from the same residue computation. Parameter decomposition, comparisons, sampling and bit runtime are outside that counter.

namespace CLRSnamespace Chapter31

Fermat's test (CLRS Theorem 31.31). For a prime p and a coprime to p, a^(p−1) ≡ 1 (mod p). Consequently, if some a coprime to n satisfies a^(n−1) ≢ 1 (mod n), then n is not prime.

theorem fermat_test {p a : ℕ} (hp : Nat.Prime p) (hcop : Nat.Coprime a p) : a ^ (p - 1) ≡ 1 [MOD p] := by rw [← Nat.totient_prime hp] exact Nat.ModEq.pow_totient hcop

n is a Fermat pseudoprime to base a: it is composite yet passes the Fermat test a^(n−1) ≡ 1 (mod n).

def fermatPseudoprime (n a : ℕ) : Prop := ¬ Nat.Prime n ∧ a ^ (n - 1) ≡ 1 [MOD n]

PSEUDOPRIME (CLRS §31.8). The executable test: n passes when 2^(n−1) ≡ 1 (mod n).

def pseudoprime (n : ℕ) : Bool := decide (2 ^ (n - 1) ≡ 1 [MOD n])

PSEUDOPRIME is correct on odd primes: every prime n ≠ 2 passes the test. (The base case n = 2 is a known special case handled separately.)

theorem pseudoprime_correct {n : ℕ} (hn : Nat.Prime n) (hn2 : n ≠ 2) : pseudoprime n = true := by unfold pseudoprime rw [decide_eq_true] have hnot : ¬ 2 ∣ n := by intro h2dvd rcases (Nat.Prime.eq_one_or_self_of_dvd hn 2 h2dvd) with h1 | h2 · omega · exact hn2 h2.symm have hcop : Nat.Coprime 2 n := (Nat.Prime.coprime_iff_not_dvd Nat.prime_two).2 hnot exact fermat_test hn hcop

n is a Carmichael number if it is composite and passes the Fermat test a^(n−1) ≡ 1 (mod n) for every a coprime to n (CLRS §31.8). Such numbers fool the Fermat primality test for every base, so the test alone cannot certify primality.

def isCarmichael (n : ℕ) : Prop := ¬ Nat.Prime n ∧ 1 < n ∧ ∀ a : ℕ, Nat.Coprime a n → a ^ (n - 1) ≡ 1 [MOD n]

A Carmichael number is composite.

theorem carmichael_not_prime {n : ℕ} (h : isCarmichael n) : ¬ Nat.Prime n := h.1

A Carmichael number is larger than one.

theorem carmichael_gt_one {n : ℕ} (h : isCarmichael n) : 1 < n := h.2.1

A Carmichael number passes the Fermat test for every base coprime to it.

theorem carmichael_passes_fermat {n a : ℕ} (h : isCarmichael n) (hcop : Nat.Coprime a n) : a ^ (n - 1) ≡ 1 [MOD n] := h.2.2 a hcop

A Carmichael number is a Fermat pseudoprime to every base coprime to it.

theorem carmichael_fermatPseudoprime {n a : ℕ} (h : isCarmichael n) (hcop : Nat.Coprime a n) : fermatPseudoprime n a := ⟨carmichael_not_prime h, carmichael_passes_fermat h hcop⟩

Combining congruences under coprime moduli. If a ≡ b (mod m) and a ≡ b (mod n) with m, n coprime, then a ≡ b (mod m·n). This is the "glue" needed to lift a congruence from pairwise-coprime prime factors to their product (used to verify that 561 is a Carmichael number).

theorem modeq_of_coprime_mul {a b m n : ℕ} (hcop : Nat.Coprime m n) (hm : a ≡ b [MOD m]) (hn : a ≡ b [MOD n]) : a ≡ b [MOD m * n] := by rcases chinese_remainder (n := m) (m := n) (a := b) (b := b) hcop with ⟨x, hx₁, hx₂, huniq⟩ have hxb : x ≡ b [MOD m * n] := huniq b (Nat.ModEq.refl b) (Nat.ModEq.refl b) have hax : x ≡ a [MOD m * n] := huniq a hm hn exact hax.symm.trans hxb

561 is a Carmichael number (CLRS §31.8). Minimality is not asserted by this theorem. It shows the Fermat test can be fooled by a composite integer for every base coprime to it, so PSEUDOPRIME cannot certify primality.

theorem isCarmichael_561 : isCarmichael 561 := by constructor · intro hp have hdiv3 : 3 ∣ 561 := by norm_num rcases (Nat.Prime.eq_one_or_self_of_dvd hp 3 hdiv3) with h1 | h561 · norm_num at h1 · norm_num at h561 · constructor · norm_num · intro a hcop have hcop3 : Nat.Coprime a 3 := (Nat.Coprime.of_dvd_left (by norm_num : 3 ∣ 561) hcop.symm).symm have hcop11 : Nat.Coprime a 11 := (Nat.Coprime.of_dvd_left (by norm_num : 11 ∣ 561) hcop.symm).symm have hcop17 : Nat.Coprime a 17 := (Nat.Coprime.of_dvd_left (by norm_num : 17 ∣ 561) hcop.symm).symm have h3 : a ^ 560 ≡ 1 [MOD 3] := by simpa [← pow_mul] using (fermat_test (p := 3) (by norm_num : Nat.Prime 3) hcop3).pow 280 have h11 : a ^ 560 ≡ 1 [MOD 11] := by simpa [← pow_mul] using (fermat_test (p := 11) (by norm_num : Nat.Prime 11) hcop11).pow 56 have h17 : a ^ 560 ≡ 1 [MOD 17] := by simpa [← pow_mul] using (fermat_test (p := 17) (by norm_num : Nat.Prime 17) hcop17).pow 35 have h33 : a ^ 560 ≡ 1 [MOD 3 * 11] := modeq_of_coprime_mul (by norm_num) h3 h11 have hfull : a ^ 560 ≡ 1 [MOD 3 * 11 * 17] := modeq_of_coprime_mul (by norm_num) h33 h17 simpa [show 3 * 11 * 17 = 561 by norm_num] using hfull

STRONG-PSEUDOPRIME parameters (CLRS §31.8). Write n−1 = 2^s · d with d odd: s is the exponent of 2 in the prime factorization of n−1, and d is the odd part.

def strongTestParams (n : ℕ) : ℕ × ℕ := (Nat.factorization (n - 1) 2, (n - 1) / 2 ^ Nat.factorization (n - 1) 2)

STRONG-PSEUDOPRIME (CLRS §31.8). n is a strong probable prime to base a if, writing n−1 = 2^s·d with d odd, either a^d ≡ 1 (mod n) or a^(2^i·d) ≡ −1 (mod n) for some i < s.

def strongPseudoprime (n a : ℕ) : Prop := let s := (strongTestParams n).1 let d := (strongTestParams n).2 a ^ d ≡ 1 [MOD n] ∨ ∃ i : Fin s, a ^ (2 ^ (i : ℕ) * d) ≡ n - 1 [MOD n]

WITNESS (CLRS §31.8). A base a witnesses that n is composite when n fails the strong-pseudoprime test to base a. A witness certifies ¬ Nat.Prime n.

def Witness (n a : ℕ) : Prop := ¬ strongPseudoprime n a
instance instDecidableStrongPseudoprime (n a : ℕ) : Decidable (strongPseudoprime n a) := by unfold strongPseudoprime infer_instance

MILLER-RABIN (single base, CLRS §31.8). The executable decision procedure returning whether n is a strong probable prime to base a.

def millerRabin (n a : ℕ) : Bool := decide (strongPseudoprime n a)

STRONG-PSEUDOPRIME decomposition (CLRS §31.8). For n ≠ 0, n = 2^s · (n / 2^s) where s is the exponent of 2 in n — i.e. the odd part of n times 2^s recovers n.

theorem strongTestParams_spec (n : ℕ) (hn : n ≠ 0) : n = 2 ^ (n.factorization 2) * (n / 2 ^ (n.factorization 2)) := by have hdvd : 2 ^ (n.factorization 2) ∣ n := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) hn).2 le_rfl exact (Nat.mul_div_cancel' hdvd).symm

Roots of unity modulo a prime. If b^2 ≡ 1 (mod p) for a prime p and b ≢ 1 (mod p), then b ≡ −1 (mod p). This is the only fact about prime moduli needed by the Miller-Rabin correctness proof: repeatedly squaring a root of unity can only reach 1 through −1.

theorem modeq_neg_one_of_sq_eq_one {p b : ℕ} (hp : Nat.Prime p) (hb : 1 ≤ b) (hb2 : b ^ 2 ≡ 1 [MOD p]) (hbne : ¬ b ≡ 1 [MOD p]) : b ≡ p - 1 [MOD p] := by have hb2' : 1 ≡ b ^ 2 [MOD p] := hb2.symm have hb2_dvd : p ∣ b ^ 2 - 1 := by exact (Nat.modEq_iff_dvd' (by nlinarith : 1 ≤ b ^ 2)).mp hb2' have hsq : b ^ 2 - 1 = (b - 1) * (b + 1) := by have hsub := Nat.sq_sub_sq b 1 simpa [mul_comm] using hsub rw [hsq] at hb2_dvd have hdvd_or := (hp.dvd_mul).1 hb2_dvd rcases hdvd_or with hd1 | hd2 · exfalso exact hbne ((Nat.modEq_iff_dvd' hb).mpr hd1).symm · have h0 : b + 1 ≡ 0 [MOD p] := hd2.modEq_zero_nat have h1 : (p - 1) + 1 ≡ 0 [MOD p] := by rw [Nat.sub_add_cancel (Nat.one_le_of_lt hp.pos)] exact Nat.modEq_zero_iff_dvd.mpr (dvd_refl p) have h : b + 1 ≡ (p - 1) + 1 [MOD p] := h0.trans h1.symm exact Nat.ModEq.add_right_cancel (Nat.ModEq.refl 1) h

Miller-Rabin correctness, prime direction. If n is prime and a is coprime to n, then n is a strong probable prime to base a — i.e. a prime has no witness. Writing n−1 = 2^s·d, the sequence a^d, a^{2d}, …, a^{2^s·d} ends at 1 (Fermat); at the first index where it reaches 1, the previous value is a square root of 1 that is not 1, hence −1 (mod a prime).

theorem strongPseudoprime_of_prime {n a : ℕ} (hn : Nat.Prime n) (hcop : Nat.Coprime a n) : strongPseudoprime n a := by have hn2 : 2 ≤ n := hn.two_le have hnm1 : n - 1 ≠ 0 := by omega have hdecomp : n - 1 = 2 ^ (strongTestParams n).1 * (strongTestParams n).2 := by have h := strongTestParams_spec (n - 1) hnm1 simpa [strongTestParams] using h have hfermat : a ^ (n - 1) ≡ 1 [MOD n] := fermat_test hn hcop let s := (strongTestParams n).1 let d := (strongTestParams n).2 have h1 : a ^ (2 ^ s * d) ≡ 1 [MOD n] := by rw [← hdecomp] exact hfermat let P : ℕ → Prop := fun i => a ^ (2 ^ i * d) ≡ 1 [MOD n] have hP_s : P s := by simpa [P, hdecomp] using hfermat let i0 := Nat.find ⟨s, hP_s⟩ have hP0 : P i0 := by exact Nat.find_spec ⟨s, hP_s⟩ have hmin : ∀ m, m < i0 → ¬ P m := by intro m hm exact Nat.find_min ⟨s, hP_s⟩ (by simpa [i0] using hm) have hile : i0 ≤ s := by simpa [i0] using (Nat.find_le (h := ⟨s, hP_s⟩) hP_s) by_cases hi00 : i0 = 0 · left have : P 0 := by simpa [hi00] using hP0 simpa [P] using this · right let j := i0 - 1 have hj : j < s := by omega have hjlt_i0 : j < i0 := by omega have hj1 : j + 1 = i0 := by omega have hb : a ^ (2 ^ j * d) ≡ n - 1 [MOD n] := by have hb2 : (a ^ (2 ^ j * d)) ^ 2 ≡ 1 [MOD n] := by have hsq : (a ^ (2 ^ j * d)) ^ 2 = a ^ (2 ^ i0 * d) := by rw [← pow_mul] congr 1 rw [mul_assoc, mul_comm d 2, ← mul_assoc, ← pow_succ, hj1] rw [hsq] simpa [P] using hP0 have hbne : ¬ a ^ (2 ^ j * d) ≡ 1 [MOD n] := by exact hmin j hjlt_i0 have ha1 : 1 ≤ a := by have ha_ne : a ≠ 0 := by intro ha have hg : Nat.gcd a n = 1 := hcop have hn1 : n = 1 := by rw [ha, Nat.gcd_zero_left] at hg exact hg omega exact Nat.succ_le_of_lt (Nat.pos_of_ne_zero ha_ne) have hb1 : 1 ≤ a ^ (2 ^ j * d) := by exact one_le_pow₀ ha1 exact modeq_neg_one_of_sq_eq_one hn hb1 hb2 hbne exact ⟨⟨j, hj⟩, hb⟩

A prime has no witness: for n prime and a coprime to n, a does not witness compositeness of n.

theorem not_witness_of_prime {n a : ℕ} (hn : Nat.Prime n) (hcop : Nat.Coprime a n) : ¬ Witness n a := by intro hw exact hw (strongPseudoprime_of_prime hn hcop)

A witness certifies compositeness: for n and a coprime base a, if a is a witness then n is not prime.

theorem witness_not_prime {n a : ℕ} (hcop : Nat.Coprime a n) (hw : Witness n a) : ¬ Nat.Prime n := by intro hn exact (not_witness_of_prime hn hcop) hw

(n−1)² ≡ 1 (mod n) for n ≠ 0: n divides (n−1)² − 1 = n·(n−2).

theorem modeq_pow_two_sub_one (hn : n ≠ 0) : (n - 1) ^ 2 ≡ 1 [MOD n] := by rw [Nat.ModEq] by_cases h2 : 2 ≤ n · have hsq2 := Nat.sq_sub_sq (n - 1) 1 have hsq2' : (n - 1) ^ 2 - 1 = (n - 1 + 1) * (n - 1 - 1) := by simpa using hsq2 have hform : (n - 1) ^ 2 - 1 = n * (n - 2) := by rw [hsq2'] have h1 : (n - 1) + 1 = n := by omega have h2' : (n - 1) - 1 = n - 2 := by omega rw [h1, h2', mul_comm] have hdvd : n ∣ (n - 1) ^ 2 - 1 := by rw [hform] exact dvd_mul_right n (n - 2) rcases hdvd with ⟨k, hk⟩ have hle : 1 ≤ (n - 1) ^ 2 := by have ht : 1 ≤ n - 1 := by omega nlinarith have hk' : (n - 1) ^ 2 = n * k + 1 := by rw [← hk] omega rw [hk'] simp · have hn1 : n = 1 := by omega simp [hn1]

A strong probable prime satisfies Fermat's congruence. If n is a strong pseudoprime to base a, then a^(n−1) ≡ 1 (mod n). Indeed n−1 = 2^s·d, and either a^d ≡ 1 or a^(2^i·d) ≡ −1 for some i < s; in the second case a^(n−1) = (a^(2^i·d))^(2^(s−i)) ≡ (−1)^even = 1. This is the first step of the Miller-Rabin error-bound proof: every strong liar lies in the kernel of a ↦ a^(n−1).

theorem strongPseudoprime_pow {n a : ℕ} (h : strongPseudoprime n a) : a ^ (n - 1) ≡ 1 [MOD n] := by by_cases hn0 : n - 1 = 0 · simpa [hn0] using (Nat.ModEq.refl 1) · have hdecomp : n - 1 = 2 ^ (strongTestParams n).1 * (strongTestParams n).2 := by have hd := strongTestParams_spec (n - 1) hn0 simpa [strongTestParams] using hd unfold strongPseudoprime at h rw [hdecomp] rcases h with hd1 | ⟨i, hi⟩ · have hpow := hd1.pow (2 ^ (strongTestParams n).1) simpa [← pow_mul, mul_comm] using hpow · have hpow := hi.pow (2 ^ ((strongTestParams n).1 - (i : ℕ))) have hsum : (i : ℕ) + ((strongTestParams n).1 - (i : ℕ)) = (strongTestParams n).1 := by exact Nat.add_sub_of_le (Nat.le_of_lt i.isLt) have hpowsum : 2 ^ (i : ℕ) * 2 ^ ((strongTestParams n).1 - (i : ℕ)) = 2 ^ (strongTestParams n).1 := by rw [← pow_add, hsum] have hexp : (2 ^ (i : ℕ) * (strongTestParams n).2) * 2 ^ ((strongTestParams n).1 - (i : ℕ)) = 2 ^ (strongTestParams n).1 * (strongTestParams n).2 := by rw [mul_assoc] rw [mul_comm (strongTestParams n).2 (2 ^ ((strongTestParams n).1 - (i : ℕ)))] rw [← mul_assoc, hpowsum] have h1 : a ^ (2 ^ (strongTestParams n).1 * (strongTestParams n).2) = (a ^ (2 ^ (i : ℕ) * (strongTestParams n).2)) ^ (2 ^ ((strongTestParams n).1 - (i : ℕ))) := by rw [← pow_mul] congr 1 exact hexp.symm rw [h1] have h2 : (n - 1) ^ (2 ^ ((strongTestParams n).1 - (i : ℕ))) ≡ 1 [MOD n] := by have hsmi_pos : 0 < (strongTestParams n).1 - (i : ℕ) := by exact Nat.sub_pos_of_lt i.isLt have heven : 2 ^ ((strongTestParams n).1 - (i : ℕ)) = 2 * 2 ^ (((strongTestParams n).1 - (i : ℕ)) - 1) := by rw [mul_comm, ← pow_succ, Nat.sub_add_cancel hsmi_pos] rw [heven] rw [pow_mul] have hsq := modeq_pow_two_sub_one (n := n) (by omega : n ≠ 0) simpa using (hsq.pow (2 ^ (((strongTestParams n).1 - (i : ℕ)) - 1))) exact hpow.trans h2

Error bound: the good subgroup S(n) (Rabin–Monier)

The Miller-Rabin error bound (Theorem 31.39; sharpened to (n−1)/4 by Rabin and Monier) states that for odd composite n, at most (n−1)/4 of the bases are strong liars. The proof embeds the liars into a subgroup S(n) of the units modulo n and bounds |S(n)|. This section develops the infrastructure: the units of ZMod n, the cyclicity of prime-power unit groups (from Mathlib), the 2-adic valuation ν(n), and the good subgroup S(n).

The reduction (ZMod n)ˣ → (ZMod p)ˣ of units when p ∣ n.

def zmodUnitReduction {n p : ℕ} (hp : p ∣ n) : (ZMod n)ˣ →* (ZMod p)ˣ := Units.map (ZMod.castHom hp (ZMod p)).toMonoidHom

The unit group of ZMod n has φ(n) elements.

theorem units_card_eq_totient {n : ℕ} [NeZero n] : Nat.card (ZMod n)ˣ = Nat.totient n := by rw [Nat.card_eq_fintype_card] exact ZMod.card_units_eq_totient n

For a prime p, the unit group of ZMod p has p − 1 elements.

theorem units_card_prime {p : ℕ} (hp : Nat.Prime p) : Nat.card (ZMod p)ˣ = p - 1 := by haveI : NeZero p := ⟨hp.ne_zero⟩ rw [Nat.card_eq_fintype_card, ZMod.card_units_eq_totient p, Nat.totient_prime hp]

The subgroup {1, −1} of the units modulo n.

def negOneTwoSubgroup {n : ℕ} [NeZero n] : Subgroup (ZMod n)ˣ where carrier := {x : (ZMod n)ˣ | x = 1 ∨ x = -1} one_mem' := by simp mul_mem' := by intro a b ha hb rcases ha with ha1 | ham · rcases hb with hb1 | hbm · simp [ha1, hb1] · simp [ha1, hbm] · rcases hb with hb1 | hbm · simp [ham, hb1] · simp [ham, hbm] inv_mem' := by intro a ha rcases ha with ha1 | ham · simp [ha1] · simp [ham]

The power homomorphism x ↦ x^m on the units modulo n.

def unitPowHom {n : ℕ} [NeZero n] (m : ℕ) : (ZMod n)ˣ →* (ZMod n)ˣ where toFun x := x ^ m map_one' := by simp map_mul' := by intro a b simp [mul_pow]

The good set S_m = {x ∈ (ZMod n)ˣ : x^m ∈ {1, −1}}, a subgroup of the units (the preimage of {1, −1} under the power map x ↦ x^m).

def goodSet {n : ℕ} [NeZero n] (m : ℕ) : Subgroup (ZMod n)ˣ := negOneTwoSubgroup.comap (unitPowHom (n := n) m)

Membership in goodSet as a power-of-a-unit condition.

theorem mem_goodSet_iff {n : ℕ} [NeZero n] (m : ℕ) (x : (ZMod n)ˣ) : x ∈ goodSet (n := n) m ↔ x ^ m = 1 ∨ x ^ m = -1 := by rfl

Membership in goodSet as a power condition in ZMod n.

theorem mem_goodSet_zmod {n : ℕ} [NeZero n] (m : ℕ) (x : (ZMod n)ˣ) : x ∈ goodSet (n := n) m ↔ (x : ZMod n) ^ m = 1 ∨ (x : ZMod n) ^ m = -1 := by rw [mem_goodSet_iff] constructor · rintro (h | h) · left exact congrArg (fun u : (ZMod n)ˣ => (u : ZMod n)) h · right exact congrArg (fun u : (ZMod n)ˣ => (u : ZMod n)) h · rintro (h | h) · left exact Units.ext h · right exact Units.ext h

ν(n): the minimum over prime factors p of n of v₂(p − 1), the 2-adic valuation of p − 1. This is the index bound used to define the good subgroup for the Miller-Rabin error bound.

noncomputable def nu (n : ℕ) : ℕ := if h : n.primeFactors.Nonempty then (n.primeFactors.image (fun p => (p - 1).factorization 2)).min' (h.image (fun p => (p - 1).factorization 2)) else 0

For every prime factor p | n, ν(n) ≤ v₂(p − 1).

theorem nu_le_v2 (hn : n.primeFactors.Nonempty) {p : ℕ} (hp : p ∈ n.primeFactors) : nu n ≤ (p - 1).factorization 2 := by unfold nu rw [dif_pos hn] have hmem : (p - 1).factorization 2 ∈ (n.primeFactors.image (fun p => (p - 1).factorization 2)) := by exact Finset.mem_image.mpr ⟨p, hp, rfl⟩ exact (Finset.isLeast_min' _ _).2 (by simpa using hmem)

2^ν(n) divides p − 1 for every prime factor p | n.

theorem two_pow_nu_dvd_prime_sub_one (hn : n.primeFactors.Nonempty) {p : ℕ} (hp : p ∈ n.primeFactors) : 2 ^ nu n ∣ p - 1 := by have hp' : Nat.Prime p := Nat.prime_of_mem_primeFactors hp have hppos : p - 1 ≠ 0 := by have := hp'.two_le omega exact (Nat.Prime.pow_dvd_iff_le_factorization (p := 2) (k := nu n) (n := p - 1) (by decide : Nat.Prime 2) hppos).2 (nu_le_v2 hn hp)

ν(n) ≥ 1 for odd n > 1: every prime factor is odd, so p − 1 is even.

theorem nu_pos {n : ℕ} (hn_odd : Odd n) (hn1 : 1 < n) : 1 ≤ nu n := by have hnne : n.primeFactors.Nonempty := (Nat.nonempty_primeFactors).2 hn1 unfold nu rw [dif_pos hnne] apply Finset.le_min' intro y hy rcases Finset.mem_image.mp hy with ⟨p, hp, rfl⟩ have hp' : Nat.Prime p := Nat.prime_of_mem_primeFactors hp have hpdvd : p ∣ n := Nat.dvd_of_mem_primeFactors hp have hpodd : Odd p := by have hne_even_n : ¬ Even n := (Nat.not_even_iff_odd.mpr hn_odd) have hne_even_p : ¬ Even p := by intro hep rcases hep with ⟨k, hk⟩ rcases hpdvd with ⟨m, hm⟩ refine hne_even_n ⟨k * m, ?_⟩ rw [hm, hk] ring exact (Nat.not_even_iff_odd.mp hne_even_p) have h2dvd : 2 ∣ p - 1 := by rcases (odd_iff_exists_bit1.mp hpodd) with ⟨k, rfl⟩ exact ⟨k, by omega⟩ have hppos : p - 1 ≠ 0 := by have := hp'.two_le omega exact (Nat.Prime.pow_dvd_iff_le_factorization (p := 2) (k := 1) (n := p - 1) (by decide : Nat.Prime 2) hppos).1 h2dvd

The odd part a / 2^(v₂(a)) of a nonzero a is odd.

lemma oddPart_odd (a : ℕ) (ha : a ≠ 0) : Odd (a / 2 ^ a.factorization 2) := by have hdvd : 2 ^ a.factorization 2 ∣ a := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) ha).2 le_rfl have hfac : (a / 2 ^ a.factorization 2).factorization 2 = 0 := by rw [Nat.factorization_div hdvd] simp [This simp argument is unused: Nat.factorization_pow_self (by decide : Nat.Prime 2) Hint: Omit it from the simp argument list. simp [̵N̵a̵t̵.̵f̵a̵c̵t̵o̵r̵i̵z̵a̵t̵i̵o̵n̵_̵p̵o̵w̵_̵s̵e̵l̵f̵[̲N̲a̲t̲.̲P̲r̲i̲m̲e̲.̲f̲a̲c̲t̲o̲r̲i̲z̲a̲t̲i̲o̲n̲_̲s̲e̲l̲f̲ (by decide : Nat.Prime 2̵)̵,̵ ̵ ̵ ̵ ̵ ̵ ̵ ̵N̵a̵t̵.̵P̵r̵i̵m̵e̵.̵f̵a̵c̵t̵o̵r̵i̵z̵a̵t̵i̵o̵n̵_̵s̵e̵l̵f̵ ̵(̵b̵y̵ ̵d̵e̵c̵i̵d̵e̵ ̵:̵ ̵N̵a̵t̵.̵P̵r̵i̵m̵e̵ ̵2)] Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Nat.factorization_pow_self (by decide : Nat.Prime 2), Nat.Prime.factorization_self (by decide : Nat.Prime 2)] have htpos : 0 < a / 2 ^ a.factorization 2 := by exact Nat.div_pos (Nat.le_of_dvd (Nat.pos_of_ne_zero ha) hdvd) (pow_pos (by norm_num) _) have h2not : ¬ 2 ∣ a / 2 ^ a.factorization 2 := by intro h2 have hle1 : 1 ≤ (a / 2 ^ a.factorization 2).factorization 2 := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) htpos.ne').1 h2 omega exact (Nat.not_even_iff_odd.mp (even_iff_two_dvd.not.mpr h2not))

The odd part t = (n−1)/2^(v₂(n−1)) of n−1 is odd.

theorem strongTestParams_odd {n : ℕ} (hn1 : 1 < n) : Odd (strongTestParams n).2 := by unfold strongTestParams exact oddPart_odd (n - 1) (by omega)

The good subgroup S(n) (Rabin–Monier). Writing n−1 = 2^s·t with t odd and ν = ν(n), S(n) = {x ∈ (ZMod n)ˣ : x^(2^(ν−1)·t) ∈ {1, −1}}. Every strong liar lies in S(n), and |S(n)| ≤ (n−1)/4.

noncomputable def goodUnits {n : ℕ} [NeZero n] : Subgroup (ZMod n)ˣ := goodSet (2 ^ (nu n - 1) * (strongTestParams n).2)

Parity lemma. If a^(2^i·d) has order 2 in a finite group and d is odd, then 2^(i+1) divides the order of a. Indeed orderOf a = 2·gcd(orderOf a, 2^i·d), whose 2-adic valuation forces v₂(orderOf a) = i+1.

theorem two_pow_succ_dvd_orderOf {G : Type*} [Group G] [Finite G] {a : G} {i d : ℕ} (hd : Odd d) (hord : orderOf (a ^ (2 ^ i * d)) = 2) : 2 ^ (i + 1) ∣ orderOf a := by let e := orderOf a let g := e.gcd (2 ^ i * d) have hpow := orderOf_pow (x := a) (n := 2 ^ i * d) have hdiv : e / g = 2 := by dsimp [e, g] rw [← hpow] exact hord have he : e = 2 * g := by have hg : g ∣ e := by dsimp [g] exact Nat.gcd_dvd_left _ _ have h := Nat.mul_div_cancel' hg rw [hdiv] at h rw [mul_comm] at h exact h.symm have epos : e ≠ 0 := by dsimp [e] exact (orderOf_pos a).ne' have dpos : d ≠ 0 := by intro hz rw [hz] at hd norm_num at hd have h2dpos : 2 ^ i * d ≠ 0 := mul_ne_zero (pow_ne_zero i (by norm_num)) dpos have hgpos : g ≠ 0 := by intro hg0 have hg_dvd : g ∣ e := by dsimp [g] exact Nat.gcd_dvd_left _ _ rw [hg0] at hg_dvd exact epos (zero_dvd_iff.mp hg_dvd) have h2notd : ¬ 2 ∣ d := (even_iff_two_dvd.not.mp (Nat.not_even_iff_odd.mpr hd)) have hpowfac : (2 ^ i * d).factorization 2 = i := by rw [Nat.factorization_mul (pow_ne_zero i (by norm_num)) dpos] simp [Nat.Prime.factorization_self (by decide : Nat.Prime 2), Nat.factorization_eq_zero_of_not_dvd h2notd] have hgfac : g.factorization 2 = min (e.factorization 2) i := by dsimp [g] rw [Nat.factorization_gcd epos h2dpos] simp [hpowfac] have hfac : e.factorization 2 = 1 + min (e.factorization 2) i := by calc e.factorization 2 = (2 * g).factorization 2 := by rw [he] _ = (2 : ℕ).factorization 2 + g.factorization 2 := by rw [Nat.factorization_mul (by norm_num) hgpos] simp _ = 1 + g.factorization 2 := by rw [Nat.Prime.factorization_self (by decide : Nat.Prime 2)] _ = 1 + min (e.factorization 2) i := by rw [hgfac] have hi_lt : i < e.factorization 2 := by by_contra hnot have hle : e.factorization 2 ≤ i := Nat.not_lt.mp hnot have hz : min (e.factorization 2) i = e.factorization 2 := min_eq_left hle have hcontra : e.factorization 2 = 1 + e.factorization 2 := by calc e.factorization 2 = 1 + min (e.factorization 2) i := hfac _ = 1 + e.factorization 2 := by rw [hz] omega change 2 ^ (i + 1) ∣ e exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) epos).2 (Nat.succ_le_of_lt hi_lt)

In (ZMod p)ˣ for an odd prime p, the element −1 has order 2.

theorem orderOf_neg_one {p : ℕ} (hp : Nat.Prime p) (hp2 : p ≠ 2) : orderOf (-1 : (ZMod p)ˣ) = 2 := by haveI : Fact p.Prime := ⟨hp⟩ exact orderOf_eq_prime (by simp) (by intro h have hz : (-1 : ZMod p) = 1 := by simpa using congrArg Units.val h haveI : Fact (2 < p) := ⟨lt_of_le_of_ne hp.two_le (Ne.symm hp2)⟩ exact ZMod.neg_one_ne_one hz)

For an odd prime p, if x^(2^i·d) = −1 with d odd, then 2^(i+1) divides p − 1: x^(2^i·d) has order 2, so by the parity lemma 2^(i+1) | orderOf x, and orderOf x | p−1 by Lagrange.

theorem two_pow_succ_dvd_prime_sub_one {p : ℕ} (hp : Nat.Prime p) (hp2 : p ≠ 2) {x : (ZMod p)ˣ} {i d : ℕ} (hd : Odd d) (hx : (x : ZMod p) ^ (2 ^ i * d) = -1) : 2 ^ (i + 1) ∣ p - 1 := by haveI : NeZero p := ⟨hp.ne_zero⟩ have hord : orderOf (x ^ (2 ^ i * d)) = 2 := by have hx' : x ^ (2 ^ i * d) = -1 := by apply Units.ext simp [hx] rw [hx'] exact orderOf_neg_one hp hp2 have hdvd_ord : 2 ^ (i + 1) ∣ orderOf x := two_pow_succ_dvd_orderOf hd hord have hcard : Nat.card (ZMod p)ˣ = p - 1 := units_card_prime hp have hord_dvd : orderOf x ∣ p - 1 := (orderOf_dvd_natCard x).trans (by rw [hcard]) exact hdvd_ord.trans hord_dvd

The strong-liar predicate on units modulo n: n is a strong probable prime to base a (equivalently strongPseudoprime n (a : ZMod n).val).

def isStrongLiar {n : ℕ} [NeZero n] (a : (ZMod n)ˣ) : Prop := (a : ZMod n) ^ (strongTestParams n).2 = 1 ∨ ∃ i : Fin (strongTestParams n).1, (a : ZMod n) ^ (2 ^ (i : ℕ) * (strongTestParams n).2) = -1

Every strong liar lies in the good subgroup. For odd n > 1, if a is a strong liar modulo n, then a^(2^(ν(n)−1)·t) ∈ {±1}, i.e. a ∈ S(n). The index bound i < ν(n) is obtained by reducing modulo each prime divisor p of n and applying the parity lemma, which forces 2^(i+1) | p−1.

theorem liar_mem_goodSet {n : ℕ} [NeZero n] (hn_odd : Odd n) (hn1 : 1 < n) {a : (ZMod n)ˣ} (hliar : isStrongLiar a) : a ∈ goodUnits := by rw [goodUnits, mem_goodSet_zmod] rcases hliar with hd1 | ⟨i, hi⟩ · left have h : (a : ZMod n) ^ (2 ^ (nu n - 1) * (strongTestParams n).2) = 1 := by rw [mul_comm] rw [pow_mul] rw [hd1] simp exact h · have hv2 : ∀ p, p ∈ n.primeFactors → 2 ^ ((i : ℕ) + 1) ∣ p - 1 := by intro p hp have hpp : Nat.Prime p := Nat.prime_of_mem_primeFactors hp have hp_dvd : p ∣ n := Nat.dvd_of_mem_primeFactors hp have hpne2 : p ≠ 2 := by intro hp2 have hne : ¬ Even n := (Nat.not_even_iff_odd.mpr hn_odd) exact hne (even_iff_two_dvd.mpr (by rw [← hp2]; exact hp_dvd)) let φ := zmodUnitReduction (n := n) (p := p) hp_dvd have hφ : ((φ a : (ZMod p)ˣ) : ZMod p) ^ (2 ^ (i : ℕ) * (strongTestParams n).2) = -1 := by have hc := congrArg (ZMod.castHom hp_dvd (ZMod p)) hi have hc' : (ZMod.castHom hp_dvd (ZMod p) (a : ZMod n)) ^ (2 ^ (i : ℕ) * (strongTestParams n).2) = -1 := by simpa [map_pow, map_neg, map_one] using hc have hval : (φ a : ZMod p) = ZMod.castHom hp_dvd (ZMod p) (a : ZMod n) := by simp [φ, zmodUnitReduction] rwa [hval] exact two_pow_succ_dvd_prime_sub_one hpp hpne2 (strongTestParams_odd (n := n) hn1) hφ have hle : (i : ℕ) + 1 ≤ nu n := by have hnne : n.primeFactors.Nonempty := (Nat.nonempty_primeFactors).2 hn1 have hfac : ∀ p, p ∈ n.primeFactors → (i : ℕ) + 1 ≤ (p - 1).factorization 2 := by intro p hp have hpp : Nat.Prime p := Nat.prime_of_mem_primeFactors hp exact (Nat.Prime.pow_dvd_iff_le_factorization (p := 2) (k := (i : ℕ) + 1) (n := p - 1) (by decide : Nat.Prime 2) (by have := hpp.two_le; omega)).1 (hv2 p hp) unfold nu rw [dif_pos hnne] apply Finset.le_min' intro y hy rcases Finset.mem_image.mp hy with ⟨p, hp, rfl⟩ exact hfac p hp have hi_lt : (i : ℕ) < nu n := by omega have hsum : (i : ℕ) + (nu n - 1 - (i : ℕ)) = nu n - 1 := by omega have hexp : (a : ZMod n) ^ (2 ^ (nu n - 1) * (strongTestParams n).2) = ((a : ZMod n) ^ (2 ^ (i : ℕ) * (strongTestParams n).2)) ^ (2 ^ (nu n - 1 - (i : ℕ))) := by have hexp_eq : 2 ^ (nu n - 1) * (strongTestParams n).2 = (2 ^ (i : ℕ) * (strongTestParams n).2) * 2 ^ (nu n - 1 - (i : ℕ)) := by rw [mul_assoc, mul_comm (strongTestParams n).2 (2 ^ (nu n - 1 - (i : ℕ))), ← mul_assoc, ← pow_add, hsum] rw [hexp_eq, pow_mul] rw [hexp, hi] by_cases h0 : nu n - 1 - (i : ℕ) = 0 · right rw [h0] simp · left have hpos : 0 < nu n - 1 - (i : ℕ) := Nat.pos_of_ne_zero h0 have hk : 2 ^ (nu n - 1 - (i : ℕ)) = 2 * 2 ^ ((nu n - 1 - (i : ℕ)) - 1) := by rw [mul_comm, ← pow_succ, Nat.sub_add_cancel hpos] rw [hk, pow_mul] simp

Counting |S(n)| — cyclic torsion counts

The Rabin-Monier bound needs the cardinality of S(n), which is a product of per-prime-power counts of solutions to x^m ≡ ±1. Each such count is a gcd in a cyclic group; the lemmas below provide that counting primitive for a finite cyclic group.

The number of multiples of d in [0, N) is N/d when d | N.

try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false` lemma card_multiples_dvd {N d : ℕ} (hd : d ∣ N) : Nat.card {i : Fin N // d ∣ (i : ℕ)} = N / d := by by_cases hN : N = 0 · subst hN simp · have hd0 : d ≠ 0 := by intro hz apply hN exact (zero_dvd_iff.mp (by simpa [hz] using hd)) have hdpos : 0 < d := Nat.pos_of_ne_zero hd0 have hN_eq : N = d * (N / d) := (Nat.mul_div_cancel' hd).symm let e : {i : Fin N // d ∣ (i : ℕ)} ≃ Fin (N / d) := { toFun := fun i => ⟨i.val.val / d, by have hx : i.val.val < N := i.val.isLt have hx' : i.val.val < d * (N / d) := lt_of_lt_of_eq hx hN_eq exact (Nat.div_lt_iff_lt_mul (k := d) (x := i.val.val) (y := N / d) hdpos).mpr (by simpa [mul_comm] using hx')⟩ invFun := fun k => ⟨⟨k.val * d, by have hk : k.val < N / d := k.isLt have hk' : k.val * d < (N / d) * d := Nat.mul_lt_mul_of_pos_right hk hdpos have : (N / d) * d = N := by simpa [mul_comm] using hN_eq.symm rw [this] at hk' exact hk'⟩, by try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [mul_comm] using (dvd_mul_right d k.val : d ∣ d * k.val)⟩ left_inv := by intro i apply Subtype.ext apply Fin.ext simp only [This simp argument is unused: Fin.val_mk Hint: Omit it from the simp argument list. simp only ̵[̵F̵i̵n̵.̵v̵a̵l̵_̵m̵k̵]̵ Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`Fin.val_mk] exact (mul_comm (i.val.val / d) d).trans (Nat.mul_div_cancel' i.property) right_inv := by intro k apply Fin.ext simpa [mul_comm] using (Nat.mul_div_right k.val hdpos) } rw [Nat.card_congr e] simp

The number of i < N with N | i·n is gcd(N, n).

lemma card_fin_dvd_mul {N n : ℕ} (hN : N ≠ 0) : Nat.card {i : Fin N // N ∣ (i : ℕ) * n} = N.gcd n := by let g := N.gcd n have hgN : g ∣ N := Nat.gcd_dvd_left N n have hgn : g ∣ n := Nat.gcd_dvd_right N n have hg0 : g ≠ 0 := by intro hz apply hN exact (Nat.gcd_eq_zero_iff.mp hz).1 have hgpos : 0 < g := Nat.pos_of_ne_zero hg0 have hN_eq : N = g * (N / g) := (Nat.mul_div_cancel' hgN).symm have hn_eq : n = g * (n / g) := (Nat.mul_div_cancel' hgn).symm have hcop : Nat.Coprime (N / g) (n / g) := by unfold Nat.Coprime rw [Nat.gcd_div hgN hgn] dsimp [g] nth_rw 1 [← Nat.mul_one (N.gcd n)] exact Nat.mul_div_right 1 (Nat.pos_of_ne_zero hg0) have hiff : ∀ i : Fin N, (N ∣ (i : ℕ) * n) ↔ ((N / g) ∣ (i : ℕ)) := by intro i constructor · intro h have hN : g * (N / g) = N := Nat.mul_div_cancel' hgN have hn' : g * ((i : ℕ) * (n / g)) = (i : ℕ) * n := by calc g * ((i : ℕ) * (n / g)) = (i : ℕ) * (g * (n / g)) := by ring _ = (i : ℕ) * n := by rw [← hn_eq] have h1 : g * (N / g) ∣ g * ((i : ℕ) * (n / g)) := by rwa [hN, hn'] have hdiv : N / g ∣ (i : ℕ) * (n / g) := (Nat.mul_dvd_mul_iff_left hgpos).mp h1 exact hcop.dvd_of_dvd_mul_left (by simpa [mul_comm] using hdiv) · intro h rcases h with ⟨c, hc⟩ have hNdvd : N ∣ (N / g) * n := by use n / g calc (N / g) * n = (N / g) * (g * (n / g)) := by rw [← hn_eq] _ = ((N / g) * g) * (n / g) := by rw [mul_assoc] _ = N * (n / g) := by have : (N / g) * g = N := by rw [mul_comm] exact Nat.mul_div_cancel' hgN rw [this] have hc' : (i : ℕ) * n = c * ((N / g) * n) := by rw [hc] ring rw [hc'] simpa [mul_assoc, mul_comm, mul_left_comm] using (dvd_mul_of_dvd_left hNdvd c) have hcount : Nat.card {i : Fin N // (N / g) ∣ (i : ℕ)} = N / (N / g) := card_multiples_dvd (d := N / g) (by use g rw [mul_comm] exact (Nat.mul_div_cancel' hgN).symm) have hquot : N / (N / g) = g := by have hNg : (N / g) * g = N := by rw [mul_comm] exact Nat.mul_div_cancel' hgN have hpos : 0 < N / g := by exact Nat.div_pos (Nat.le_of_dvd (Nat.pos_of_ne_zero hN) hgN) hgpos calc N / (N / g) = ((N / g) * g) / (N / g) := by nth_rw 1 [← hNg] _ = g := by exact Nat.mul_div_right g hpos let e : {i : Fin N // N ∣ (i : ℕ) * n} ≃ {i : Fin N // N / g ∣ (i : ℕ)} := { toFun := fun i => ⟨i.1, (hiff i.1).mp i.2⟩ invFun := fun i => ⟨i.1, (hiff i.1).mpr i.2⟩ left_inv := by intro i; apply Subtype.ext; rfl right_inv := by intro i; apply Subtype.ext; rfl } rw [Nat.card_congr e] rw [hcount, hquot]

In a finite cyclic group of order N, the number of elements with x^n = 1 is gcd(n, N). This is the counting primitive for the Rabin-Monier bound: the number of solutions to x^m ≡ 1 modulo a prime power.

theorem card_pow_eq_one_cyclic {α : Type*} [Group α] [Fintype α] [DecidableEq α] [IsCyclic α] (n : ℕ) : Nat.card {x : α // x ^ n = 1} = Nat.gcd n (Fintype.card α) := by let N := Nat.card α obtain ⟨g, hg⟩ := IsCyclic.exists_generator (α := α) have horder : orderOf g = N := by dsimp [N] exact orderOf_eq_card_of_forall_mem_zpowers hg let f : Fin N → α := fun i => g ^ (i : ℕ) have hf_inj : Function.Injective f := by intro i j hij apply Fin.ext have hmod : (i : ℕ) ≡ (j : ℕ) [MOD orderOf g] := by exact (pow_eq_pow_iff_modEq).1 hij rw [horder] at hmod have hi : (i : ℕ) < N := i.isLt have hj : (j : ℕ) < N := j.isLt exact Nat.ModEq.eq_of_lt_of_lt hmod hi hj have hf_surj : Function.Surjective f := by intro x have hx : x ∈ (Finset.range N).image (fun i : ℕ => g ^ i) := by dsimp [N] rw [IsCyclic.image_range_card hg] exact Finset.mem_univ x rcases Finset.mem_image.mp hx with ⟨k, hk, rfl⟩ exact ⟨⟨k, (Finset.mem_range.mp hk)⟩, rfl⟩ let ef := Equiv.ofBijective f ⟨hf_inj, hf_surj⟩ have hN : N ≠ 0 := by dsimp [N] exact Nat.card_pos.ne' have h1 : Nat.card {i : Fin N // (f i) ^ n = 1} = Nat.card {x : α // x ^ n = 1} := by exact Nat.card_congr (Equiv.subtypeEquiv ef (by intro i; rfl)) have h2 : Nat.card {i : Fin N // (f i) ^ n = 1} = Nat.card {i : Fin N // N ∣ (i : ℕ) * n} := by apply Nat.card_congr apply Equiv.subtypeEquiv (Equiv.refl (Fin N)) intro i have hpow : (f i) ^ n = g ^ ((i : ℕ) * n) := by simp [f, ← pow_mul, mul_comm, This simp argument is unused: mul_left_comm Hint: Omit it from the simp argument list. simp [f, ← pow_mul, mul_comm,̵ ̵m̵u̵l̵_̵l̵e̵f̵t̵_̵c̵o̵m̵m̵] Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`mul_left_comm] constructor · intro h rw [hpow] at h have hdvd : orderOf g ∣ (i : ℕ) * n := orderOf_dvd_iff_pow_eq_one.mpr h rw [horder] at hdvd exact hdvd · intro h rw [hpow] have hdvd : orderOf g ∣ (i : ℕ) * n := by rw [horder] exact h exact orderOf_dvd_iff_pow_eq_one.mp hdvd rw [← h1, h2] rw [card_fin_dvd_mul hN] simp [N, Nat.gcd_comm]

The per-prime-power count: the number of solutions to x^m = 1 in the unit group of ZMod (p^e) is gcd(m, φ(p^e)), since that group is cyclic for odd prime powers.

lemma card_pow_eq_one_prime_pow {p e m : ℕ} (hp : Nat.Prime p) (hp2 : p ≠ 2) : Nat.card {x : (ZMod (p ^ e))ˣ // x ^ m = 1} = Nat.gcd m (Nat.totient (p ^ e)) := by classical haveI : NeZero (p ^ e) := ⟨pow_ne_zero e hp.ne_zero⟩ have hcyc : IsCyclic (ZMod (p ^ e))ˣ := ZMod.isCyclic_units_of_prime_pow p hp hp2 e have hcard : Fintype.card (ZMod (p ^ e))ˣ = Nat.totient (p ^ e) := by exact ZMod.card_units_eq_totient (p ^ e) have h := CLRS.Chapter31.card_pow_eq_one_cyclic (α := (ZMod (p ^ e))ˣ) m rw [hcard] at h exact h

In a commutative finite group, the fiber of the m-th power map over any element in its image has the same size as the kernel (the m-torsion).

lemma card_pow_eq_c_of_exists {α : Type*} [CommGroup α] [Fintype α] [DecidableEq α] {m : ℕ} {c : α} (hc : ∃ x, x ^ m = c) : Nat.card {x : α // x ^ m = c} = Nat.card {x : α // x ^ m = 1} := by classical rcases hc with ⟨x₀, hx₀⟩ let e : {x : α // x ^ m = 1} ≃ {x : α // x ^ m = c} := { toFun := fun k => ⟨x₀ * k.1, by rw [mul_pow, hx₀, k.2] simp⟩ invFun := fun x => ⟨x₀⁻¹ * x.1, by rw [mul_pow, inv_pow, hx₀, x.2] simp⟩ left_inv := by intro k apply Subtype.ext simp [This simp argument is unused: mul_assoc Hint: Omit it from the simp argument list. simp ̵[̵m̵u̵l̵_̵a̵s̵s̵o̵c̵]̵ Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`mul_assoc] right_inv := by intro x apply Subtype.ext simp [This simp argument is unused: mul_assoc Hint: Omit it from the simp argument list. simp ̵[̵m̵u̵l̵_̵a̵s̵s̵o̵c̵]̵ Note: This linter can be disabled with `set_option linter.unusedSimpArgs false`mul_assoc] } exact (Nat.card_congr e).symm

The fiber of the m-th power map over any element is no larger than the kernel (it is empty, or a coset of the kernel).

lemma card_pow_le_card_pow_eq_one {α : Type*} [CommGroup α] [Fintype α] [DecidableEq α] {m : ℕ} {c : α} : Nat.card {x : α // x ^ m = c} ≤ Nat.card {x : α // x ^ m = 1} := by classical by_cases hc : ∃ x, x ^ m = c · rw [card_pow_eq_c_of_exists hc] · have : Nat.card {x : α // x ^ m = c} = 0 := by rw [Nat.card_eq_fintype_card] rw [Fintype.card_eq_zero_iff] exact ⟨fun x => hc ⟨x.1, x.2⟩⟩ rw [this] exact Nat.zero_le _

CRT decomposition of the m-torsion. The number of units x modulo n with x^m = 1 is the product over the prime factors p of n of the number of units modulo p^(e_p) (e_p = v_p(n)) with x^m = 1.

lemma card_pow_eq_one_crt {n m : ℕ} (hn : n ≠ 0) : Nat.card {x : (ZMod n)ˣ // x ^ m = 1} = ∏ p : n.primeFactors, Nat.card {x : (ZMod (p ^ (n.factorization p)))ˣ // x ^ m = 1} := by classical let e : (ZMod n)ˣ ≃* Π p : n.primeFactors, (ZMod (p ^ (n.factorization p)))ˣ := (Units.mapEquiv (ZMod.equivPi n hn : ZMod n ≃* Π p : n.primeFactors, ZMod (p ^ n.factorization p))).trans (MulEquiv.piUnits) have h1 : Nat.card {x : (ZMod n)ˣ // x ^ m = 1} = Nat.card {f : Π p : n.primeFactors, (ZMod (p ^ (n.factorization p)))ˣ // f ^ m = 1} := by apply Nat.card_congr refine (Equiv.subtypeEquiv e ?_) intro x constructor · intro hx calc (e x) ^ m = e (x ^ m) := (map_pow e x m).symm _ = 1 := by simp [hx] · intro hx apply e.injective calc e (x ^ m) = (e x) ^ m := map_pow e x m _ = 1 := hx _ = e 1 := by simp have h2 : Nat.card {f : Π p : n.primeFactors, (ZMod (p ^ (n.factorization p)))ˣ // f ^ m = 1} = ∏ p : n.primeFactors, Nat.card {b : (ZMod (p ^ (n.factorization p)))ˣ // b ^ m = 1} := by have h3 : Nat.card {f : Π p : n.primeFactors, (ZMod (p ^ (n.factorization p)))ˣ // f ^ m = 1} = Nat.card {f : Π p : n.primeFactors, (ZMod (p ^ (n.factorization p)))ˣ // ∀ p, (f p) ^ m = 1} := by apply Nat.card_congr refine (Equiv.subtypeEquiv (Equiv.refl _) ?_) intro f constructor · intro hf p have := congrFun hf p simpa using this · intro hf funext p simpa using hf p rw [h3] rw [Nat.card_congr (Equiv.subtypePiEquivPi (p := fun p y => y ^ m = 1))] exact Nat.card_pi rw [h1, h2]

An odd divisor of 2·u divides u.

lemma odd_dvd_of_dvd_mul_two {d u : ℕ} (hd2 : d ∣ 2 * u) (hdodd : Odd d) : d ∣ u := by have hcop : d.Coprime 2 := Nat.coprime_two_right.mpr hdodd exact hcop.dvd_of_dvd_mul_right (by simpa [mul_comm] using hd2)

A divisor of an odd number is odd.

lemma odd_of_dvd_odd {d t : ℕ} (ht : Odd t) (hdt : d ∣ t) : Odd d := by by_contra hd have heven : Even d := (Nat.not_odd_iff_even.mp hd) have h2d : 2 ∣ d := even_iff_two_dvd.mp heven exact ht.not_two_dvd_nat (h2d.trans hdt)

For m = 2^(ν−1)·t with t odd and 2^ν | p−1, the greatest common divisor of m and p−1 is at most (p−1)/2: its 2-adic valuation is at most ν−1 and its odd part divides the odd part (p−1)/2^ν.

lemma gcd_pow_mul_le_half {p ν t : ℕ} (hp : 0 < p - 1) (hν : 1 ≤ ν) (ht : Odd t) (hdvd : 2 ^ ν ∣ p - 1) : Nat.gcd (2 ^ (ν - 1) * t) (p - 1) ≤ (p - 1) / 2 := by rcases hdvd with ⟨u, hu⟩ have hu_pos : 0 < u := by rw [hu] at hp apply Nat.pos_of_ne_zero intro hu0 rw [hu0] at hp norm_num at hp have hu2 : 2 ^ ν * u = 2 ^ (ν - 1) * (2 * u) := by rw [← mul_assoc] congr 1 rw [← pow_succ, Nat.sub_add_cancel hν] rw [hu, hu2] have hg : (2 ^ (ν - 1) * t).gcd (2 ^ (ν - 1) * (2 * u)) = 2 ^ (ν - 1) * t.gcd (2 * u) := by change gcd (2 ^ (ν - 1) * t) (2 ^ (ν - 1) * (2 * u)) = 2 ^ (ν - 1) * gcd t (2 * u) rw [gcd_mul_left] simp rw [hg] have hgu : Nat.gcd t (2 * u) ≤ u := by have hd2u : Nat.gcd t (2 * u) ∣ 2 * u := Nat.gcd_dvd_right _ _ have hdodd : Odd (Nat.gcd t (2 * u)) := odd_of_dvd_odd (d := Nat.gcd t (2 * u)) ht (Nat.gcd_dvd_left _ _) have hdu : Nat.gcd t (2 * u) ∣ u := odd_dvd_of_dvd_mul_two (d := Nat.gcd t (2 * u)) hd2u hdodd exact Nat.le_of_dvd hu_pos hdu have hdiv : (2 ^ (ν - 1) * (2 * u)) / 2 = 2 ^ (ν - 1) * u := by rw [show 2 ^ (ν - 1) * (2 * u) = 2 * (2 ^ (ν - 1) * u) by ring] exact Nat.mul_div_right (2 ^ (ν - 1) * u) (by norm_num) rw [hdiv] exact Nat.mul_le_mul_left (2 ^ (ν - 1)) hgu

|{x : α // p x}| equals the number of elements of α satisfying p.

try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false` lemma card_subtype_filter {α : Type*} [Fintype α] (p : α → Prop) [DecidablePred p] : Nat.card {x : α // p x} = (Finset.univ.filter p).card := by rw [Nat.card_eq_fintype_card] try 'simp' instead of 'simpa' Note: This linter can be disabled with `set_option linter.unnecessarySimpa false`simpa [Fintype.card_subtype]

The number of elements satisfying p ∨ q is at most the sum of the numbers satisfying p and q separately.

lemma card_or_le {α : Type*} [Fintype α] [DecidableEq α] (p q : α → Prop) [DecidablePred p] [DecidablePred q] : Nat.card {x : α // p x ∨ q x} ≤ Nat.card {x : α // p x} + Nat.card {x : α // q x} := by classical rw [card_subtype_filter (fun x => p x ∨ q x), card_subtype_filter p, card_subtype_filter q] have hset : (Finset.univ.filter p) ∪ (Finset.univ.filter q) = Finset.univ.filter (fun x => p x ∨ q x) := by ext x simp [Finset.mem_union, Finset.mem_filter] rw [← hset] exact Finset.card_union_le _ _

The good set {x : x^m ∈ {±1}} has size at most twice the m-torsion: it splits into the x^m = 1 and x^m = −1 parts, and the latter is no larger than the former (a fiber of the power map).

lemma goodSet_card_le {n : ℕ} [NeZero n] (m : ℕ) : Nat.card {x : (ZMod n)ˣ // x ∈ goodSet m} ≤ 2 * Nat.card {x : (ZMod n)ˣ // x ^ m = 1} := by classical have h1 : Nat.card {x : (ZMod n)ˣ // x ∈ goodSet m} = Nat.card {x : (ZMod n)ˣ // x ^ m = 1 ∨ x ^ m = -1} := by apply Nat.card_congr refine (Equiv.subtypeEquiv (Equiv.refl _) ?_) intro x exact mem_goodSet_iff m x rw [h1] have h2 := card_or_le (α := (ZMod n)ˣ) (fun x => x ^ m = 1) (fun x => x ^ m = -1) nlinarith [h2, card_pow_le_card_pow_eq_one (α := (ZMod n)ˣ) (m := m) (c := -1)]

The m-torsion of (ZMod n)ˣ is the product over the prime factors p of n of gcd(m, φ(p^(e_p))), where e_p = v_p(n).

lemma mTorsion_eq_prod {n m : ℕ} (hn : n ≠ 0) (hn_odd : Odd n) : Nat.card {x : (ZMod n)ˣ // x ^ m = 1} = ∏ p : n.primeFactors, Nat.gcd m (Nat.totient (p ^ n.factorization p)) := by rw [card_pow_eq_one_crt hn] apply Finset.prod_congr rfl intro p hp have hpp : Nat.Prime (p : ℕ) := Nat.prime_of_mem_primeFactors p.2 have hpne2 : (p : ℕ) ≠ 2 := by intro h2 have hdvd : (2 : ℕ) ∣ n := by rw [← h2] exact Nat.dvd_of_mem_primeFactors p.2 have hne : ¬ Even n := (Nat.not_even_iff_odd.mpr hn_odd) exact hne (even_iff_two_dvd.mpr hdvd) exact card_pow_eq_one_prime_pow (p := (p : ℕ)) (e := n.factorization (p : ℕ)) (m := m) hpp hpne2

ν(n) ≤ v₂(n−1): every prime factor is ≡ 1 (mod 2^ν), so n ≡ 1 and 2^ν | n−1.

lemma nu_le_v2_nat_sub_one {n : ℕ} (hn1 : 1 < n) (Variable name `hn_odd` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hn_odd : Odd n) : nu n ≤ (n - 1).factorization 2 := by have hnne : n.primeFactors.Nonempty := (Nat.nonempty_primeFactors).2 hn1 have hmod : ∀ p ∈ n.primeFactors, p ≡ 1 [MOD 2 ^ nu n] := by intro p hp exact ((Nat.modEq_iff_dvd' (a := 1) (b := p) (Nat.succ_le_of_lt (Nat.pos_of_mem_primeFactors hp))).mpr (two_pow_nu_dvd_prime_sub_one (n := n) hnne hp)).symm have hprod : ∏ p : n.primeFactors, (p : ℕ) ^ n.factorization (p : ℕ) ≡ 1 [MOD 2 ^ nu n] := by exact Nat.ModEq.prod_one (s := Finset.univ) (f := fun p : n.primeFactors => (p : ℕ) ^ n.factorization (p : ℕ)) (by intro p hp simpa using (hmod p p.2).pow (n.factorization p)) have hn_mod : n ≡ 1 [MOD 2 ^ nu n] := by conv_lhs => rw [Nat.prod_pow_primeFactors_factorization (by omega : n ≠ 0)] exact hprod have hdvd : 2 ^ nu n ∣ n - 1 := (Nat.modEq_iff_dvd' (a := 1) (b := n) (by omega)).mp hn_mod.symm exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by omega)).1 hdvd

m = 2^(ν−1)·t divides n−1 (since t·2^s = n−1 and ν ≤ s).

lemma mExp_dvd {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) : 2 ^ (nu n - 1) * (strongTestParams n).2 ∣ n - 1 := by have hs : (strongTestParams n).2 * 2 ^ (strongTestParams n).1 = n - 1 := by unfold strongTestParams have hdvd : 2 ^ (n - 1).factorization 2 ∣ n - 1 := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by omega)).2 le_rfl rw [Nat.mul_comm] exact Nat.mul_div_cancel' hdvd have hν_le_s : nu n - 1 ≤ (strongTestParams n).1 := by have hν_s : nu n ≤ (n - 1).factorization 2 := nu_le_v2_nat_sub_one hn1 hn_odd unfold strongTestParams omega have hpow_dvd : 2 ^ (nu n - 1) ∣ 2 ^ (strongTestParams n).1 := by exact pow_dvd_pow 2 hν_le_s have hmul : 2 ^ (nu n - 1) * (strongTestParams n).2 ∣ 2 ^ (strongTestParams n).1 * (strongTestParams n).2 := by exact Nat.mul_dvd_mul hpow_dvd (dvd_refl _) rw [← hs] simpa [mul_comm] using hmul

m is coprime to every prime factor p of n: m | n−1 and p | n.

lemma mExp_coprime_prime {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) {p : ℕ} (hp : p ∈ n.primeFactors) : (2 ^ (nu n - 1) * (strongTestParams n).2).Coprime p := by have hmn : 2 ^ (nu n - 1) * (strongTestParams n).2 ∣ n - 1 := mExp_dvd (n := n) hn1 hn_odd have hpn : p ∣ n := Nat.dvd_of_mem_primeFactors hp have hcop : (n - 1).Coprime n := by exact ((Nat.coprime_self_sub_right (m := 1) (n := n) (by omega)).mpr (by simp)).symm exact (hcop.of_dvd_left hmn).of_dvd_right hpn

gcd a (b·c) = gcd a c when a is coprime to b.

lemma gcd_eq_gcd_of_coprime {a b c : ℕ} (h : a.Coprime b) : a.gcd (b * c) = a.gcd c := by apply Nat.dvd_antisymm · apply Nat.dvd_gcd · exact Nat.gcd_dvd_left _ _ · have hd2 : a.gcd (b * c) ∣ b * c := Nat.gcd_dvd_right _ _ have hcop : (a.gcd (b * c)).Coprime b := h.of_dvd_left (Nat.gcd_dvd_left _ _) exact hcop.dvd_of_dvd_mul_right (by simpa [mul_comm] using hd2) · apply Nat.dvd_gcd · exact Nat.gcd_dvd_left _ _ · exact (Nat.gcd_dvd_right _ _).trans (dvd_mul_left c b)

gcd(m, φ(p^e)) = gcd(m, p−1) when gcd(m, p) = 1 and 0 < e.

lemma gcd_totient_eq_gcd_prime {p e m : ℕ} (hp : Nat.Prime p) (he : 0 < e) (hpm : m.Coprime p) : m.gcd (Nat.totient (p ^ e)) = m.gcd (p - 1) := by rw [Nat.totient_prime_pow hp he] have hpm' : m.Coprime (p ^ (e - 1)) := hpm.pow_right (e - 1) exact gcd_eq_gcd_of_coprime (a := m) (b := p ^ (e - 1)) (c := p - 1) hpm'

The m-torsion is at most the product over the prime factors of (p−1)/2, where m = 2^(ν(n)−1)·t.

lemma mTorsion_le_prod_half {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) : Nat.card {x : (ZMod n)ˣ // x ^ (2 ^ (nu n - 1) * (strongTestParams n).2) = 1} ≤ ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by rw [mTorsion_eq_prod (by omega) hn_odd] have hν : 1 ≤ nu n := nu_pos hn_odd hn1 have ht : Odd (strongTestParams n).2 := strongTestParams_odd hn1 have hnne : n.primeFactors.Nonempty := (Nat.nonempty_primeFactors).2 hn1 apply Finset.prod_le_prod' intro p hp have hpp : Nat.Prime (p : ℕ) := Nat.prime_of_mem_primeFactors p.2 have hppos : 0 < (p : ℕ) - 1 := by have := hpp.two_le omega have hpm : (2 ^ (nu n - 1) * (strongTestParams n).2).Coprime (p : ℕ) := mExp_coprime_prime (n := n) hn1 hn_odd p.2 have hdvd : 2 ^ nu n ∣ (p : ℕ) - 1 := two_pow_nu_dvd_prime_sub_one (n := n) hnne p.2 have he : 0 < n.factorization (p : ℕ) := by exact hpp.factorization_pos_of_dvd (by omega) (Nat.dvd_of_mem_primeFactors p.2) calc (2 ^ (nu n - 1) * (strongTestParams n).2).gcd (Nat.totient ((p : ℕ) ^ n.factorization (p : ℕ))) = (2 ^ (nu n - 1) * (strongTestParams n).2).gcd ((p : ℕ) - 1) := by exact gcd_totient_eq_gcd_prime hpp he hpm _ ≤ ((p : ℕ) - 1) / 2 := gcd_pow_mul_le_half (p := (p : ℕ)) (ν := nu n) (t := (strongTestParams n).2) hppos hν ht hdvd

For an odd prime p, 2 ∣ p − 1.

lemma two_dvd_prime_sub_one_of_odd {p : ℕ} (hp : Nat.Prime p) (hp_odd : Odd p) : 2 ∣ p - 1 := by rcases (odd_iff_exists_bit1.mp hp_odd) with ⟨k, rfl⟩ exact ⟨k, by omega⟩

A product over the single prime factor p of n is just the value at p.

lemma prod_primeFactors_singleton {n p : ℕ} (hpf : n.primeFactors = ({p} : Finset ℕ)) (f : ℕ → ℕ) : ∏ q : n.primeFactors, f (q : ℕ) = f p := by rw [hpf] simp

The product over the prime factors p | n of p−1 equals 2^k times the product of (p−1)/2, since each odd prime factor satisfies 2 | p−1.

lemma prod_prime_sub_one_eq_two_mul {n : ℕ} (hn_odd : Odd n) (Variable name `hn1` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hn1 : 1 < n) : ∏ p : n.primeFactors, ((p : ℕ) - 1) = 2 ^ n.primeFactors.card * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by calc ∏ p : n.primeFactors, ((p : ℕ) - 1) = ∏ p : n.primeFactors, 2 * (((p : ℕ) - 1) / 2) := by apply Finset.prod_congr rfl intro p hp have hpp : Nat.Prime (p : ℕ) := Nat.prime_of_mem_primeFactors p.2 have hpdvd : (p : ℕ) ∣ n := Nat.dvd_of_mem_primeFactors p.2 have hpodd : Odd (p : ℕ) := odd_of_dvd_odd hn_odd hpdvd have h2 : 2 ∣ (p : ℕ) - 1 := two_dvd_prime_sub_one_of_odd hpp hpodd rw [mul_comm] exact (Nat.div_mul_cancel h2).symm _ = (∏ p : n.primeFactors, 2) * ∏ p : n.primeFactors, (((p : ℕ) - 1) / 2) := by rw [Finset.prod_mul_distrib] _ = 2 ^ n.primeFactors.card * ∏ p : n.primeFactors, (((p : ℕ) - 1) / 2) := by simp [Finset.prod_const]

∏_{p|n}(p−1) ≤ n−1: each p−1 < p, and ∏ p ≤ n.

lemma prod_prime_sub_one_le {n : ℕ} (hn1 : 1 < n) : ∏ p : n.primeFactors, ((p : ℕ) - 1) ≤ n - 1 := by have hlt : ∏ p : n.primeFactors, ((p : ℕ) - 1) < ∏ p : n.primeFactors, (p : ℕ) := by refine Finset.prod_lt_prod (s := (Finset.univ : Finset n.primeFactors)) ?_ ?_ ?_ · intro p hp have hpp : Nat.Prime (p : ℕ) := Nat.prime_of_mem_primeFactors p.2 have h2 : 2 ≤ (p : ℕ) := hpp.two_le omega · intro p hp exact Nat.sub_le _ _ · rcases (Nat.nonempty_primeFactors).2 hn1 with ⟨p0, hp0⟩ refine ⟨⟨p0, hp0⟩, by simp, ?_⟩ have hpp : Nat.Prime p0 := Nat.prime_of_mem_primeFactors hp0 have h2 : 2 ≤ p0 := hpp.two_le change p0 - 1 < p0 omega have hle : ∏ p : n.primeFactors, (p : ℕ) ≤ n := by calc ∏ p : n.primeFactors, (p : ℕ) ≤ ∏ p : n.primeFactors, (p : ℕ) ^ (n.factorization (p : ℕ)) := by refine Finset.prod_le_prod' (s := (Finset.univ : Finset n.primeFactors)) (f := fun p : n.primeFactors => (p : ℕ)) (g := fun p : n.primeFactors => (p : ℕ) ^ (n.factorization (p : ℕ))) ?_ intro p hp have hpp : Nat.Prime (p : ℕ) := Nat.prime_of_mem_primeFactors p.2 have hpos : 0 < n.factorization (p : ℕ) := hpp.factorization_pos_of_dvd (by omega) (Nat.dvd_of_mem_primeFactors p.2) exact le_self_pow (by have h2 := hpp.two_le; omega) hpos.ne' _ = n := (Nat.prod_pow_primeFactors_factorization (by omega : n ≠ 0)).symm omega

The good subgroup S(n) has size at most 2·∏_{p|n}(p−1)/2: the union bound goodSet_card_le splits off the factor 2, and the m-torsion is at most ∏(p−1)/2 by mTorsion_le_prod_half.

lemma goodUnits_card_le_prodHalf {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 2 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by unfold goodUnits calc Nat.card {x : (ZMod n)ˣ // x ∈ goodSet (2 ^ (nu n - 1) * (strongTestParams n).2)} ≤ 2 * Nat.card {x : (ZMod n)ˣ // x ^ (2 ^ (nu n - 1) * (strongTestParams n).2) = 1} := goodSet_card_le (n := n) (2 ^ (nu n - 1) * (strongTestParams n).2) _ ≤ 2 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by exact Nat.mul_le_mul_left 2 (mTorsion_le_prod_half (n := n) hn1 hn_odd)

For n with at least three prime factors, |S(n)| ≤ (n−1)/4.

lemma goodUnits_card_le_of_ge_three {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hk : 3 ≤ n.primeFactors.card) : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ (n - 1) / 4 := by have h8 : 8 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 ≤ n - 1 := by calc 8 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 ≤ 2 ^ n.primeFactors.card * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by exact Nat.mul_le_mul_right _ (by have hpow : 2 ^ 3 ≤ 2 ^ n.primeFactors.card := pow_le_pow_right₀ (by norm_num) hk norm_num at hpow exact hpow) _ = ∏ p : n.primeFactors, ((p : ℕ) - 1) := (prod_prime_sub_one_eq_two_mul (n := n) hn_odd hn1).symm _ ≤ n - 1 := prod_prime_sub_one_le (n := n) hn1 have h4 : 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ n - 1 := by calc 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 4 * (2 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2) := Nat.mul_le_mul_left 4 (goodUnits_card_le_prodHalf (n := n) hn1 hn_odd) _ = 8 * ∏ p : n.primeFactors, ((p : ℕ) - 1) / 2 := by ring _ ≤ n - 1 := h8 exact (Nat.le_div_iff_mul_le (by norm_num : 0 < 4)).mpr (by simpa [mul_comm] using h4)

For n = p^e a prime power (composite), |S(n)| ≤ (n−1)/4.

lemma goodUnits_card_le_prime_power {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hn_comp : ¬ Nat.Prime n) (hk : n.primeFactors.card = 1) : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ (n - 1) / 4 := by rcases Finset.card_eq_one.mp hk with ⟨p, hpf⟩ have hp_mem : p ∈ n.primeFactors := by rw [hpf]; simp have hpp : Nat.Prime p := Nat.prime_of_mem_primeFactors hp_mem have hpodd : Odd p := odd_of_dvd_odd hn_odd (Nat.dvd_of_mem_primeFactors hp_mem) have hne : n = p ^ (n.factorization p) := by conv_lhs => rw [Nat.prod_pow_primeFactors_factorization (by omega : n ≠ 0)] exact prod_primeFactors_singleton (n := n) (p := p) hpf (fun x : ℕ => x ^ (n.factorization x)) have he1 : 1 ≤ n.factorization p := hpp.factorization_pos_of_dvd (by omega) (Nat.dvd_of_mem_primeFactors hp_mem) have hne1 : n.factorization p ≠ 1 := by intro h1 have : n = p := by simpa [h1] using hne exact hn_comp (by simpa [this] using hpp) have he : 2 ≤ n.factorization p := by omega have hp_ne2 : p ≠ 2 := by intro h2 rw [h2] at hpodd norm_num at hpodd have hp2_le_n : p ^ 2 ≤ n := by have hp2d : p ^ 2 ∣ n := (hpp.pow_dvd_iff_le_factorization (by omega : n ≠ 0)).2 he exact Nat.le_of_dvd (by omega : 0 < n) hp2d have hsp : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ p - 1 := by calc Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 2 * ∏ q : n.primeFactors, ((q : ℕ) - 1) / 2 := goodUnits_card_le_prodHalf (n := n) hn1 hn_odd _ = 2 * ((p - 1) / 2) := by rw [hpf] simp _ = p - 1 := by have h2 : 2 ∣ p - 1 := two_dvd_prime_sub_one_of_odd hpp hpodd calc 2 * ((p - 1) / 2) = ((p - 1) / 2) * 2 := by omega _ = p - 1 := Nat.div_mul_cancel h2 have hp3 : 3 ≤ p := by exact Nat.succ_le_of_lt (lt_of_le_of_ne hpp.two_le (Ne.symm hp_ne2)) have hsq : 4 * (p - 1) ≤ p ^ 2 - 1 := by have hsqf : p ^ 2 - 1 = (p - 1) * (p + 1) := by simpa [mul_comm] using (Nat.sq_sub_sq p 1) rw [hsqf] rw [mul_comm 4 (p - 1)] exact Nat.mul_le_mul_left (p - 1) (by omega : 4 ≤ p + 1) have h4p : 4 * (p - 1) ≤ n - 1 := by exact hsq.trans (by omega) have h4 : 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ n - 1 := by calc 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 4 * (p - 1) := Nat.mul_le_mul_left 4 hsp _ ≤ n - 1 := h4p exact (Nat.le_div_iff_mul_le (by norm_num : 0 < 4)).mpr (by simpa [mul_comm] using h4)

An odd divisor g of an odd d is at most d/3.

lemma odd_divisor_le_div_three {d g : ℕ} (hd : Odd d) (hg : g ∣ d) (hlt : g ≠ d) : g ≤ d / 3 := by rcases hg with ⟨c, rfl⟩ have hc_ne1 : c ≠ 1 := by intro h1 apply hlt rw [h1, mul_one] have hc_ne0 : c ≠ 0 := by intro hc0 rw [hc0, mul_zero] at hd norm_num at hd have hodd_c : Odd c := (Nat.odd_mul.mp hd).2 have hc3 : 3 ≤ c := by rcases (odd_iff_exists_bit1.mp hodd_c) with ⟨k, rfl⟩ have hk0 : 1 ≤ k := by by_contra hk have hk0' : k = 0 := by omega rw [hk0'] at hc_ne1 norm_num at hc_ne1 omega rw [Nat.le_div_iff_mul_le (by norm_num : 0 < 3)] exact Nat.mul_le_mul_left g hc3

For a prime p, p−1 = 2^s·((p−1)/2^s) where s = v₂(p−1).

lemma prime_sub_one_decomp {p : ℕ} (hp : Nat.Prime p) : p - 1 = 2 ^ (p - 1).factorization 2 * ((p - 1) / 2 ^ (p - 1).factorization 2) := by have hppos : p - 1 ≠ 0 := by have h2 := hp.two_le omega have hdvd : 2 ^ (p - 1).factorization 2 ∣ p - 1 := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) hppos).2 le_rfl exact (Nat.mul_div_cancel' hdvd).symm

The odd part t of n−1 divides n−1.

lemma strongTestParams_snd_dvd {n : ℕ} (hn1 : 1 < n) : (strongTestParams n).2 ∣ n - 1 := by unfold strongTestParams change (n - 1) / 2 ^ (n - 1).factorization 2 ∣ n - 1 have hdvd : 2 ^ (n - 1).factorization 2 ∣ n - 1 := by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by omega)).2 le_rfl refine ⟨2 ^ (n - 1).factorization 2, by rw [mul_comm]; exact (Nat.mul_div_cancel' hdvd).symm⟩

ν(p·q) = min (v₂(p−1)) (v₂(q−1)) for distinct primes p, q.

lemma nu_semiprime {p q : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) : nu (p * q) = min ((p - 1).factorization 2) ((q - 1).factorization 2) := by let s := (p - 1).factorization 2 let r := (q - 1).factorization 2 have hpq_fac : (p * q).primeFactors = ({p, q} : Finset ℕ) := by rw [Nat.primeFactors_mul hp.ne_zero hq.ne_zero] rw [Nat.Prime.primeFactors hp, Nat.Prime.primeFactors hq] ext x simp unfold nu rw [dif_pos (by rw [hpq_fac]; simp)] simp [hpq_fac]

gcd(2^(ν−1)·t, 2^s·d) = 2^(ν−1)·gcd(t, d) when 1 ≤ ν ≤ s and t is odd.

lemma gcd_pow_mul_oddPart {t s ν d : ℕ} (hν1 : 1 ≤ ν) (hνs : ν ≤ s) (ht : Odd t) : (2 ^ (ν - 1) * t).gcd (2 ^ s * d) = 2 ^ (ν - 1) * t.gcd d := by have hpow : 2 ^ s = 2 ^ (ν - 1) * 2 ^ (s - ν + 1) := by rw [← pow_add] congr 1 omega rw [hpow] rw [mul_assoc] rw [Nat.gcd_mul_left] have hcop : t.Coprime (2 ^ (s - ν + 1)) := by exact (Nat.coprime_two_right.mpr ht).pow_right (s - ν + 1) rw [gcd_eq_gcd_of_coprime hcop]

(2^s·d_p)·(2^r·d_q) = 2^(s+r)·(d_p·d_q).

lemma pow_mul_mul {s r d_p d_q : ℕ} : (2 ^ s * d_p) * (2 ^ r * d_q) = 2 ^ (s + r) * (d_p * d_q) := by rw [pow_add] ring

For primes p, q ≥ 3, (p−1)(q−1) ≤ p·q−1.

lemma prod_sub_one_le {p q : ℕ} (Variable name `hp` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hp : 3 ≤ p) (hq : 3 ≤ q) : (p - 1) * (q - 1) ≤ p * q - 1 := by calc (p - 1) * (q - 1) ≤ (p - 1) * q := Nat.mul_le_mul_left (p - 1) (Nat.sub_le _ _) _ = p * q - q := by rw [Nat.mul_sub_right_distrib] simp _ ≤ p * q - 1 := by omega

A product over n.primeFactors = {p, q} is f p·f q.

lemma prod_primeFactors_pair {n p q : ℕ} (hpf : n.primeFactors = ({p, q} : Finset ℕ)) (hpq : p ≠ q) (f : ℕ → ℕ) : ∏ x : n.primeFactors, f (x : ℕ) = f p * f q := by rw [hpf] change (({p, q} : Finset ℕ).attach.prod (fun x : {x // x ∈ ({p, q} : Finset ℕ)} => f (x : ℕ))) = f p * f q simp only [Finset.prod_attach] rw [show ({p, q} : Finset ℕ) = insert p {q} by ext x; simp] rw [Finset.prod_insert] · simp · simp [hpq]

Key lemma. For n = p·q, if d_p | t and d_q | t then d_p = d_q.

lemma semiprime_key_lemma {p q : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) (Variable name `hpq` is not explicitly referenced. The binding can be removed (if unused) or named `_` (if used implicitly). Note: This linter can be disabled with `set_option linter.unusedVariables false`hpq : p ≠ q) : ((p - 1) / 2 ^ (p - 1).factorization 2) ∣ (strongTestParams (p * q)).2 → ((q - 1) / 2 ^ (q - 1).factorization 2) ∣ (strongTestParams (p * q)).2 → (p - 1) / 2 ^ (p - 1).factorization 2 = (q - 1) / 2 ^ (q - 1).factorization 2 := by classical intro hdp_t hdq_t let s := (p - 1).factorization 2 let r := (q - 1).factorization 2 let d_p : ℕ := (p - 1) / 2 ^ s let d_q : ℕ := (q - 1) / 2 ^ r let t : ℕ := (strongTestParams (p * q)).2 have hp2 : 2 ≤ p := hp.two_le have hq2 : 2 ≤ q := hq.two_le have hpq1 : 1 < p * q := by have h4 : 4 ≤ p * q := Nat.mul_le_mul hp2 hq2 omega have hdp : d_p ∣ t := by simpa [d_p, s, t] using hdp_t have hdq : d_q ∣ t := by simpa [d_q, r, t] using hdq_t have htpq : t ∣ p * q - 1 := by simpa [t] using strongTestParams_snd_dvd (n := p * q) hpq1 have hdppq : d_p ∣ p * q - 1 := hdp.trans htpq have hp_eq : p = 2 ^ s * d_p + 1 := by have h1 : p - 1 = 2 ^ s * d_p := by dsimp [d_p] simpa [s] using (prime_sub_one_decomp hp) omega have hq_eq : q = 2 ^ r * d_q + 1 := by have h1 : q - 1 = 2 ^ r * d_q := by dsimp [d_q] simpa [r] using (prime_sub_one_decomp hq) omega have hfac : p * q - 1 = d_p * (2 ^ (s + r) * d_q + 2 ^ s) + 2 ^ r * d_q := by rw [hp_eq, hq_eq] calc (2 ^ s * d_p + 1) * (2 ^ r * d_q + 1) - 1 = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q + 1 - 1 := by have h : (2 ^ s * d_p + 1) * (2 ^ r * d_q + 1) = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q + 1 := by ring rw [h] _ = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q := by rw [Nat.add_sub_cancel] _ = 2 ^ (s + r) * (d_p * d_q) + 2 ^ s * d_p + 2 ^ r * d_q := by rw [pow_mul_mul (s := s) (r := r)] _ = d_p * (2 ^ (s + r) * d_q + 2 ^ s) + 2 ^ r * d_q := by ring have hdvd_rest : d_p ∣ 2 ^ r * d_q := by rw [hfac] at hdppq have hdpk : d_p ∣ d_p * (2 ^ (s + r) * d_q + 2 ^ s) := dvd_mul_right _ _ exact (Nat.dvd_add_iff_left hdpk).mpr (by simpa [add_comm] using hdppq) have hdvd_dq : d_p ∣ d_q := by have hodd_dp : Odd d_p := by dsimp [d_p] have hp1 : p - 1 ≠ 0 := by omega simpa [s] using oddPart_odd (p - 1) hp1 have hcop : d_p.Coprime (2 ^ r) := (Nat.coprime_two_right.mpr hodd_dp).pow_right r exact hcop.dvd_of_dvd_mul_right (by simpa [mul_comm] using hdvd_rest) have hdvd_dp : d_q ∣ d_p := by have hdqq : d_q ∣ p * q - 1 := hdq.trans htpq have hfac2 : p * q - 1 = d_q * (2 ^ (s + r) * d_p + 2 ^ r) + 2 ^ s * d_p := by rw [hp_eq, hq_eq] calc (2 ^ s * d_p + 1) * (2 ^ r * d_q + 1) - 1 = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q + 1 - 1 := by have h : (2 ^ s * d_p + 1) * (2 ^ r * d_q + 1) = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q + 1 := by ring rw [h] _ = (2 ^ s * d_p) * (2 ^ r * d_q) + 2 ^ s * d_p + 2 ^ r * d_q := by rw [Nat.add_sub_cancel] _ = 2 ^ (s + r) * (d_p * d_q) + 2 ^ s * d_p + 2 ^ r * d_q := by rw [pow_mul_mul (s := s) (r := r)] _ = d_q * (2 ^ (s + r) * d_p + 2 ^ r) + 2 ^ s * d_p := by ring have hdvd_rest2 : d_q ∣ 2 ^ s * d_p := by rw [hfac2] at hdqq have hdqk : d_q ∣ d_q * (2 ^ (s + r) * d_p + 2 ^ r) := dvd_mul_right _ _ exact (Nat.dvd_add_iff_left hdqk).mpr (by simpa [add_comm] using hdqq) have hodd_dq : Odd d_q := by dsimp [d_q] have hq1 : q - 1 ≠ 0 := by omega simpa [r] using oddPart_odd (q - 1) hq1 have hcop : d_q.Coprime (2 ^ s) := (Nat.coprime_two_right.mpr hodd_dq).pow_right s exact hcop.dvd_of_dvd_mul_right (by simpa [mul_comm] using hdvd_rest2) exact Nat.dvd_antisymm hdvd_dq hdvd_dp

8·gcd(m,p−1)·gcd(m,q−1) ≤ n−1 when v₂(p−1) ≤ v₂(q−1).

lemma semiprime_gcd_bound_sle {p q : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hp2 : p ≠ 2) (hq2 : q ≠ 2) (hle : (p - 1).factorization 2 ≤ (q - 1).factorization 2) : 8 * ((2 ^ (nu (p * q) - 1) * (strongTestParams (p * q)).2).gcd (p - 1)) * ((2 ^ (nu (p * q) - 1) * (strongTestParams (p * q)).2).gcd (q - 1)) ≤ p * q - 1 := by classical let s := (p - 1).factorization 2 let r := (q - 1).factorization 2 let ν := nu (p * q) let t := (strongTestParams (p * q)).2 let d_p := (p - 1) / 2 ^ s let d_q := (q - 1) / 2 ^ r let g_p := t.gcd d_p let g_q := t.gcd d_q have hν : ν = s := by have hmin : nu (p * q) = min s r := by simpa [s, r] using nu_semiprime hp hq omega have hp3 : 3 ≤ p := by have h2 := hp.two_le omega have hq3 : 3 ≤ q := by have h2 := hq.two_le omega have hpodd : Odd p := (hp.odd_iff).mpr hp3 have hqodd : Odd q := (hq.odd_iff).mpr hq3 have ht : Odd t := by dsimp [t] exact strongTestParams_odd (n := p * q) (by have h4 : 4 ≤ p * q := Nat.mul_le_mul hp.two_le hq.two_le omega) have hν1 : 1 ≤ ν := by have hs1 : 1 ≤ s := by dsimp [s] exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by have h2 := hp.two_le; omega)).1 (by simpa using two_dvd_prime_sub_one_of_odd hp hpodd) omega have hodd_dp : Odd d_p := by dsimp [d_p] have hpp1 : p - 1 ≠ 0 := by have h2 := hp.two_le omega simpa [s] using oddPart_odd (p - 1) hpp1 have hodd_dq : Odd d_q := by dsimp [d_q] have hqq1 : q - 1 ≠ 0 := by have h2 := hq.two_le omega simpa [r] using oddPart_odd (q - 1) hqq1 have hdp : p - 1 = 2 ^ s * d_p := by dsimp [d_p] simpa [s] using prime_sub_one_decomp hp have hdq : q - 1 = 2 ^ r * d_q := by dsimp [d_q] simpa [r] using prime_sub_one_decomp hq have hgs : ν ≤ s := by omega have hgr : ν ≤ r := by omega have hga : (2 ^ (ν - 1) * t).gcd (p - 1) = 2 ^ (ν - 1) * g_p := by rw [hdp] dsimp [g_p] exact gcd_pow_mul_oddPart (t := t) (s := s) (ν := ν) (d := d_p) hν1 hgs ht have hgb : (2 ^ (ν - 1) * t).gcd (q - 1) = 2 ^ (ν - 1) * g_q := by rw [hdq] dsimp [g_q] exact gcd_pow_mul_oddPart (t := t) (s := r) (ν := ν) (d := d_q) hν1 hgr ht have hmain : 8 * (2 ^ (ν - 1) * g_p) * (2 ^ (ν - 1) * g_q) ≤ (p - 1) * (q - 1) := by have hpow8 : 8 * (2 ^ (ν - 1) * g_p) * (2 ^ (ν - 1) * g_q) = 2 ^ (2 * ν + 1) * g_p * g_q := by calc 8 * (2 ^ (ν - 1) * g_p) * (2 ^ (ν - 1) * g_q) = 8 * 2 ^ (ν - 1) * g_p * 2 ^ (ν - 1) * g_q := by ring _ = 2 ^ 3 * 2 ^ (ν - 1) * 2 ^ (ν - 1) * g_p * g_q := by rw [show 8 = 2 ^ 3 by norm_num] ring _ = 2 ^ (3 + (ν - 1) + (ν - 1)) * g_p * g_q := by rw [← pow_add] rw [← pow_add] _ = 2 ^ (2 * ν + 1) * g_p * g_q := by have hexp : 3 + (ν - 1) + (ν - 1) = 2 * ν + 1 := by omega rw [hexp] rw [hpow8] have htarget : 2 ^ (2 * ν + 1) * g_p * g_q ≤ 2 ^ (s + r) * (d_p * d_q) := by rw [hν] by_cases hsr : s = r · have h2g : 2 * g_p * g_q ≤ d_p * d_q := by have hnot : ¬ (g_p = d_p ∧ g_q = d_q) := by intro hboth rcases hboth with ⟨hg1, hg2⟩ have hd1 : d_p ∣ t := by rw [← hg1] dsimp [g_p] exact Nat.gcd_dvd_left t d_p have hd2 : d_q ∣ t := by rw [← hg2] dsimp [g_q] exact Nat.gcd_dvd_left t d_q have hdeq : d_p = d_q := semiprime_key_lemma hp hq hpq hd1 hd2 have hp_eq : p = q := by have h1 : p - 1 = 2 ^ s * d_p := by dsimp [d_p] simpa [s] using prime_sub_one_decomp hp have h2 : q - 1 = 2 ^ r * d_q := by dsimp [d_q] simpa [r] using prime_sub_one_decomp hq have hsub : p - 1 = q - 1 := by calc p - 1 = 2 ^ s * d_p := h1 _ = 2 ^ r * d_q := by rw [hsr, hdeq] _ = q - 1 := h2.symm omega exact hpq hp_eq have hprop : g_p ≠ d_p ∨ g_q ≠ d_q := by by_contra hc push Not at hc exact hnot hc have hg_p_dvd : g_p ∣ d_p := by dsimp [g_p]; exact Nat.gcd_dvd_right t d_p have hg_q_dvd : g_q ∣ d_q := by dsimp [g_q]; exact Nat.gcd_dvd_right t d_q have hd_p_pos : 0 < d_p := by dsimp [d_p] have hpp : 0 < p - 1 := by have h2 := hp.two_le omega exact Nat.div_pos (Nat.le_of_dvd hpp (by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by have h2 := hp.two_le; omega)).2 le_rfl)) (by norm_num : 0 < 2 ^ s) have hd_q_pos : 0 < d_q := by dsimp [d_q] have hqq : 0 < q - 1 := by have h2 := hq.two_le omega exact Nat.div_pos (Nat.le_of_dvd hqq (by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by have h2 := hq.two_le; omega)).2 le_rfl)) (by norm_num : 0 < 2 ^ r) rcases hprop with hne1 | hne2 · have hlt1 : g_p < d_p := lt_of_le_of_ne (Nat.le_of_dvd hd_p_pos hg_p_dvd) hne1 have h13 : g_p ≤ d_p / 3 := odd_divisor_le_div_three hodd_dp hg_p_dvd hne1 have hq_le : g_q ≤ d_q := Nat.le_of_dvd hd_q_pos hg_q_dvd have h23 : 2 * (d_p / 3) * d_q ≤ d_p * d_q := by have h2 : 2 * (d_p / 3) ≤ d_p := by omega exact Nat.mul_le_mul_right d_q h2 calc 2 * g_p * g_q ≤ 2 * (d_p / 3) * d_q := by exact Nat.mul_le_mul (Nat.mul_le_mul_left 2 h13) hq_le _ ≤ d_p * d_q := h23 · have h23 : g_q ≤ d_q / 3 := odd_divisor_le_div_three hodd_dq hg_q_dvd hne2 have hp_le : g_p ≤ d_p := Nat.le_of_dvd hd_p_pos hg_p_dvd have h23' : 2 * (d_q / 3) * d_p ≤ d_q * d_p := by have h2 : 2 * (d_q / 3) ≤ d_q := by omega exact Nat.mul_le_mul_right d_p h2 calc 2 * g_p * g_q = 2 * g_q * g_p := by ring _ ≤ 2 * (d_q / 3) * d_p := by exact Nat.mul_le_mul (Nat.mul_le_mul_left 2 h23) hp_le _ ≤ d_q * d_p := h23' _ = d_p * d_q := by ring calc 2 ^ (2 * s + 1) * g_p * g_q = 2 ^ (s + r) * (2 * g_p * g_q) := by rw [hsr] rw [pow_add] norm_num ring _ ≤ 2 ^ (s + r) * (d_p * d_q) := by exact Nat.mul_le_mul_left (2 ^ (s + r)) h2g · have hlt : s < r := lt_of_le_of_ne hle hsr have hpow_le : 2 ^ (2 * s + 1) ≤ 2 ^ (s + r) := by apply pow_le_pow_right₀ (by norm_num) omega have hg_p_le : g_p ≤ d_p := by have hpp : 0 < d_p := by dsimp [d_p] have hpp' : 0 < p - 1 := by have h2 := hp.two_le omega exact Nat.div_pos (Nat.le_of_dvd hpp' (by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by have h2 := hp.two_le; omega)).2 le_rfl)) (by norm_num : 0 < 2 ^ s) exact Nat.le_of_dvd hpp (by dsimp [g_p]; exact Nat.gcd_dvd_right t d_p) have hg_q_le : g_q ≤ d_q := by have hqq : 0 < d_q := by dsimp [d_q] have hqq' : 0 < q - 1 := by have h2 := hq.two_le omega exact Nat.div_pos (Nat.le_of_dvd hqq' (by exact (Nat.Prime.pow_dvd_iff_le_factorization (by decide : Nat.Prime 2) (by have h2 := hq.two_le; omega)).2 le_rfl)) (by norm_num : 0 < 2 ^ r) exact Nat.le_of_dvd hqq (by dsimp [g_q]; exact Nat.gcd_dvd_right t d_q) calc 2 ^ (2 * s + 1) * g_p * g_q ≤ 2 ^ (s + r) * d_p * d_q := by exact Nat.mul_le_mul (Nat.mul_le_mul hpow_le hg_p_le) hg_q_le _ = 2 ^ (s + r) * (d_p * d_q) := by ring calc 2 ^ (2 * ν + 1) * g_p * g_q ≤ 2 ^ (s + r) * (d_p * d_q) := htarget _ = (p - 1) * (q - 1) := by rw [hdp, hdq] rw [pow_mul_mul (s := s) (r := r)] calc 8 * ((2 ^ (ν - 1) * t).gcd (p - 1)) * ((2 ^ (ν - 1) * t).gcd (q - 1)) = 8 * (2 ^ (ν - 1) * g_p) * (2 ^ (ν - 1) * g_q) := by rw [hga, hgb] _ ≤ (p - 1) * (q - 1) := hmain _ ≤ p * q - 1 := prod_sub_one_le hp3 hq3

8·gcd(m,p−1)·gcd(m,q−1) ≤ n−1 for n with n.primeFactors = {p,q}.

lemma semiprime_gcd_bound {p q : ℕ} (hp : Nat.Prime p) (hq : Nat.Prime q) (hpq : p ≠ q) (hp2 : p ≠ 2) (hq2 : q ≠ 2) : 8 * ((2 ^ (nu (p * q) - 1) * (strongTestParams (p * q)).2).gcd (p - 1)) * ((2 ^ (nu (p * q) - 1) * (strongTestParams (p * q)).2).gcd (q - 1)) ≤ p * q - 1 := by by_cases hle : (p - 1).factorization 2 ≤ (q - 1).factorization 2 · exact semiprime_gcd_bound_sle hp hq hpq hp2 hq2 hle · have hle' : (q - 1).factorization 2 ≤ (p - 1).factorization 2 := by omega have hmain := semiprime_gcd_bound_sle (p := q) (q := p) hq hp (Ne.symm hpq) hq2 hp2 hle' simpa [mul_comm, mul_left_comm, mul_assoc] using hmain

For primes p, q ≥ 3, 2·(p−1)(q−1) ≤ p²q−1.

lemma crude_bound {p q : ℕ} (hp3 : 3 ≤ p) (hq3 : 3 ≤ q) : 2 * (p - 1) * (q - 1) ≤ p ^ 2 * q - 1 := by have hpq_le : (p - 1) * (q - 1) ≤ p * q := Nat.mul_le_mul (Nat.sub_le _ _) (Nat.sub_le _ _) have h2pq : 2 * (p * q) ≤ p ^ 2 * q - 1 := by have h1pq : 1 ≤ (p - 2) * (p * q) := by have hp2 : 1 ≤ p - 2 := by omega have hpqpos : 1 ≤ p * q := by have h9 : 9 ≤ p * q := Nat.mul_le_mul hp3 hq3 omega have hmul : 1 * 1 ≤ (p - 2) * (p * q) := Nat.mul_le_mul hp2 hpqpos simpa using hmul have h'' : 1 ≤ p ^ 2 * q - 2 * (p * q) := by rw [Nat.mul_sub_right_distrib] at h1pq simpa [pow_two, Nat.mul_assoc] using h1pq have hgoal : 2 * (p * q) + 1 ≤ p ^ 2 * q := by omega omega calc 2 * (p - 1) * (q - 1) = 2 * ((p - 1) * (q - 1)) := by ring _ ≤ 2 * (p * q) := Nat.mul_le_mul_left 2 hpq_le _ ≤ p ^ 2 * q - 1 := h2pq

For primes p, q ≥ 3, 2·(p−1)(q−1) ≤ p·q²−1.

lemma crude_bound' {p q : ℕ} (hp3 : 3 ≤ p) (hq3 : 3 ≤ q) : 2 * (p - 1) * (q - 1) ≤ p * q ^ 2 - 1 := by have hpq_le : (p - 1) * (q - 1) ≤ p * q := Nat.mul_le_mul (Nat.sub_le _ _) (Nat.sub_le _ _) have h2pq : 2 * (p * q) ≤ p * q ^ 2 - 1 := by have h1pq : 1 ≤ (q - 2) * (p * q) := by have hq2 : 1 ≤ q - 2 := by omega have hpqpos : 1 ≤ p * q := by have h9 : 9 ≤ p * q := Nat.mul_le_mul hp3 hq3 omega have hmul : 1 * 1 ≤ (q - 2) * (p * q) := Nat.mul_le_mul hq2 hpqpos simpa using hmul have h'' : 1 ≤ p * q ^ 2 - 2 * (p * q) := by rw [Nat.mul_sub_right_distrib] at h1pq simpa [pow_two, Nat.mul_assoc, Nat.mul_comm, Nat.mul_left_comm] using h1pq have hgoal : 2 * (p * q) + 1 ≤ p * q ^ 2 := by omega omega calc 2 * (p - 1) * (q - 1) = 2 * ((p - 1) * (q - 1)) := by ring _ ≤ 2 * (p * q) := Nat.mul_le_mul_left 2 hpq_le _ ≤ p * q ^ 2 - 1 := h2pq

For n with exactly two distinct prime factors, |S(n)| ≤ (n−1)/4.

lemma goodUnits_card_le_semiprime {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hk : n.primeFactors.card = 2) : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ (n - 1) / 4 := by rcases Finset.card_eq_two.mp hk with ⟨p, q, hpq_ne, hpq⟩ have hp_mem : p ∈ n.primeFactors := by rw [hpq]; simp have hq_mem : q ∈ n.primeFactors := by rw [hpq]; simp have hp : Nat.Prime p := Nat.prime_of_mem_primeFactors hp_mem have hq : Nat.Prime q := Nat.prime_of_mem_primeFactors hq_mem have hp_ne2 : p ≠ 2 := by intro h2 have hpdvd : p ∣ n := Nat.dvd_of_mem_primeFactors hp_mem have hpodd : Odd p := odd_of_dvd_odd hn_odd hpdvd rw [h2] at hpodd norm_num at hpodd have hq_ne2 : q ≠ 2 := by intro h2 have hqdvd : q ∣ n := Nat.dvd_of_mem_primeFactors hq_mem have hqodd : Odd q := odd_of_dvd_odd hn_odd hqdvd rw [h2] at hqodd norm_num at hqodd have hne : n = p ^ (n.factorization p) * q ^ (n.factorization q) := by conv_lhs => rw [Nat.prod_pow_primeFactors_factorization (by omega : n ≠ 0)] exact prod_primeFactors_pair hpq hpq_ne (fun x : ℕ => x ^ (n.factorization x)) let a := n.factorization p let b := n.factorization q have ha : 1 ≤ a := hp.factorization_pos_of_dvd (by omega) (Nat.dvd_of_mem_primeFactors hp_mem) have hb : 1 ≤ b := hq.factorization_pos_of_dvd (by omega) (Nat.dvd_of_mem_primeFactors hq_mem) have hp2le : 2 ≤ p := hp.two_le have hq2le : 2 ≤ q := hq.two_le have hp3 : 3 ≤ p := by omega have hq3 : 3 ≤ q := by omega by_cases hsq : a = 1 ∧ b = 1 · have hnpq : n = p * q := by rw [hne] change p ^ a * q ^ b = p * q rw [hsq.1, hsq.2] simp let m := 2 ^ (nu n - 1) * (strongTestParams n).2 have hmTorsion : Nat.card {x : (ZMod n)ˣ // x ^ m = 1} = m.gcd (p - 1) * m.gcd (q - 1) := by rw [mTorsion_eq_prod (n := n) (by omega : n ≠ 0) hn_odd] rw [prod_primeFactors_pair hpq hpq_ne (fun x : ℕ => m.gcd (Nat.totient (x ^ (n.factorization x))))] have hpm : m.Coprime p := by simpa [m] using mExp_coprime_prime (n := n) hn1 hn_odd hp_mem have hqm : m.Coprime q := by simpa [m] using mExp_coprime_prime (n := n) hn1 hn_odd hq_mem rw [gcd_totient_eq_gcd_prime hp (by omega : 0 < a) hpm, gcd_totient_eq_gcd_prime hq (by omega : 0 < b) hqm] have hS : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 2 * m.gcd (p - 1) * m.gcd (q - 1) := by unfold goodUnits calc Nat.card {x : (ZMod n)ˣ // x ∈ goodSet m} ≤ 2 * Nat.card {x : (ZMod n)ˣ // x ^ m = 1} := goodSet_card_le (n := n) m _ = 2 * (m.gcd (p - 1) * m.gcd (q - 1)) := by rw [hmTorsion] _ = 2 * m.gcd (p - 1) * m.gcd (q - 1) := by ring have hb8 : 8 * m.gcd (p - 1) * m.gcd (q - 1) ≤ n - 1 := by have hb' := semiprime_gcd_bound (p := p) (q := q) hp hq hpq_ne hp_ne2 hq_ne2 rw [← hnpq] at hb' simpa [m] using hb' have h4 : 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ n - 1 := by calc 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 4 * (2 * m.gcd (p - 1) * m.gcd (q - 1)) := by exact Nat.mul_le_mul_left 4 hS _ = 8 * m.gcd (p - 1) * m.gcd (q - 1) := by ring _ ≤ n - 1 := hb8 exact (Nat.le_div_iff_mul_le (by norm_num : 0 < 4)).mpr (by simpa [mul_comm] using h4) · have hnsq : 2 ≤ a ∨ 2 ≤ b := by omega have hS' : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 2 * ((p - 1) / 2) * ((q - 1) / 2) := by calc Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 2 * ∏ x : n.primeFactors, ((x : ℕ) - 1) / 2 := goodUnits_card_le_prodHalf (n := n) hn1 hn_odd _ = 2 * (((p - 1) / 2) * ((q - 1) / 2)) := by rw [prod_primeFactors_pair hpq hpq_ne (fun x : ℕ => (x - 1) / 2)] _ = 2 * ((p - 1) / 2) * ((q - 1) / 2) := by ring have h8 : 8 * ((p - 1) / 2) * ((q - 1) / 2) ≤ n - 1 := by have hp_even : 2 ∣ p - 1 := two_dvd_prime_sub_one_of_odd hp (by have hpdvd : p ∣ n := Nat.dvd_of_mem_primeFactors hp_mem exact odd_of_dvd_odd hn_odd hpdvd) have hq_even : 2 ∣ q - 1 := two_dvd_prime_sub_one_of_odd hq (by have hqdvd : q ∣ n := Nat.dvd_of_mem_primeFactors hq_mem exact odd_of_dvd_odd hn_odd hqdvd) have h8eq : 8 * ((p - 1) / 2) * ((q - 1) / 2) = 2 * (p - 1) * (q - 1) := by calc 8 * ((p - 1) / 2) * ((q - 1) / 2) = 2 * (2 * ((p - 1) / 2)) * (2 * ((q - 1) / 2)) := by ring _ = 2 * (p - 1) * (q - 1) := by rw [Nat.mul_div_cancel' hp_even, Nat.mul_div_cancel' hq_even] rw [h8eq] rcases hnsq with ha2 | hb2 · have hp2_n : p ^ 2 * q ≤ n := by rw [hne] have hp2a : p ^ 2 ≤ p ^ a := pow_le_pow_right₀ (by omega) ha2 have hq1b : q ≤ q ^ b := le_self_pow (by omega) (by omega : b ≠ 0) exact Nat.mul_le_mul hp2a hq1b exact (crude_bound hp3 hq3).trans (by omega) · have hq2_n : p * q ^ 2 ≤ n := by rw [hne] have hp1a : p ≤ p ^ a := le_self_pow (by omega) (by omega : a ≠ 0) have hq2b : q ^ 2 ≤ q ^ b := pow_le_pow_right₀ (by omega) hb2 exact Nat.mul_le_mul hp1a hq2b exact (crude_bound' hp3 hq3).trans (by omega) have h4 : 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ n - 1 := by calc 4 * Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ 4 * (2 * ((p - 1) / 2) * ((q - 1) / 2)) := by exact Nat.mul_le_mul_left 4 hS' _ = 8 * ((p - 1) / 2) * ((q - 1) / 2) := by ring _ ≤ n - 1 := h8 exact (Nat.le_div_iff_mul_le (by norm_num : 0 < 4)).mpr (by simpa [mul_comm] using h4)

The good subgroup S(n) has size at most (n−1)/4 for odd composite n (Rabin–Monier). The proof splits on the number k of distinct prime factors: k = 1 (prime power), k = 2 (semiprime), and k ≥ 3.

theorem goodUnits_card_le {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hn_comp : ¬ Nat.Prime n) : Nat.card {x : (ZMod n)ˣ // x ∈ goodUnits} ≤ (n - 1) / 4 := by have hcard : 1 ≤ n.primeFactors.card := by have hpos : 0 < n.primeFactors.card := (Finset.card_pos).mpr ((Nat.nonempty_primeFactors).2 hn1) omega have hcases : n.primeFactors.card = 1 ∨ n.primeFactors.card = 2 ∨ 3 ≤ n.primeFactors.card := by omega rcases hcases with h1 | h2 | h3 · exact goodUnits_card_le_prime_power (n := n) hn1 hn_odd hn_comp h1 · exact goodUnits_card_le_semiprime (n := n) hn1 hn_odd h2 · exact goodUnits_card_le_of_ge_three (n := n) hn1 hn_odd h3

The Miller-Rabin error bound (Theorem 31.39; sharpened to (n−1)/4 by Rabin–Monier). For odd composite n, at most (n−1)/4 of the bases in (Z/nZ)ˣ are strong liars. Every strong liar lies in the good subgroup S(n) (liar_mem_goodSet), and |S(n)| ≤ (n−1)/4 (goodUnits_card_le).

theorem strongLiars_card_le {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hn_comp : ¬ Nat.Prime n) : Nat.card {a : (ZMod n)ˣ // isStrongLiar a} ≤ (n - 1) / 4 := by have hle : Nat.card {a : (ZMod n)ˣ // isStrongLiar a} ≤ Nat.card {a : (ZMod n)ˣ // a ∈ goodUnits} := by refine Nat.card_le_card_of_injective (fun a : {a : (ZMod n)ˣ // isStrongLiar a} => ⟨(a : (ZMod n)ˣ), liar_mem_goodSet (n := n) hn_odd hn1 a.2⟩) ?_ intro a b h exact Subtype.ext (by simpa using congrArg Subtype.val h) exact hle.trans (goodUnits_card_le (n := n) hn1 hn_odd hn_comp)

Random-witness analysis (the MILLER-RABIN error bound)

A strong pseudoprime base is coprime to the modulus.

theorem strongPseudoprime_coprime {n a : ℕ} (hn1 : 1 < n) (h : strongPseudoprime n a) : Nat.Coprime a n := by have hpow : a ^ (n - 1) ≡ 1 [MOD n] := strongPseudoprime_pow h have hmul : a * a ^ (n - 2) ≡ 1 [MOD n] := by have hn' : n - 1 = 1 + (n - 2) := by omega rw [hn'] at hpow rw [pow_add, pow_one] at hpow exact hpow exact Nat.coprime_of_mul_modEq_one (a ^ (n - 2)) hmul

A natural strong liar lifts to a strong liar in the unit group.

theorem isStrongLiar_of_strongPseudoprime {n : ℕ} [NeZero n] {a : ℕ} (hcop : Nat.Coprime a n) (h : strongPseudoprime n a) : isStrongLiar (ZMod.unitOfCoprime a hcop) := by rw [isStrongLiar] rw [ZMod.coe_unitOfCoprime] unfold strongPseudoprime at h rcases h with h1 | ⟨i, hi⟩ · left simpa [Nat.cast_pow] using (ZMod.natCast_eq_natCast_iff (a ^ (strongTestParams n).2) 1 n).mpr h1 · right refine ⟨i, ?_⟩ rw [← Nat.cast_pow] have hnat : ((a ^ (2 ^ (i : ℕ) * (strongTestParams n).2) : ℕ) : ZMod n) = ((n - 1 : ℕ) : ZMod n) := (ZMod.natCast_eq_natCast_iff (a ^ (2 ^ (i : ℕ) * (strongTestParams n).2)) (n - 1) n).mpr hi rw [hnat] have hneg : ((n - 1 : ℕ) : ZMod n) = -1 := by have hnpos : 0 < n := NeZero.pos n rw [Nat.cast_sub (show 1 ≤ n by omega)] rw [ZMod.natCast_self] simp exact hneg

Random-witness count bound. For odd composite n, at most (n-1)/4 of the bases 1, …, n-1 are strong liars. This is the sampling interpretation of the Miller-Rabin error bound (Theorem 31.39).

theorem strongLiars_nat_card_le {n : ℕ} [NeZero n] (hn1 : 1 < n) (hn_odd : Odd n) (hn_comp : ¬ Nat.Prime n) : Nat.card {a : Fin (n - 1) // strongPseudoprime n (a.val + 1)} ≤ (n - 1) / 4 := by let f : {a : Fin (n - 1) // strongPseudoprime n (a.val + 1)} → {u : (ZMod n)ˣ // isStrongLiar u} := fun a => let hcop : Nat.Coprime (a.val.val + 1) n := strongPseudoprime_coprime hn1 a.2 ⟨ZMod.unitOfCoprime (a.val.val + 1) hcop, isStrongLiar_of_strongPseudoprime hcop a.2⟩ have hfinj : Function.Injective f := by intro a b hab apply Subtype.ext apply Fin.ext have h : ZMod.unitOfCoprime (a.val.val + 1) (strongPseudoprime_coprime hn1 a.2) = ZMod.unitOfCoprime (b.val.val + 1) (strongPseudoprime_coprime hn1 b.2) := by change (f a).val = (f b).val exact congrArg Subtype.val hab have hcoef : ((ZMod.unitOfCoprime (a.val.val + 1) (strongPseudoprime_coprime hn1 a.2) : ZMod n)) = ((ZMod.unitOfCoprime (b.val.val + 1) (strongPseudoprime_coprime hn1 b.2) : ZMod n)) := by exact congrArg (fun x : (ZMod n)ˣ => (x : ZMod n)) h rw [ZMod.coe_unitOfCoprime, ZMod.coe_unitOfCoprime] at hcoef have hmod : a.val.val + 1 ≡ b.val.val + 1 [MOD n] := (ZMod.natCast_eq_natCast_iff (a.val.val + 1) (b.val.val + 1) n).mp hcoef have ha : a.val.val + 1 < n := by have := a.val.isLt; omega have hb : b.val.val + 1 < n := by have := b.val.isLt; omega rw [Nat.ModEq] at hmod rw [Nat.mod_eq_of_lt ha, Nat.mod_eq_of_lt hb] at hmod omega calc Nat.card {a : Fin (n - 1) // strongPseudoprime n (a.val + 1)} ≤ Nat.card {u : (ZMod n)ˣ // isStrongLiar u} := Nat.card_le_card_of_injective f hfinj _ ≤ (n - 1) / 4 := strongLiars_card_le (n := n) hn1 hn_odd hn_comp

Legacy multi-base semantic loop and detached budget. The first component is true exactly when all supplied bases pass millerRabin. The second sums the counters of independent exponentiations with exponent n - 1; those residues do not determine these Boolean decisions. This is an analysis budget, not an attached execution-cost theorem. The Execution companion provides the counted residue-based checker and a lazy multi-base loop.

def millerRabinLoop (n : ℕ) : List ℕ → Bool × ℕ | [] => (true, 0) | a :: as => let (pass, cost) := millerRabinLoop n as (millerRabin n a && pass, (modExpWithCount a n (n - 1)).2 + cost)

The loop reports "probably prime" exactly when every base passes (CLRS §31.8).

theorem millerRabinLoop_fst_iff (n : ℕ) (bases : List ℕ) : (millerRabinLoop n bases).1 = true ↔ ∀ a ∈ bases, millerRabin n a = true := by induction bases with | nil => simp [millerRabinLoop] | cons a as ih => simp [millerRabinLoop, ih, Bool.and_eq_true, List.mem_cons]

The legacy detached exponentiation budget is at most 2 · |bases| · Nat.size (n−1). This does not measure the decision procedure.

theorem millerRabinLoop_count_le (n : ℕ) (bases : List ℕ) : (millerRabinLoop n bases).2 ≤ bases.length * (2 * Nat.size (n - 1)) := by induction bases with | nil => simp [millerRabinLoop] | cons a as ih => simp [millerRabinLoop] have h := modExpWithCount_count_le a n (n - 1) calc (modExpWithCount a n (n - 1)).2 + (millerRabinLoop n as).2 ≤ 2 * Nat.size (n - 1) + as.length * (2 * Nat.size (n - 1)) := Nat.add_le_add h ih _ = (as.length + 1) * (2 * Nat.size (n - 1)) := by ring
end Chapter31end CLRS

Definitions and proofs

CLRSLean.FourthEdition.Chapter_31.Section_31_8_Primality_Testing.Execution

Miller–Rabin with counted residue execution

The single-base decision inspects the residue returned by counted modular exponentiation and the residues of subsequent counted modular squares. The multi-base loop stops at the first rejection and never evaluates its suffix. The counter measures modular multiplications, including the initial exponentiation; it is not bit complexity or total runtime. Parameter decomposition, comparisons, and supplied-base construction/sampling are outside this counter. In particular, this module does not charge the implementation of strongTestParams. The probability companion specifies independent uniform input samples separately. The kernel preserves the existing strong-probable-prime predicate on all naturals; its primality interpretation uses the existing admissible-input hypotheses.

namespace CLRS.Chapter31.MillerRabinExecution

Inspect the current residue, then at most fuel further modular squares.

def squareSearch (n target : Nat) : Nat → Nat → Bool × Nat | 0, x => (decide (x % n = target % n), 0) | k + 1, x => if x % n = target % n then (true, 0) else let next := squareSearch n target k ((x * x) % n) (next.1, next.2 + 1)
private theorem square_pow (x n i : Nat) : (((x * x) % n) ^ (2 ^ i)) % n = (x ^ (2 ^ (i + 1))) % n := by rw [Nat.pow_mod, Nat.mod_mod, ← Nat.pow_mod] rw [← pow_two, ← pow_mul] congr 2 rw [pow_succ] omegatheorem squareSearch_spec (n target fuel x : Nat) : (squareSearch n target fuel x).1 = true ↔ ∃ i ≤ fuel, x ^ (2 ^ i) % n = target % n := by induction fuel generalizing x with | zero => simp [squareSearch] | succ k ih => by_cases hx : x % n = target % n · simp only [squareSearch, hx, if_true] exact ⟨fun _ => ⟨0, by omega, by simpa using hx⟩, fun _ => trivial⟩ · simp only [squareSearch, hx, if_false, ih] constructor · rintro ⟨i, hi, he⟩ exact ⟨i + 1, by omega, by simpa only [square_pow] using he⟩ · rintro ⟨i, hi, he⟩ cases i with | zero => simp [hx] at he | succ i => exact ⟨i, by omega, by simpa only [square_pow] using he⟩theorem squareSearch_count_le (n target fuel x : Nat) : (squareSearch n target fuel x).2 ≤ fuel := by induction fuel generalizing x with | zero => simp [squareSearch] | succ k ih => simp only [squareSearch] split · simp · exact Nat.add_le_add_right (ih _) 1

Rejection exhausts exactly the allowed squarings; successful scans may stop earlier.

theorem squareSearch_reject_count (n target fuel x : Nat) (h : (squareSearch n target fuel x).1 = false) : (squareSearch n target fuel x).2 = fuel := by induction fuel generalizing x with | zero => rfl | succ k ih => by_cases hx : x % n = target % n · simp [squareSearch, hx] at h · simp only [squareSearch, hx, if_false] at h ⊢ rw [ih _ h]

The result and count come from the same exponentiation and residue scan.

def trial (n a : Nat) : Bool × Nat := let params := strongTestParams n let initial := modExpWithCount a n params.2 if initial.1 = 1 % n then (true, initial.2) else match params.1 with | 0 => (false, initial.2) | s + 1 => let scanned := squareSearch n (n - 1) s initial.1 (scanned.1, initial.2 + scanned.2)
private theorem initial_pow (a n d i : Nat) : ((a ^ d % n) ^ (2 ^ i)) % n = a ^ (2 ^ i * d) % n := by rw [Nat.pow_mod, Nat.mod_mod, ← Nat.pow_mod, ← pow_mul, Nat.mul_comm d]theorem trial_spec (n a : Nat) : (trial n a).1 = true ↔ strongPseudoprime n a := by unfold strongPseudoprime simp only [trial, modExpWithCount_spec] split · rename_i he constructor · intro _; exact Or.inl he · intro _; rfl · rename_i he cases hs : (strongTestParams n).1 with | zero => simp [Nat.ModEq, he] | succ s => simp only [squareSearch_spec, initial_pow, Nat.ModEq] constructor · rintro ⟨i, hi, h⟩ exact Or.inr ⟨⟨i, by omega⟩, h⟩ · rintro (h | ⟨i, hi⟩) · exact (he h).elim · exact ⟨i.val, by omega, hi⟩theorem trial_eq_millerRabin (n a : Nat) : (trial n a).1 = millerRabin n a := by apply Bool.eq_iff_iff.mpr simp [trial_spec, millerRabin]theorem trial_prime {n a : Nat} (hn : Nat.Prime n) (hcop : Nat.Coprime a n) : (trial n a).1 = true := (trial_spec n a).mpr (strongPseudoprime_of_prime hn hcop)theorem trial_count_le (n a : Nat) : (trial n a).2 ≤ 2 * Nat.size (strongTestParams n).2 + (strongTestParams n).1 := by have hc := modExpWithCount_count_le a n (strongTestParams n).2 simp only [trial] split · omega · cases hs : (strongTestParams n).1 with | zero => simp only; omega | succ s => have hsq := squareSearch_count_le n (n - 1) s (modExpWithCount a n (strongTestParams n).2).1 simp only omegatheorem trial_count_le_size (n a : Nat) : (trial n a).2 ≤ 3 * Nat.size (n - 1) := by have hc := trial_count_le n a have hd : Nat.size (strongTestParams n).2 ≤ Nat.size (n - 1) := Nat.size_le_size (Nat.div_le_self _ _) have hs : (strongTestParams n).1 ≤ Nat.size (n - 1) := Nat.factorization_le_of_le_pow (Nat.size_le.mp (Nat.le_refl (Nat.size (n - 1)))).le omega

Supplied bases are processed left to right, stopping before the suffix on rejection.

def run (n : Nat) : List Nat → Bool × Nat | [] => (true, 0) | a :: bases => let tested := trial n a if tested.1 then let rest := run n bases (rest.1, tested.2 + rest.2) else (false, tested.2)
theorem run_spec (n : Nat) (bases : List Nat) : (run n bases).1 = true ↔ ∀ a ∈ bases, strongPseudoprime n a := by induction bases with | nil => simp [run] | cons a bases ih => simp only [run] cases ht : (trial n a).1 with | false => have hn : ¬ strongPseudoprime n a := by intro h; have := (trial_spec n a).mpr h; simp [ht] at this simp [hn] | true => simp [ih, (trial_spec n a).mp ht]theorem run_reject (n a : Nat) (bases : List Nat) (h : (trial n a).1 = false) : run n (a :: bases) = (false, (trial n a).2) := by simp [run, h]theorem run_count_le (n : Nat) (bases : List Nat) : (run n bases).2 ≤ bases.length * (3 * Nat.size (n - 1)) := by induction bases with | nil => simp [run] | cons a bases ih => have ht := trial_count_le_size n a simp only [run] split · simp only [List.length_cons]; nlinarith · simp only [List.length_cons]; nlinarith theorem run_eq_legacy_value (n : Nat) (bases : List Nat) : (run n bases).1 = (millerRabinLoop n bases).1 := by apply Bool.eq_iff_iff.mpr rw [run_spec, millerRabinLoop_fst_iff] simp [millerRabin]end CLRS.Chapter31.MillerRabinExecution

CLRSLean.FourthEdition.Chapter_31.Section_31_8_Primality_Testing.Probability

Independent uniform Miller–Rabin trials

The sample space is the full finite product of nonzero natural residues, with replacement. Acceptance is defined by the counted execution, including its early rejection behavior. A coordinatewise bijection counts accepting tuples exactly. Uniform probability is the rational ratio of accepting samples to all samples; the one-base strong-liar theorem then gives the repeated-trial error bound. This is a finite probability model, not an implementation of random sampling.

namespace CLRS.Chapter31.MillerRabinExecution

The full product space: each coordinate is a uniform base in 1,...,n-1.

abbrev Samples (n rounds : Nat) := Fin rounds → Fin (n - 1)
def sampleBases (ω : Samples n rounds) : List Nat := List.ofFn (fun i => (ω i).val + 1)

Acceptance is the Boolean returned by the counted, short-circuiting execution.

abbrev Accepted (n rounds : Nat) := {ω : Samples n rounds // (run n (sampleBases ω)).1 = true}
abbrev Liars (n : Nat) := {a : Fin (n - 1) // strongPseudoprime n (a.val + 1)}theorem sample_accept_iff (ω : Samples n rounds) : (run n (sampleBases ω)).1 = true ↔ ∀ i, strongPseudoprime n ((ω i).val + 1) := by simp [run_spec, sampleBases]

Coordinatewise restriction is a bijection, so no independence premise is supplied by callers.

def acceptedEquiv (n rounds : Nat) : Accepted n rounds ≃ (Fin rounds → Liars n) where toFun ω i := ⟨ω.val i, (sample_accept_iff ω.val).mp ω.property i⟩ invFun f := ⟨fun i => (f i).val, (sample_accept_iff _).mpr (fun i => (f i).property)⟩ left_inv ω := by rfl right_inv f := by rfl
theorem accepted_card (n rounds : Nat) : Nat.card (Accepted n rounds) = Nat.card (Liars n) ^ rounds := by rw [Nat.card_congr (acceptedEquiv n rounds), Nat.card_fun] simptheorem samples_card (n rounds : Nat) : Nat.card (Samples n rounds) = (n - 1) ^ rounds := by simp [Samples]

Uniform probability on the explicitly finite product space, expressed as a rational count ratio. No random generator, entropy source, or sampling-operation runtime is asserted.

noncomputable def uniformError (n rounds : Nat) : ℚ := (Nat.card (Accepted n rounds) : ℚ) / Nat.card (Samples n rounds)
theorem uniformError_eq (n rounds : Nat) : uniformError n rounds = ((Nat.card (Liars n) : ℚ) / (n - 1 : Nat)) ^ rounds := by rw [uniformError, accepted_card, samples_card] push_cast rw [div_pow]

Actual executed false acceptance under independent uniform bases is at most 4^-rounds.

theorem uniform_error_le {n : Nat} (hn : 1 < n) (hodd : Odd n) (hcomp : ¬ Nat.Prime n) (rounds : Nat) : uniformError n rounds ≤ (1 / 4 : ℚ) ^ rounds := by letI : NeZero n := ⟨by omega⟩ have hc := strongLiars_nat_card_le hn hodd hcomp have hprod : Nat.card (Liars n) * 4 ≤ n - 1 := (Nat.mul_le_mul_right 4 hc).trans (Nat.div_mul_le_self (n - 1) 4) have hnpos : (0 : ℚ) < (n - 1 : Nat) := by exact_mod_cast (show 0 < n - 1 by omega) have hratio : (Nat.card (Liars n) : ℚ) / (n - 1 : Nat) ≤ 1 / 4 := by apply (div_le_div_iff₀ hnpos (by norm_num)).mpr simp only [one_mul] exact_mod_cast hprod rw [uniformError_eq] exact pow_le_pow_left₀ (by positivity) hratio rounds
@[simp] theorem uniformError_zero (n : Nat) : uniformError n 0 = 1 := by simp [uniformError_eq]end CLRS.Chapter31.MillerRabinExecution

Scope and implementation notes

Imports
import CLRSLean.Chapter_31 import CLRSLean.FourthEdition.Chapter_31.Section_31_1_Elementary_Number_Theory import CLRSLean.FourthEdition.Chapter_31.Section_31_2_Greatest_Common_Divisor import CLRSLean.FourthEdition.Chapter_31.Section_31_3_Modular_Arithmetic import CLRSLean.FourthEdition.Chapter_31.Section_31_4_Solving_Modular_Linear_Equations import CLRSLean.FourthEdition.Chapter_31.Section_31_5_Chinese_Remainder_Theorem import CLRSLean.FourthEdition.Chapter_31.Section_31_6_Powers_Of_An_Element import CLRSLean.FourthEdition.Chapter_31.Section_31_7_RSA import CLRSLean.FourthEdition.Chapter_31.Section_31_8_Primality_Testing import CLRSLean.FourthEdition.Chapter_31.Section_31_2_Greatest_Common_Divisor.Execution import CLRSLean.FourthEdition.Chapter_31.Section_31_7_RSA.KeyRoundTrip import CLRSLean.FourthEdition.Chapter_31.Section_31_8_Primality_Testing.Probability

Current source

Sections 31.1--31.8 are native fourth-edition sections (elementary number-theoretic notions, the greatest common divisor, modular arithmetic, solving modular linear equations, the Chinese remainder theorem, powers of an element, the RSA public-key cryptosystem, and primality testing), imported directly under Chapter 31. The integer-factorization development (legacy Section 31.9) is retained as supplementary online material (reachable through CLRSLean.OnlineMaterial). Declarations keep their current namespaces; the third-edition-numbered imports CLRSLean.Chapter_31 and CLRSLean.Chapter_31.Section_31_* forward to these sources.

Implementation details

The supporting implementation pages remain available outside the main sidebar:

Coverage boundary

The native sections supply the represented fourth-edition number-theoretic sections (§31.1--31.8), including the executable cost layers: CLRS.Chapter31.modularLinearEquationSolver (§31.4), CLRS.Chapter31.modExpWithCount (§31.6), CLRS.Chapter31.rsaKeyGen and CLRS.Chapter31.rsaEncrypt/CLRS.Chapter31.rsaDecrypt (§31.7), and CLRS.Chapter31.MillerRabinExecution.run (§31.8). §31.1 also carries the least-common-multiple layer (CLRS.Chapter31.gcd_mul_lcm_eq, CLRS.Chapter31.lcm_eq_mul_of_coprime), and §31.5 packages the Chinese remainder theorem as the ring isomorphism CLRS.Chapter31.zmod_chineseRemainder.

The counted Euclid companion refines the public first-argument recursion and proves that its actual division count is euclidDivisions b a for the call euclid a b. Generated RSA keys from supplied distinct primes round-trip all messages modulo the generated modulus; in-range messages recover exactly. No key-prime generation or RSA key-assembly runtime is claimed.

Miller–Rabin now decides from the residues of the same counted exponentiation and squaring execution. It stops at the first rejecting base and uses at most 3 * bases.length * Nat.size (n - 1) modular multiplications. Parameter decomposition, comparisons, supplied-base sampling and bit runtime are outside that counter. The legacy loop retains its detached exponentiation budget. The finite product sample space proves actual executed false acceptance is at most (1 / 4)^rounds for odd composite inputs greater than one, under independent uniform sampling with replacement from residues 1,...,n-1. The Carmichael theorem for 561 establishes membership, not minimality.

See docs/clrs-fourth-edition-map.csv for the section-level mapping and docs/migrations/clrs4.md for compatibility and deprecation policy.

CLRS, fourth edition · Chapter 31 of 35