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 Scott 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 mersenne_pos {p : ℕ} (h : 0 < p) : 0 < mersenne p

theorem one_lt_mersenne {p : ℕ} (hp : 1 < p) : 1 < mersenne 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 (Nat.succ i) = LucasLehmer.s i ^ 2 - 2

def LucasLehmer.sZMod (p : ℕ) : ℕ → 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 (Nat.succ i) = 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 (Nat.succ i) = (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 ≤ LucasLehmer.sMod p i

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

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

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

theorem LucasLehmer.Int.coe_nat_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 : ℕ) : LucasLehmer.sZMod p i = ↑(LucasLehmer.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 ↔ LucasLehmer.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 = { val := Nat.minFac (mersenne p), property := (_ : 0 < Nat.minFac (mersenne p)) }

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.instInhabitedX {q : ℕ+} : Inhabited (LucasLehmer.X q)

Equations

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

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

Equations

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

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

Equations

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

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

Equations

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

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

@[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 : LucasLehmer.X q) (y : LucasLehmer.X q) : (x + y).1 = x.1 + y.1

@[simp]

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

@[simp]

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

@[simp]

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

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

Equations

• LucasLehmer.X.instMulX = { 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 : LucasLehmer.X q) (y : LucasLehmer.X q) : (x * y).1 = x.1 * y.1 + 3 * x.2 * y.2

@[simp]

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

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

Equations

• LucasLehmer.X.instOneX = { 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.instMonoidX {q : ℕ+} : Monoid (LucasLehmer.X q)

Equations

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

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

Equations

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

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

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

Equations

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

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

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

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

Equations

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

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

Equations

• LucasLehmer.X.instCommRingX = let src := inferInstanceAs (Ring (LucasLehmer.X q)); CommRing.mk (_ : ∀ (x x_1 : LucasLehmer.X q), x * x_1 = x_1 * x)

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

Equations

• (_ : Nontrivial (LucasLehmer.X q)) = (_ : Nontrivial (LucasLehmer.X q))

@[simp]

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

@[simp]

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

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

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

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

The cardinality of X is q^2.

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

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

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

We define ω = 2 + √3.

Equations

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

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

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

Equations

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

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

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

theorem LucasLehmer.X.closed_form {q : ℕ+} (i : ℕ) : ↑(LucasLehmer.s i) = LucasLehmer.X.ω ^ 2 ^ i + LucasLehmer.X.ω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 < LucasLehmer.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 : ℤ), LucasLehmer.X.ω ^ 2 ^ (p' + 1) = ↑k * ↑(mersenne (p' + 2)) * LucasLehmer.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) : LucasLehmer.X.ω ^ 2 ^ (p' + 1) = -1

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

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

ω as an element of the group of units.

Equations

• LucasLehmer.ωUnit p = { val := LucasLehmer.X.ω, inv := LucasLehmer.X.ωb, val_inv := (_ : LucasLehmer.X.ω * LucasLehmer.X.ωb = 1), inv_val := (_ : LucasLehmer.X.ωb * LucasLehmer.X.ω = 1) }

@[simp]

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

theorem LucasLehmer.order_ω (p' : ℕ) (h : lucasLehmerResidue (p' + 2) = 0) : orderOf (LucasLehmer.ω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) < ↑(LucasLehmer.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.sMod' (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.sMod' q 0 = 4 % q
• LucasLehmer.norm_num_ext.sMod' q (Nat.succ i) = (LucasLehmer.norm_num_ext.sMod' q i ^ 2 + (q - 2)) % q

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

theorem LucasLehmer.norm_num_ext.testTrueHelper (p : ℕ) (hp : Nat.blt 1 p = true) (h : LucasLehmer.norm_num_ext.sMod' (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 (LucasLehmer.norm_num_ext.sMod' (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]