# mathlib3documentation

number_theory.lucas_lehmer

# The Lucas-Lehmer test for Mersenne primes. #

THIS FILE IS SYNCHRONIZED WITH MATHLIB4. Any changes to this file require a corresponding PR to mathlib4.

We define lucas_lehmer_residue : Π p : ℕ, zmod (2^p - 1), and prove lucas_lehmer_residue p = 0 → prime (mersenne p).

We construct a tactic lucas_lehmer.run_test, which iteratively certifies the arithmetic required to calculate the residue, and enables us to prove

example : prime (mersenne 127) :=
lucas_lehmer_sufficiency _ (by norm_num) (by lucas_lehmer.run_test)


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

def mersenne (p : ) :

The Mersenne numbers, 2^p - 1.

Equations
theorem mersenne_pos {p : } (h : 0 < p) :
0 <
@[simp]
theorem succ_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 lucas_lehmer.s  :

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

Equations
def lucas_lehmer.s_zmod (p : ) :
zmod (2 ^ p - 1)

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

Equations
def lucas_lehmer.s_mod (p : ) :

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

Equations
theorem lucas_lehmer.mersenne_int_ne_zero (p : ) (w : 0 < p) :
2 ^ p - 1 0
theorem lucas_lehmer.s_mod_nonneg (p : ) (w : 0 < p) (i : ) :
0
theorem lucas_lehmer.s_mod_mod (p i : ) :
% (2 ^ p - 1) =
theorem lucas_lehmer.s_mod_lt (p : ) (w : 0 < p) (i : ) :
< 2 ^ p - 1
theorem lucas_lehmer.int.coe_nat_pow_pred (b p : ) (w : 0 < b) :
(b ^ p - 1) = b ^ p - 1
theorem lucas_lehmer.int.coe_nat_two_pow_pred (p : ) :
(2 ^ p - 1) = 2 ^ p - 1
theorem lucas_lehmer.s_zmod_eq_s_mod (p i : ) :
= i)

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

Equations
theorem lucas_lehmer.residue_eq_zero_iff_s_mod_eq_zero (p : ) (w : 1 < p) :
(p - 2) = 0
@[protected, instance]

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
Instances for lucas_lehmer.lucas_lehmer_test
def lucas_lehmer.q (p : ) :

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

Equations
@[protected, instance]
@[protected, instance]
def lucas_lehmer.X (q : ℕ+) :

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

Equations
Instances for lucas_lehmer.X
@[protected, instance]
@[protected, instance]
@[ext]
theorem lucas_lehmer.X.ext {q : ℕ+} {x y : lucas_lehmer.X q} (h₁ : x.fst = y.fst) (h₂ : x.snd = y.snd) :
x = y
@[simp]
theorem lucas_lehmer.X.add_fst {q : ℕ+} (x y : lucas_lehmer.X q) :
(x + y).fst = x.fst + y.fst
@[simp]
theorem lucas_lehmer.X.add_snd {q : ℕ+} (x y : lucas_lehmer.X q) :
(x + y).snd = x.snd + y.snd
@[simp]
theorem lucas_lehmer.X.neg_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(-x).fst = -x.fst
@[simp]
theorem lucas_lehmer.X.neg_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(-x).snd = -x.snd
@[protected, instance]
Equations
@[simp]
theorem lucas_lehmer.X.mul_fst {q : ℕ+} (x y : lucas_lehmer.X q) :
(x * y).fst = x.fst * y.fst + 3 * x.snd * y.snd
@[simp]
theorem lucas_lehmer.X.mul_snd {q : ℕ+} (x y : lucas_lehmer.X q) :
(x * y).snd = x.fst * y.snd + x.snd * y.fst
@[protected, instance]
Equations
@[simp]
theorem lucas_lehmer.X.one_fst {q : ℕ+} :
1.fst = 1
@[simp]
theorem lucas_lehmer.X.one_snd {q : ℕ+} :
1.snd = 0
@[simp]
theorem lucas_lehmer.X.bit0_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(bit0 x).fst = bit0 x.fst
@[simp]
theorem lucas_lehmer.X.bit0_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(bit0 x).snd = bit0 x.snd
@[simp]
theorem lucas_lehmer.X.bit1_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(bit1 x).fst = bit1 x.fst
@[simp]
theorem lucas_lehmer.X.bit1_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(bit1 x).snd = bit0 x.snd
@[protected, instance]
Equations
@[protected, instance]
Equations
theorem lucas_lehmer.X.left_distrib {q : ℕ+} (x y z : lucas_lehmer.X q) :
x * (y + z) = x * y + x * z
theorem lucas_lehmer.X.right_distrib {q : ℕ+} (x y z : lucas_lehmer.X q) :
(x + y) * z = x * z + y * z
@[protected, instance]
def lucas_lehmer.X.ring {q : ℕ+} :
Equations
@[protected, instance]
Equations
@[protected, instance]
def lucas_lehmer.X.nontrivial {q : ℕ+} [fact (1 < q)] :
@[simp]
theorem lucas_lehmer.X.nat_coe_fst {q : ℕ+} (n : ) :
@[simp]
theorem lucas_lehmer.X.nat_coe_snd {q : ℕ+} (n : ) :
n.snd = 0
@[simp]
theorem lucas_lehmer.X.int_coe_fst {q : ℕ+} (n : ) :
@[simp]
theorem lucas_lehmer.X.int_coe_snd {q : ℕ+} (n : ) :
n.snd = 0
@[norm_cast]
theorem lucas_lehmer.X.coe_mul {q : ℕ+} (n m : ) :
(n * m) = n * m
@[norm_cast]
theorem lucas_lehmer.X.coe_nat {q : ℕ+} (n : ) :
theorem lucas_lehmer.X.X_card {q : ℕ+} :
= q ^ 2

The cardinality of X is q^2.

theorem lucas_lehmer.X.units_card {q : ℕ+} (w : 1 < q) :
< q ^ 2

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

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

We define ω = 2 + √3.

Equations
def lucas_lehmer.X.ωb {q : ℕ+} :

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

Equations
theorem lucas_lehmer.X.closed_form {q : ℕ+} (i : ) :

A closed form for the recurrence relation.

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

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

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

theorem lucas_lehmer.ω_pow_formula (p' : ) (h : = 0) :
(k : ), lucas_lehmer.X.ω ^ 2 ^ (p' + 1) = k * (mersenne (p' + 2)) * lucas_lehmer.X.ω ^ 2 ^ p' - 1
theorem lucas_lehmer.mersenne_coe_X (p : ) :

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

theorem lucas_lehmer.ω_pow_eq_neg_one (p' : ) (h : = 0) :
lucas_lehmer.X.ω ^ 2 ^ (p' + 1) = -1
theorem lucas_lehmer.ω_pow_eq_one (p' : ) (h : = 0) :
lucas_lehmer.X.ω ^ 2 ^ (p' + 2) = 1
def lucas_lehmer.ω_unit (p : ) :

ω as an element of the group of units.

Equations
@[simp]
theorem lucas_lehmer.ω_unit_coe (p : ) :
theorem lucas_lehmer.order_ω (p' : ) (h : = 0) :
order_of (lucas_lehmer.ω_unit (p' + 2)) = 2 ^ (p' + 2)

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

theorem lucas_lehmer.order_ineq (p' : ) (h : = 0) :
2 ^ (p' + 2) < (lucas_lehmer.q (p' + 2)) ^ 2
theorem lucas_lehmer_sufficiency (p : ) (w : 1 < p) :
theorem lucas_lehmer.s_mod_succ {p : } {a : } {i : } {b c : } (h1 : 2 ^ p - 1 = a) (h2 : = b) (h3 : (b * b - 2) % a = c) :
(i + 1) = c
meta def lucas_lehmer.run_test  :

Given a goal of the form lucas_lehmer_test p, attempt to do the calculation using norm_num to certify each step.

This implementation works successfully to prove (2^127 - 1).prime, and all the Mersenne primes up to this point appear in [archive/examples/mersenne_primes.lean].

(2^127 - 1).prime takes about 5 minutes to run (depending on your CPU!), and unfortunately the next Mersenne prime (2^521 - 1), which was the first "computer era" prime, is out of reach with the current implementation.

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]