mathlib documentation

number_theory.lucas_lehmer

The Lucas-Lehmer test for Mersenne primes.

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

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  :

The Mersenne numbers, 2^p - 1.

Equations
theorem mersenne_pos {p : } :
0 < p0 < 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 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  :

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

Equations
theorem lucas_lehmer.mersenne_int_ne_zero (p : ) :
0 < p2 ^ p - 1 0

theorem lucas_lehmer.s_mod_nonneg (p : ) (w : 0 < p) (i : ) :

theorem lucas_lehmer.s_mod_lt (p : ) (w : 0 < p) (i : ) :

theorem lucas_lehmer.int.coe_nat_pow_pred (b p : ) :
0 < b(b ^ p - 1) = b ^ p - 1

theorem lucas_lehmer.int.coe_nat_two_pow_pred (p : ) :
(2 ^ p - 1) = 2 ^ p - 1

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

Equations

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
def lucas_lehmer.q  :
ℕ+

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

Equations
@[instance]

def lucas_lehmer.X  :
ℕ+ → Type

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

Equations
@[ext]
theorem lucas_lehmer.X.ext {q : ℕ+} {x y : lucas_lehmer.X q} :
x.fst = y.fstx.snd = y.sndx = 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

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

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

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

@[instance]

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

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

theorem lucas_lehmer.X.coe_nat {q : ℕ+} (n : ) :

The cardinality of X is q^2.

theorem lucas_lehmer.X.units_card {q : ℕ+} :

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

We define ω = 2 + √3.

Equations

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

Equations

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' : ) :
lucas_lehmer_residue (p' + 2) = 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' : ) :
lucas_lehmer_residue (p' + 2) = 0lucas_lehmer.X.ω ^ 2 ^ (p' + 1) = -1

theorem lucas_lehmer.ω_pow_eq_one (p' : ) :
lucas_lehmer_residue (p' + 2) = 0lucas_lehmer.X.ω ^ 2 ^ (p' + 2) = 1

ω as an element of the group of units.

Equations
theorem lucas_lehmer.order_ω (p' : ) :
lucas_lehmer_residue (p' + 2) = 0order_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' : ) :
lucas_lehmer_residue (p' + 2) = 02 ^ (p' + 2) < (lucas_lehmer.q (p' + 2)) ^ 2

@[instance]

@[instance]

theorem lucas_lehmer.s_mod_succ {p : } {a : } {i : } {b c : } :
2 ^ p - 1 = alucas_lehmer.s_mod p i = b(b * b - 2) % a = clucas_lehmer.s_mod p (i + 1) = c

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]