The Lucas-Lehmer test for Mersenne primes.

We define lucasLehmerResidue : Π p : ℕ, ZMod (2^p - 1), and prove lucasLehmerResidue p = 0 → Prime (mersenne p).

We construct a norm_num extension to calculate this residue to certify primality of Mersenne primes using lucas_lehmer_sufficiency.

TODO

• Show reverse implication.
• Speed up the calculations using n ≡ (n % 2^p) + (n / 2^p) [MOD 2^p - 1].
• Find some bigger primes!

History

This development began as a student project by Ainsley Pahljina, and was then cleaned up for mathlib by Kim Morrison. The tactic for certified computation of Lucas-Lehmer residues was provided by Mario Carneiro. This tactic was ported by Thomas Murrills to Lean 4, and then it was converted to a norm_num extension and made to use kernel reductions by Kyle Miller.

def mersenne (p : ℕ) : ℕ

The Mersenne numbers, 2^p - 1.

Equations

• mersenne p = 2 ^ p - 1

theorem strictMono_mersenne : StrictMono mersenne

@[simp]

theorem mersenne_lt_mersenne {p q : ℕ} : mersenne p < mersenne q ↔ p < q

theorem GCongr.mersenne_lt_mersenne {p q : ℕ} : p < q → mersenne p < mersenne q

Alias of the reverse direction of mersenne_lt_mersenne.

@[simp]

theorem mersenne_le_mersenne {p q : ℕ} : mersenne p ≤ mersenne q ↔ p ≤ q

theorem GCongr.mersenne_le_mersenne {p q : ℕ} : p ≤ q → mersenne p ≤ mersenne q

Alias of the reverse direction of mersenne_le_mersenne.

@[simp]

theorem mersenne_zero : mersenne 0 = 0

@[simp]

theorem mersenne_odd {p : ℕ} : Odd (mersenne p) ↔ p ≠ 0

@[simp]

theorem mersenne_pos {p : ℕ} : 0 < mersenne p ↔ 0 < p

theorem Mathlib.Meta.Positivity.mersenne_pos_of_pos {p : ℕ} : 0 < p → 0 < mersenne p

Alias of the reverse direction of mersenne_pos.

def Mathlib.Meta.Positivity.evalMersenne : PositivityExt

Extension for the positivity tactic: mersenne.

Equations

• One or more equations did not get rendered due to their size.

@[simp]

theorem one_lt_mersenne {p : ℕ} : 1 < mersenne p ↔ 1 < p

@[simp]

theorem succ_mersenne (k : ℕ) : mersenne k + 1 = 2 ^ k

We now define three(!) different versions of the recurrence s (i+1) = (s i)^2 - 2.

These versions take values either in ℤ, in ZMod (2^p - 1), or in ℤ but applying % (2^p - 1) at each step.

They are each useful at different points in the proof, so we take a moment setting up the lemmas relating them.

def LucasLehmer.s : ℕ → ℤ

The recurrence s (i+1) = (s i)^2 - 2 in ℤ.

Equations

• LucasLehmer.s 0 = 4
• LucasLehmer.s i.succ = LucasLehmer.s i ^ 2 - 2

def LucasLehmer.sZMod (p a✝ : ℕ) : ZMod (2 ^ p - 1)

The recurrence s (i+1) = (s i)^2 - 2 in ZMod (2^p - 1).

Equations

• LucasLehmer.sZMod p 0 = 4
• LucasLehmer.sZMod p i.succ = LucasLehmer.sZMod p i ^ 2 - 2

def LucasLehmer.sMod (p : ℕ) : ℕ → ℤ

The recurrence s (i+1) = ((s i)^2 - 2) % (2^p - 1) in ℤ.

Equations

• LucasLehmer.sMod p 0 = 4 % (2 ^ p - 1)
• LucasLehmer.sMod p i.succ = (LucasLehmer.sMod p i ^ 2 - 2) % (2 ^ p - 1)

theorem LucasLehmer.mersenne_int_pos {p : ℕ} (hp : p ≠ 0) : 0 < 2 ^ p - 1

theorem LucasLehmer.mersenne_int_ne_zero (p : ℕ) (hp : p ≠ 0) : 2 ^ p - 1 ≠ 0

theorem LucasLehmer.sMod_nonneg (p : ℕ) (hp : p ≠ 0) (i : ℕ) : 0 ≤ sMod p i

theorem LucasLehmer.sMod_mod (p i : ℕ) : sMod p i % (2 ^ p - 1) = sMod p i

theorem LucasLehmer.sMod_lt (p : ℕ) (hp : p ≠ 0) (i : ℕ) : sMod p i < 2 ^ p - 1

theorem LucasLehmer.sZMod_eq_s (p' i : ℕ) : sZMod (p' + 2) i = ↑(s i)

theorem LucasLehmer.Int.natCast_pow_pred (b p : ℕ) (w : 0 < b) : ↑(b ^ p - 1) = ↑b ^ p - 1

theorem LucasLehmer.Int.coe_nat_two_pow_pred (p : ℕ) : ↑(2 ^ p - 1) = 2 ^ p - 1

theorem LucasLehmer.sZMod_eq_sMod (p i : ℕ) : sZMod p i = ↑(sMod p i)

def LucasLehmer.lucasLehmerResidue (p : ℕ) : ZMod (2 ^ p - 1)

The Lucas-Lehmer residue is s p (p-2) in ZMod (2^p - 1).

Equations

• lucasLehmerResidue p = LucasLehmer.sZMod p (p - 2)

theorem LucasLehmer.residue_eq_zero_iff_sMod_eq_zero (p : ℕ) (w : 1 < p) : lucasLehmerResidue p = 0 ↔ sMod p (p - 2) = 0

def LucasLehmer.LucasLehmerTest (p : ℕ) : Prop

Lucas-Lehmer Test: a Mersenne number 2^p-1 is prime if and only if the Lucas-Lehmer residue s p (p-2) % (2^p - 1) is zero.

Equations

• LucasLehmerTest p = (lucasLehmerResidue p = 0)

def LucasLehmer.q (p : ℕ) : ℕ+

q is defined as the minimum factor of mersenne p, bundled as an ℕ+.

Equations

• LucasLehmer.q p = ⟨(mersenne p).minFac, ⋯⟩

def LucasLehmer.X (q : ℕ+) : Type

We construct the ring X q as ℤ/qℤ + √3 ℤ/qℤ.

Equations

• LucasLehmer.X q = (ZMod ↑q × ZMod ↑q)

instance LucasLehmer.X.instInhabited {q : ℕ+} : Inhabited (X q)

Equations

• LucasLehmer.X.instInhabited = inferInstanceAs (Inhabited (ZMod ↑q × ZMod ↑q))

instance LucasLehmer.X.instFintype {q : ℕ+} : Fintype (X q)

Equations

• LucasLehmer.X.instFintype = inferInstanceAs (Fintype (ZMod ↑q × ZMod ↑q))

instance LucasLehmer.X.instDecidableEq {q : ℕ+} : DecidableEq (X q)

Equations

• LucasLehmer.X.instDecidableEq = inferInstanceAs (DecidableEq (ZMod ↑q × ZMod ↑q))

instance LucasLehmer.X.instAddCommGroup {q : ℕ+} : AddCommGroup (X q)

Equations

• LucasLehmer.X.instAddCommGroup = inferInstanceAs (AddCommGroup (ZMod ↑q × ZMod ↑q))

theorem LucasLehmer.X.ext {q : ℕ+} {x y : X q} (h₁ : x.1 = y.1) (h₂ : x.2 = y.2) : x = y

theorem LucasLehmer.X.ext_iff {q : ℕ+} {x y : X q} : x = y ↔ x.1 = y.1 ∧ x.2 = y.2

@[simp]

theorem LucasLehmer.X.zero_fst {q : ℕ+} : 0.1 = 0

@[simp]

theorem LucasLehmer.X.zero_snd {q : ℕ+} : 0.2 = 0

@[simp]

theorem LucasLehmer.X.add_fst {q : ℕ+} (x y : X q) : (x + y).1 = x.1 + y.1

@[simp]

theorem LucasLehmer.X.add_snd {q : ℕ+} (x y : X q) : (x + y).2 = x.2 + y.2

@[simp]

theorem LucasLehmer.X.neg_fst {q : ℕ+} (x : X q) : (-x).1 = -x.1

@[simp]

theorem LucasLehmer.X.neg_snd {q : ℕ+} (x : X q) : (-x).2 = -x.2

instance LucasLehmer.X.instMul {q : ℕ+} : Mul (X q)

Equations

• LucasLehmer.X.instMul = { mul := fun (x y : LucasLehmer.X q) => (x.1 * y.1 + 3 * x.2 * y.2, x.1 * y.2 + x.2 * y.1) }

@[simp]

theorem LucasLehmer.X.mul_fst {q : ℕ+} (x y : X q) : (x * y).1 = x.1 * y.1 + 3 * x.2 * y.2

@[simp]

theorem LucasLehmer.X.mul_snd {q : ℕ+} (x y : X q) : (x * y).2 = x.1 * y.2 + x.2 * y.1

instance LucasLehmer.X.instOne {q : ℕ+} : One (X q)

Equations

• LucasLehmer.X.instOne = { one := (1, 0) }

@[simp]

theorem LucasLehmer.X.one_fst {q : ℕ+} : 1.1 = 1

@[simp]

theorem LucasLehmer.X.one_snd {q : ℕ+} : 1.2 = 0

instance LucasLehmer.X.instMonoid {q : ℕ+} : Monoid (X q)

Equations

• One or more equations did not get rendered due to their size.

instance LucasLehmer.X.instNatCast {q : ℕ+} : NatCast (X q)

Equations

• LucasLehmer.X.instNatCast = { natCast := fun (n : ℕ) => (↑n, 0) }

@[simp]

theorem LucasLehmer.X.fst_natCast {q : ℕ+} (n : ℕ) : (↑n).1 = ↑n

@[simp]

theorem LucasLehmer.X.snd_natCast {q : ℕ+} (n : ℕ) : (↑n).2 = 0

@[simp]

theorem LucasLehmer.X.ofNat_fst {q : ℕ+} (n : ℕ) [n.AtLeastTwo] : (OfNat.ofNat n).1 = OfNat.ofNat n

@[simp]

theorem LucasLehmer.X.ofNat_snd {q : ℕ+} (n : ℕ) [n.AtLeastTwo] : (OfNat.ofNat n).2 = 0

instance LucasLehmer.X.instAddGroupWithOne {q : ℕ+} : AddGroupWithOne (X q)

Equations

• One or more equations did not get rendered due to their size.

theorem LucasLehmer.X.left_distrib {q : ℕ+} (x y z : X q) : x * (y + z) = x * y + x * z

theorem LucasLehmer.X.right_distrib {q : ℕ+} (x y z : X q) : (x + y) * z = x * z + y * z

instance LucasLehmer.X.instRing {q : ℕ+} : Ring (X q)

Equations

• One or more equations did not get rendered due to their size.

instance LucasLehmer.X.instCommRing {q : ℕ+} : CommRing (X q)

Equations

• LucasLehmer.X.instCommRing = { toRing := inferInstanceAs (Ring (LucasLehmer.X q)), mul_comm := ⋯ }

instance LucasLehmer.X.instNontrivialOfFactLtNatOfNatVal {q : ℕ+} [Fact (1 < ↑q)] : Nontrivial (X q)

@[simp]

theorem LucasLehmer.X.fst_intCast {q : ℕ+} (n : ℤ) : (↑n).1 = ↑n

@[simp]

theorem LucasLehmer.X.snd_intCast {q : ℕ+} (n : ℤ) : (↑n).2 = 0

theorem LucasLehmer.X.coe_mul {q : ℕ+} (n m : ℤ) : ↑(n * m) = ↑n * ↑m

theorem LucasLehmer.X.coe_natCast {q : ℕ+} (n : ℕ) : ↑↑n = ↑n

theorem LucasLehmer.X.card_eq {q : ℕ+} : Fintype.card (X q) = ↑q ^ 2

The cardinality of X is q^2.

theorem LucasLehmer.X.card_units_lt {q : ℕ+} (w : 1 < q) : Fintype.card (X q)ˣ < ↑q ^ 2

There are strictly fewer than q^2 units, since 0 is not a unit.

def LucasLehmer.X.ω {q : ℕ+} : X q

We define ω = 2 + √3.

Equations

• LucasLehmer.X.ω = (2, 1)

def LucasLehmer.X.ωb {q : ℕ+} : X q

We define ωb = 2 - √3, which is the inverse of ω.

Equations

• LucasLehmer.X.ωb = (2, -1)

theorem LucasLehmer.X.ω_mul_ωb (q : ℕ+) : ω * ωb = 1

theorem LucasLehmer.X.ωb_mul_ω (q : ℕ+) : ωb * ω = 1

theorem LucasLehmer.X.closed_form {q : ℕ+} (i : ℕ) : ↑(s i) = ω ^ 2 ^ i + ωb ^ 2 ^ i

A closed form for the recurrence relation.

Here and below, we introduce p' = p - 2, in order to avoid using subtraction in ℕ.

theorem LucasLehmer.two_lt_q (p' : ℕ) : 2 < q (p' + 2)

If 1 < p, then q p, the smallest prime factor of mersenne p, is more than 2.

theorem LucasLehmer.ω_pow_formula (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : ∃ (k : ℤ), X.ω ^ 2 ^ (p' + 1) = ↑k * ↑(mersenne (p' + 2)) * X.ω ^ 2 ^ p' - 1

theorem LucasLehmer.mersenne_coe_X (p : ℕ) : ↑(mersenne p) = 0

q is the minimum factor of mersenne p, so M p = 0 in X q.

theorem LucasLehmer.ω_pow_eq_neg_one (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : X.ω ^ 2 ^ (p' + 1) = -1

theorem LucasLehmer.ω_pow_eq_one (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : X.ω ^ 2 ^ (p' + 2) = 1

def LucasLehmer.ωUnit (p : ℕ) : (X (q p))ˣ

ω as an element of the group of units.

Equations

• LucasLehmer.ωUnit p = { val := LucasLehmer.X.ω, inv := LucasLehmer.X.ωb, val_inv := ⋯, inv_val := ⋯ }

@[simp]

theorem LucasLehmer.ωUnit_coe (p : ℕ) : ↑(ωUnit p) = X.ω

theorem LucasLehmer.order_ω (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : orderOf (ωUnit (p' + 2)) = 2 ^ (p' + 2)

The order of ω in the unit group is exactly 2^p.

theorem LucasLehmer.order_ineq (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : 2 ^ (p' + 2) < ↑(q (p' + 2)) ^ 2

theorem lucas_lehmer_sufficiency (p : ℕ) (w : 1 < p) : LucasLehmerTest p → Nat.Prime (mersenne p)

norm_num extension

Next we define a norm_num extension that calculates LucasLehmerTest p for 1 < p. It makes use of a version of sMod that is specifically written to be reducible by the Lean 4 kernel, which has the capability of efficiently reducing natural number expressions. With this reduction in hand, it's a simple matter of applying the lemma LucasLehmer.residue_eq_zero_iff_sMod_eq_zero.

See [Archive/Examples/MersennePrimes.lean] for certifications of all Mersenne primes up through mersenne 4423.

def LucasLehmer.norm_num_ext.sModNat (q : ℕ) : ℕ → ℕ

Version of sMod that is ℕ-valued. One should have q = 2 ^ p - 1. This can be reduced by the kernel.

Equations

• LucasLehmer.norm_num_ext.sModNat q 0 = 4 % q
• LucasLehmer.norm_num_ext.sModNat q i.succ = (LucasLehmer.norm_num_ext.sModNat q i ^ 2 + (q - 2)) % q

theorem LucasLehmer.norm_num_ext.sModNat_eq_sMod (p k : ℕ) (hp : 2 ≤ p) : ↑(sModNat (2 ^ p - 1) k) = sMod p k

def LucasLehmer.norm_num_ext.sModNatTR (q k : ℕ) : ℕ

Tail-recursive version of sModNat.

Equations

• LucasLehmer.norm_num_ext.sModNatTR q k = LucasLehmer.norm_num_ext.sModNatTR.go q k (4 % q)

def LucasLehmer.norm_num_ext.sModNatTR.go (q : ℕ) : ℕ → ℕ → ℕ

Helper function for sMod''.

Equations

• LucasLehmer.norm_num_ext.sModNatTR.go q 0 x✝ = x✝
• LucasLehmer.norm_num_ext.sModNatTR.go q n.succ x✝ = LucasLehmer.norm_num_ext.sModNatTR.go q n ((x✝ ^ 2 + (q - 2)) % q)

def LucasLehmer.norm_num_ext.sModNat_aux (b q : ℕ) : ℕ → ℕ

Generalization of sModNat with arbitrary base case, useful for proving sModNatTR and sModNat agree.

Equations

• LucasLehmer.norm_num_ext.sModNat_aux b q 0 = b
• LucasLehmer.norm_num_ext.sModNat_aux b q i.succ = (LucasLehmer.norm_num_ext.sModNat_aux b q i ^ 2 + (q - 2)) % q

theorem LucasLehmer.norm_num_ext.sModNat_aux_eq (q k : ℕ) : sModNat_aux (4 % q) q k = sModNat q k

theorem LucasLehmer.norm_num_ext.sModNatTR_eq_sModNat (q i : ℕ) : sModNatTR q i = sModNat q i

theorem LucasLehmer.norm_num_ext.sModNatTR_eq_sModNat.helper (b q k : ℕ) : sModNatTR.go q k b = sModNat_aux b q k

theorem LucasLehmer.norm_num_ext.testTrueHelper (p : ℕ) (hp : Nat.blt 1 p = true) (h : sModNatTR (2 ^ p - 1) (p - 2) = 0) : LucasLehmerTest p

theorem LucasLehmer.norm_num_ext.testFalseHelper (p : ℕ) (hp : Nat.blt 1 p = true) (h : Nat.ble 1 (sModNatTR (2 ^ p - 1) (p - 2)) = true) : ¬LucasLehmerTest p

theorem LucasLehmer.norm_num_ext.isNat_lucasLehmerTest {p np : ℕ} : Mathlib.Meta.NormNum.IsNat p np → LucasLehmerTest np → LucasLehmerTest p

theorem LucasLehmer.norm_num_ext.isNat_not_lucasLehmerTest {p np : ℕ} : Mathlib.Meta.NormNum.IsNat p np → ¬LucasLehmerTest np → ¬LucasLehmerTest p

def LucasLehmer.norm_num_ext.evalLucasLehmerTest : Mathlib.Meta.NormNum.NormNumExt

Calculate LucasLehmer.LucasLehmerTest p for 2 ≤ p by using kernel reduction for the sMod' function.

Equations

• One or more equations did not get rendered due to their size.

This implementation works successfully to prove (2^4423 - 1).Prime, and all the Mersenne primes up to this point appear in [Archive/Examples/MersennePrimes.lean]. These can be calculated nearly instantly, and (2^9689 - 1).Prime only fails due to deep recursion.

(Note by kmill: the following notes were for the Lean 3 version. They seem like they could still be useful, so I'm leaving them here.)

There's still low hanging fruit available to do faster computations based on the formula

n ≡ (n % 2^p) + (n / 2^p) [MOD 2^p - 1]

and the fact that % 2^p and / 2^p can be very efficient on the binary representation. Someone should do this, too!

theorem modEq_mersenne (n k : ℕ) : k ≡ k / 2 ^ n + k % 2 ^ n [MOD 2 ^ n - 1]