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

• 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.

theorem strictMono_mersenne : StrictMono mersenne

@[simp]

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

@[simp]

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

@[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 mersenne_succ (n : ℕ) : mersenne (n + 1) = 2 * mersenne n + 1

theorem Nat.Prime.of_mersenne {p : ℕ} (h : Prime (mersenne p)) : Prime p

If 2 ^ p - 1 is prime, then p is prime.

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.

@[simp]

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

@[simp]

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

theorem mersenne_mod_four {n : ℕ} (h : 2 ≤ n) : mersenne n % 4 = 3

theorem mersenne_mod_three {n : ℕ} (odd : Odd n) (h : 3 ≤ n) : mersenne n % 3 = 1

theorem mersenne_mod_eight {n : ℕ} (h : 3 ≤ n) : mersenne n % 8 = 7

theorem legendreSym_mersenne_two {p : ℕ} [Fact (Nat.Prime (mersenne p))] (hp : 3 ≤ p) : legendreSym (mersenne p) 2 = 1

If 2^p - 1 is prime then 2 is a square mod 2^p - 1.

theorem legendreSym_mersenne_three {p : ℕ} [Fact (Nat.Prime (mersenne p))] (hp : 3 ≤ p) (odd : Odd p) : legendreSym (mersenne p) 3 = -1

If 2^p - 1 is prime then 3 is not a square mod 2^p - 1.

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 ℤ.

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

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

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

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

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.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).

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.

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

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

def LucasLehmer.X (q : ℕ) : Type

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

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

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

instance LucasLehmer.X.instAddCommGroup {q : ℕ} : AddCommGroup (X 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)

@[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)

@[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)

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

@[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)

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)

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

instance LucasLehmer.X.instNontrivialOfFactLtNatOfNat {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

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

We define ω = 2 + √3.

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

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

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.

def LucasLehmer.X.α {q : ℕ} : X q

We define α = √3.

@[simp]

theorem LucasLehmer.X.α_sq {q : ℕ} : α ^ 2 = 3

@[simp]

theorem LucasLehmer.X.one_add_α_sq {q : ℕ} : (1 + α) ^ 2 = 2 * ω

theorem LucasLehmer.X.α_pow {q : ℕ} (i : ℕ) : α ^ (2 * i + 1) = 3 ^ i * α

We show that X q has characteristic q, so that we can apply the binomial theorem.

instance LucasLehmer.X.instCharP {q : ℕ} : CharP (X q) q

instance LucasLehmer.X.instCoeZMod {q : ℕ} : Coe (ZMod q) (X q)

theorem LucasLehmer.X.one_add_α_pow_q {q : ℕ} [Fact (Nat.Prime q)] (odd : Odd q) (leg3 : legendreSym q 3 = -1) : (1 + α) ^ q = 1 - α

If 3 is not a square mod q then (1 + α) ^ q = 1 - α

theorem LucasLehmer.X.one_add_α_pow_q_succ {q : ℕ} [Fact (Nat.Prime q)] (odd : Odd q) (leg3 : legendreSym q 3 = -1) : (1 + α) ^ (q + 1) = -2

If 3 is not a square then (1 + α) ^ (q + 1) = -2.

theorem LucasLehmer.X.two_mul_ω_pow {q : ℕ} [Fact (Nat.Prime q)] (odd : Odd q) (leg3 : legendreSym q 3 = -1) : (2 * ω) ^ ((q + 1) / 2) = -2

If 3 is not a square then (2 * ω) ^ ((q + 1) / 2) = -2.

theorem LucasLehmer.X.pow_ω {q : ℕ} [Fact (Nat.Prime q)] (odd : Odd q) (leg3 : legendreSym q 3 = -1) (leg2 : legendreSym q 2 = 1) : ω ^ ((q + 1) / 2) = -1

If 3 is not a square and 2 is square then $\omega^{(q+1)/2}=-1$.

theorem LucasLehmer.X.ω_pow_trace {q : ℕ} [Fact (Nat.Prime q)] (odd : Odd q) (leg3 : legendreSym q 3 = -1) (leg2 : legendreSym q 2 = 1) (hq4 : 4 ∣ q + 1) : ω ^ ((q + 1) / 4) + ωb ^ ((q + 1) / 4) = 0

The final evaluation needed to establish the Lucas-Lehmer necessity.

instance LucasLehmer.X.instFintype {q : ℕ} [NeZero q] : Fintype (X q)

theorem LucasLehmer.X.card_eq {q : ℕ} [NeZero q] : Fintype.card (X q) = q ^ 2

The cardinality of X is q^2.

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

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

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.

@[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)

theorem lucas_lehmer_necessity (p : ℕ) (w : 3 ≤ p) (hp : Nat.Prime (mersenne p)) : LucasLehmerTest p

If 2^p - 1 is prime then the Lucas-Lehmer test holds, s (p - 2) % (2^p - 1) = 0.

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.

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.

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

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

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.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.

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]