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

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.

source

def mersenne :

ℕ → ℕ

The Mersenne numbers, 2^p - 1.

Equations

mersenne p = 2 ^ p - 1

source

theorem mersenne_pos {p : ℕ} :

0 < p → 0 < mersenne p

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

lucas_lehmer.s (i + 1) = lucas_lehmer.s i ^ 2 - 2
lucas_lehmer.s 0 = 4

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

lucas_lehmer.s_zmod p (i + 1) = lucas_lehmer.s_zmod p i ^ 2 - 2
lucas_lehmer.s_zmod p 0 = 4

source

def lucas_lehmer.s_mod :

ℕ → ℕ → ℤ

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

Equations

lucas_lehmer.s_mod p (i + 1) = (lucas_lehmer.s_mod p i ^ 2 - 2) % (2 ^ p - 1)
lucas_lehmer.s_mod p 0 = 4 % (2 ^ p - 1)

source

theorem lucas_lehmer.mersenne_int_ne_zero (p : ℕ) :

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 : ℕ) :

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

lucas_lehmer.lucas_lehmer_residue p = lucas_lehmer.s_zmod p (p - 2)

source

theorem lucas_lehmer.residue_eq_zero_iff_s_mod_eq_zero (p : ℕ) :

1 < p → (lucas_lehmer.lucas_lehmer_residue p = 0 ↔ lucas_lehmer.s_mod p (p - 2) = 0)

source

def lucas_lehmer.lucas_lehmer_test :

ℕ → 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

lucas_lehmer.lucas_lehmer_test p = (lucas_lehmer.lucas_lehmer_residue p = 0)

source

@[instance]

def lucas_lehmer.lucas_lehmer_test.decidable_pred :

decidable_pred lucas_lehmer.lucas_lehmer_test

source

def lucas_lehmer.q :

ℕ → ℕ+

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

Equations

lucas_lehmer.q p = ⟨(mersenne p).min_fac, _⟩

source

@[instance]

def lucas_lehmer.fact_pnat_pos (q : ℕ+) :

fact (0 < ↑q)

Equations

_ = _

source

@[instance]

def lucas_lehmer.X.inhabited (q : ℕ+) :

inhabited (lucas_lehmer.X q)

source

def lucas_lehmer.X :

ℕ+ → Type

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

Equations

lucas_lehmer.X q = (zmod ↑q × zmod ↑q)

source

@[instance]

def lucas_lehmer.X.add_comm_group (q : ℕ+) :

add_comm_group (lucas_lehmer.X q)

source

@[instance]

def lucas_lehmer.X.decidable_eq (q : ℕ+) :

decidable_eq (lucas_lehmer.X q)

source

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

x.fst = y.fst → 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

@[instance]

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

has_mul (lucas_lehmer.X q)

Equations

lucas_lehmer.X.has_mul = {mul := λ (x y : lucas_lehmer.X q), (x.fst * y.fst + 3 * x.snd * y.snd, x.fst * y.snd + x.snd * y.fst)}

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

@[instance]

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

has_one (lucas_lehmer.X q)

Equations

lucas_lehmer.X.has_one = {one := (1, 0)}

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

@[instance]

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

monoid (lucas_lehmer.X q)

Equations

lucas_lehmer.X.monoid = {mul := has_mul.mul infer_instance, mul_assoc := _, one := (1, 0), one_mul := _, mul_one := _}

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

@[instance]

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

ring (lucas_lehmer.X q)

Equations

lucas_lehmer.X.ring = {add := add_comm_group.add infer_instance, add_assoc := _, zero := add_comm_group.zero infer_instance, zero_add := _, add_zero := _, neg := add_comm_group.neg infer_instance, add_left_neg := _, add_comm := _, mul := monoid.mul infer_instance, mul_assoc := _, one := monoid.one infer_instance, one_mul := _, mul_one := _, left_distrib := _, right_distrib := _}

source

@[instance]

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

comm_ring (lucas_lehmer.X q)

Equations

lucas_lehmer.X.comm_ring = {add := ring.add infer_instance, add_assoc := _, zero := ring.zero infer_instance, zero_add := _, add_zero := _, neg := ring.neg infer_instance, add_left_neg := _, add_comm := _, mul := ring.mul infer_instance, mul_assoc := _, one := ring.one infer_instance, one_mul := _, mul_one := _, left_distrib := _, right_distrib := _, mul_comm := _}

source

@[instance]

def lucas_lehmer.X.nontrivial {q : ℕ+} [fact (1 < ↑q)] :

nontrivial (lucas_lehmer.X q)

Equations

lucas_lehmer.X.nontrivial = _

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

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

↑(n * m) = ↑n * ↑m

source

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 : ℕ+} :

1 < q → fintype.card (units (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

lucas_lehmer.X.ω = (2, 1)

source

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

lucas_lehmer.X q

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

Equations

lucas_lehmer.X.ωb = (2, -1)

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' : ℕ) :

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' : ℕ) :

lucas_lehmer.lucas_lehmer_residue (p' + 2) = 0 → lucas_lehmer.X.ω ^ 2 ^ (p' + 1) = -1

source

theorem lucas_lehmer.ω_pow_eq_one (p' : ℕ) :

lucas_lehmer.lucas_lehmer_residue (p' + 2) = 0 → lucas_lehmer.X.ω ^ 2 ^ (p' + 2) = 1

source

def lucas_lehmer.ω_unit (p : ℕ) :

units (lucas_lehmer.X (lucas_lehmer.q p))

ω as an element of the group of units.

Equations

lucas_lehmer.ω_unit p = {val := lucas_lehmer.X.ω (lucas_lehmer.q p), inv := lucas_lehmer.X.ωb (lucas_lehmer.q p), val_inv := _, inv_val := _}

source

@[simp]

theorem lucas_lehmer.ω_unit_coe (p : ℕ) :

↑(lucas_lehmer.ω_unit p) = lucas_lehmer.X.ω

source

theorem lucas_lehmer.order_ω (p' : ℕ) :

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' : ℕ) :

lucas_lehmer.lucas_lehmer_residue (p' + 2) = 0 → 2 ^ (p' + 2) < ↑(lucas_lehmer.q (p' + 2)) ^ 2

source

theorem lucas_lehmer_sufficiency (p : ℕ) :

1 < p → lucas_lehmer.lucas_lehmer_test p → nat.prime (mersenne p)

source

@[instance]

meta def lucas_lehmer.nat_pexpr :

has_to_pexpr ℕ

source

@[instance]

meta def lucas_lehmer.int_pexpr :

has_to_pexpr ℤ

source

theorem lucas_lehmer.s_mod_succ {p : ℕ} {a : ℤ} {i : ℕ} {b c : ℤ} :

2 ^ p - 1 = a → lucas_lehmer.s_mod p i = b → (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]

mathlib documentation

number_theory.lucas_lehmer

The Lucas-Lehmer test for Mersenne primes.

TODO

History

General documentation

Additional documentation

Library

mathlib documentation

number_theory.​lucas_lehmer

The Lucas-Lehmer test for Mersenne primes.

TODO

History

General documentation

Additional documentation

Library

number_theory.lucas_lehmer