mathlib3 documentation

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 #

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.

source
def mersenne (p : ℕ) :
ℕ

The Mersenne numbers, 2^p - 1.

Equations
source
theorem mersenne_pos {p : ℕ} (h : 0 < p) :
0 < mersenne p
source
@[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.

source
def lucas_lehmer.s  :
ℕ → ℤ

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

Equations
source
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
source
def lucas_lehmer.s_mod (p : ℕ) :
ℕ → ℤ

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

Equations
source
theorem lucas_lehmer.mersenne_int_ne_zero (p : ℕ) (w : 0 < p) :
2 ^ p - 1 ≠ 0
source
theorem lucas_lehmer.s_mod_nonneg (p : ℕ) (w : 0 < p) (i : ℕ) :
0 ≤ lucas_lehmer.s_mod p i
source
theorem lucas_lehmer.s_mod_mod (p i : ℕ) :
lucas_lehmer.s_mod p i % (2 ^ p - 1) = lucas_lehmer.s_mod p i
source
theorem lucas_lehmer.s_mod_lt (p : ℕ) (w : 0 < p) (i : ℕ) :
lucas_lehmer.s_mod p i < 2 ^ p - 1
source
theorem lucas_lehmer.s_zmod_eq_s (p' i : ℕ) :
lucas_lehmer.s_zmod (p' + 2) i = ↑(lucas_lehmer.s i)
source
theorem lucas_lehmer.int.coe_nat_pow_pred (b p : ℕ) (w : 0 < b) :
↑(b ^ p - 1) = ↑b ^ p - 1
source
theorem lucas_lehmer.int.coe_nat_two_pow_pred (p : ℕ) :
↑(2 ^ p - 1) = 2 ^ p - 1
source
theorem lucas_lehmer.s_zmod_eq_s_mod (p i : ℕ) :
lucas_lehmer.s_zmod p i = ↑(lucas_lehmer.s_mod p i)
source
def lucas_lehmer.lucas_lehmer_residue (p : ℕ) :
zmod (2 ^ p - 1)

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

Equations
source
theorem lucas_lehmer.residue_eq_zero_iff_s_mod_eq_zero (p : ℕ) (w : 1 < p) :
lucas_lehmer.lucas_lehmer_residue p = 0 ↔ lucas_lehmer.s_mod p (p - 2) = 0
source
@[protected, instance]
def lucas_lehmer.lucas_lehmer_test.decidable_pred  :
decidable_pred lucas_lehmer.lucas_lehmer_test
source
def lucas_lehmer.lucas_lehmer_test (p : ℕ) :
Prop

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
source
def lucas_lehmer.q (p : ℕ) :
ℕ+

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

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

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

Equations
Instances for lucas_lehmer.X
source
@[protected, instance]
def lucas_lehmer.X.decidable_eq (q : ℕ+) :
decidable_eq (lucas_lehmer.X q)
source
@[protected, instance]
def lucas_lehmer.X.fintype (q : ℕ+) :
fintype (lucas_lehmer.X q)
source
@[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
source
@[simp]
theorem lucas_lehmer.X.add_fst {q : ℕ+} (x y : lucas_lehmer.X q) :
(x + y).fst = x.fst + y.fst
source
@[simp]
theorem lucas_lehmer.X.add_snd {q : ℕ+} (x y : lucas_lehmer.X q) :
(x + y).snd = x.snd + y.snd
source
@[simp]
theorem lucas_lehmer.X.neg_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(-x).fst = -x.fst
source
@[simp]
theorem lucas_lehmer.X.neg_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(-x).snd = -x.snd
source
@[protected, instance]
def lucas_lehmer.X.has_mul {q : ℕ+} :
has_mul (lucas_lehmer.X q)
Equations
source
@[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
source
@[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
source
@[protected, instance]
def lucas_lehmer.X.has_one {q : ℕ+} :
has_one (lucas_lehmer.X q)
Equations
source
@[simp]
theorem lucas_lehmer.X.one_fst {q : ℕ+} :
1.fst = 1
source
@[simp]
theorem lucas_lehmer.X.one_snd {q : ℕ+} :
1.snd = 0
source
@[simp]
theorem lucas_lehmer.X.bit0_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(bit0 x).fst = bit0 x.fst
source
@[simp]
theorem lucas_lehmer.X.bit0_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(bit0 x).snd = bit0 x.snd
source
@[simp]
theorem lucas_lehmer.X.bit1_fst {q : ℕ+} (x : lucas_lehmer.X q) :
(bit1 x).fst = bit1 x.fst
source
@[simp]
theorem lucas_lehmer.X.bit1_snd {q : ℕ+} (x : lucas_lehmer.X q) :
(bit1 x).snd = bit0 x.snd
source
@[protected, instance]
def lucas_lehmer.X.monoid {q : ℕ+} :
monoid (lucas_lehmer.X q)
Equations
source
@[protected, instance]
def lucas_lehmer.X.add_group_with_one {q : ℕ+} :
add_group_with_one (lucas_lehmer.X q)
Equations
source
theorem lucas_lehmer.X.left_distrib {q : ℕ+} (x y z : lucas_lehmer.X q) :
x * (y + z) = x * y + x * z
source
theorem lucas_lehmer.X.right_distrib {q : ℕ+} (x y z : lucas_lehmer.X q) :
(x + y) * z = x * z + y * z
source
@[protected, instance]
def lucas_lehmer.X.ring {q : ℕ+} :
ring (lucas_lehmer.X q)
Equations
source
@[protected, instance]
def lucas_lehmer.X.comm_ring {q : ℕ+} :
comm_ring (lucas_lehmer.X q)
Equations
source
@[protected, instance]
def lucas_lehmer.X.nontrivial {q : ℕ+} [fact (1 < ↑q)] :
nontrivial (lucas_lehmer.X q)
source
@[simp]
theorem lucas_lehmer.X.nat_coe_fst {q : ℕ+} (n : ℕ) :
↑n.fst = ↑n
source
@[simp]
theorem lucas_lehmer.X.nat_coe_snd {q : ℕ+} (n : ℕ) :
↑n.snd = 0
source
@[simp]
theorem lucas_lehmer.X.int_coe_fst {q : ℕ+} (n : ℤ) :
↑n.fst = ↑n
source
@[simp]
theorem lucas_lehmer.X.int_coe_snd {q : ℕ+} (n : ℤ) :
↑n.snd = 0
source
@[norm_cast]
theorem lucas_lehmer.X.coe_mul {q : ℕ+} (n m : ℤ) :
↑(n * m) = ↑n * ↑m
source
@[norm_cast]
theorem lucas_lehmer.X.coe_nat {q : ℕ+} (n : ℕ) :
↑↑n = ↑n
source
theorem lucas_lehmer.X.X_card {q : ℕ+} :
fintype.card (lucas_lehmer.X q) = ↑q ^ 2

The cardinality of X is q^2.

source
theorem lucas_lehmer.X.units_card {q : ℕ+} (w : 1 < q) :
fintype.card (lucas_lehmer.X q)ˣ < ↑q ^ 2

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

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

We define ω = 2 + √3.

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

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

Equations
source
theorem lucas_lehmer.X.ω_mul_ωb (q : ℕ+) :
lucas_lehmer.X.ω * lucas_lehmer.X.ωb = 1
source
theorem lucas_lehmer.X.ωb_mul_ω (q : ℕ+) :
lucas_lehmer.X.ωb * lucas_lehmer.X.ω = 1
source
theorem lucas_lehmer.X.closed_form {q : ℕ+} (i : ℕ) :
↑(lucas_lehmer.s i) = lucas_lehmer.X.ω ^ 2 ^ i + lucas_lehmer.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 ℕ.

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

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

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

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

ω as an element of the group of units.

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

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

source
theorem lucas_lehmer.order_ineq (p' : ℕ) (h : lucas_lehmer.lucas_lehmer_residue (p' + 2) = 0) :
2 ^ (p' + 2) < ↑(lucas_lehmer.q (p' + 2)) ^ 2
source
theorem lucas_lehmer_sufficiency (p : ℕ) (w : 1 < p) :
lucas_lehmer.lucas_lehmer_test p → nat.prime (mersenne p)
source
theorem lucas_lehmer.s_mod_succ {p : ℕ} {a : ℤ} {i : ℕ} {b c : ℤ} (h1 : 2 ^ p - 1 = a) (h2 : lucas_lehmer.s_mod p i = b) (h3 : (b * b - 2) % a = c) :
lucas_lehmer.s_mod p (i + 1) = c
source
meta def lucas_lehmer.run_test  :
tactic unit

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!

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