feat: port Archive.Examples.MersennePrimes (leanprover-community#5704)

Co-authored-by: thorimur <68410468+thorimur@users.noreply.github.com>
Co-authored-by: Kyle Miller <kmill31415@gmail.com>

Author: 3 people

Archive.lean

@@ -1,4 +1,5 @@
 import Archive.Arithcc
+import Archive.Examples.MersennePrimes
 import Archive.Examples.PropEncodable
 import Archive.Imo.Imo1959Q1
 import Archive.Imo.Imo1960Q1

Archive/Examples/MersennePrimes.lean

@@ -0,0 +1,92 @@
+/-
+Copyright (c) 2020 Scott Morrison. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Scott Morrison
+
+! This file was ported from Lean 3 source module examples.mersenne_primes
+! leanprover-community/mathlib commit 58581d0fe523063f5651df0619be2bf65012a94a
+! Please do not edit these lines, except to modify the commit id
+! if you have ported upstream changes.
+-/
+import Mathlib.NumberTheory.LucasLehmer
+
+/-!
+# Explicit Mersenne primes
+
+We run some Lucas-Lehmer tests to prove the first Mersenne primes are prime.
+
+See the discussion at the end of [Mathlib/NumberTheory/LucasLehmer.lean]
+for ideas about extending this to larger Mersenne primes.
+-/
+
+-- The Lucas-Lehmer test does not apply to `mersenne 2`
+example : ¬ LucasLehmerTest 2 := by norm_num
+
+example : (mersenne 2).Prime := by decide
+
+example : (mersenne 3).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 5).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 7).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 13).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 17).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 19).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+/-- 2147483647.Prime, Euler (1772) -/
+example : (mersenne 31).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+/-- Pervushin (1883), Seelhoff (1886) -/
+example : (mersenne 61).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 89).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 107).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+/-- Édouard Lucas (1876) -/
+example : (mersenne 127).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+/-- Lehmer and Robinson using SWAC computer, (1952) -/
+example : (mersenne 521).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 607).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 1279).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 2203).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 2281).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 3217).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 4253).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+example : (mersenne 4423).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+
+-- First failure ("deep recursion detected")
+/-
+example : (mersenne 9689).Prime :=
+  lucas_lehmer_sufficiency _ (by norm_num) (by norm_num)
+-/

Mathlib/NumberTheory/LucasLehmer.lean

@@ -21,10 +21,8 @@ import Mathlib.Tactic.IntervalCases
 We define `lucasLehmerResidue : Π p : ℕ, ZMod (2^p - 1)`, and
 prove `lucasLehmerResidue p = 0 → Prime (mersenne p)`.
 
-Porting note: the tactics have not been ported yet.
-
-We construct a tactic `LucasLehmer.run_test`, which iteratively certifies the arithmetic
-required to calculate the residue, and enables us to prove
+We construct a `norm_num` extension to calculate this residue to certify primality of Mersenne
+primes using `lucas_lehmer_sufficiency`.
 
 
 ## TODO
@@ -38,6 +36,8 @@ required to calculate the residue, and enables us to prove
 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.
 -/
 
 
@@ -170,8 +170,12 @@ def LucasLehmerTest (p : ℕ) : Prop :=
   lucasLehmerResidue p = 0
 #align lucas_lehmer.lucas_lehmer_test LucasLehmer.LucasLehmerTest
 
+-- Porting note: We have a fast `norm_num` extension, and we would rather use that than accidentally
+-- have `simp` use `decide`!
+/-
 instance : DecidablePred LucasLehmerTest :=
   inferInstanceAs (DecidablePred (lucasLehmerResidue · = 0))
+-/
 
 /-- `q` is defined as the minimum factor of `mersenne p`, bundled as an `ℕ+`. -/
 def q (p : ℕ) : ℕ+ :=
@@ -517,71 +521,105 @@ theorem lucas_lehmer_sufficiency (p : ℕ) (w : 1 < p) : LucasLehmerTest p → (
   exact not_lt_of_ge (Nat.sub_le _ _) h
 #align lucas_lehmer_sufficiency lucas_lehmer_sufficiency
 
--- Here we calculate the residue, very inefficiently, using `dec_trivial`. We can do much better.
-example : (mersenne 5).Prime :=
-  lucas_lehmer_sufficiency 5 (by norm_num) (by decide)
-
--- Next we use `norm_num` to calculate each `s p i`.
 namespace LucasLehmer
 
-open Tactic
+/-!
+### `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`.
+-/
 
-theorem sMod_succ {p a i b c} (h1 : (2 ^ p - 1 : ℤ) = a) (h2 : sMod p i = b)
-    (h3 : (b * b - 2) % a = c) : sMod p (i + 1) = c := by
-  substs a b c
-  rw [← sq]
+namespace norm_num_ext
+open Qq Lean Elab.Tactic Mathlib.Meta.NormNum
+
+/-- Version of `sMod` that is `ℕ`-valued. One should have `q = 2 ^ p - 1`.
+This can be reduced by the kernel. -/
+def sMod' (q : ℕ) : ℕ → ℕ
+  | 0 => 4 % q
+  | i + 1 => (sMod' q i ^ 2 + (q - 2)) % q
+
+theorem sMod'_eq_sMod (p k : ℕ) (hp : 2 ≤ p) : (sMod' (2 ^ p - 1) k : ℤ) = sMod p k := by
+  have h1 := calc
+    4 = 2 ^ 2 := by norm_num
+    _ ≤ 2 ^ p := Nat.pow_le_pow_of_le_right (by norm_num) hp
+  have h2 : 1 ≤ 2 ^ p := by linarith
+  induction k with
+  | zero =>
+    rw [sMod', sMod, Int.ofNat_emod]
+    simp [h2]
+  | succ k ih =>
+    rw [sMod', sMod, ← ih]
+    have h3 : 2 ≤ 2 ^ p - 1 := by
+      zify [h2]
+      calc
+        (2 : Int) ≤ 4 - 1 := by norm_num
+        _         ≤ 2 ^ p - 1 := by zify at h1; exact Int.sub_le_sub_right h1 _
+    zify [h2, h3]
+    rw [← add_sub_assoc, sub_eq_add_neg, add_assoc, add_comm _ (-2), ← add_assoc,
+      Int.add_emod_self, ← sub_eq_add_neg]
+
+lemma testTrueHelper (p : ℕ) (hp : Nat.blt 1 p = true) (h : sMod' (2 ^ p - 1) (p - 2) = 0) :
+    LucasLehmerTest p := by
+  rw [Nat.blt_eq] at hp
+  rw [LucasLehmerTest, LucasLehmer.residue_eq_zero_iff_sMod_eq_zero p hp, ← sMod'_eq_sMod p _ hp, h]
   rfl
-#align lucas_lehmer.s_mod_succ LucasLehmer.sMod_succ
-
--- /- ./././Mathport/Syntax/Translate/Expr.lean:330:4: warning: unsupported (TODO): `[tacs] -/
--- /- ./././Mathport/Syntax/Translate/Expr.lean:330:4: warning: unsupported (TODO): `[tacs] -/
--- /- ./././Mathport/Syntax/Translate/Expr.lean:330:4: warning: unsupported (TODO): `[tacs] -/
--- /-- Given a goal of the form `lucas_lehmer_test p`,
--- attempt to do the calculation using `norm_num` to certify each step.
--- -/
--- unsafe def run_test : tactic Unit := do
---   let q(LucasLehmerTest $(p)) ← target
---   sorry
---   sorry
---   let p ← eval_expr ℕ p
---   let-- Calculate the candidate Mersenne prime
---   M : ℤ := 2 ^ p - 1
---   let t ← to_expr ``(2 ^ $(q(p)) - 1 = $(q(M)))
---   let v ← to_expr ``((by norm_num : 2 ^ $(q(p)) - 1 = $(q(M))))
---   let w ← assertv `w t v
---   let t
---     ←-- base case
---         to_expr
---         ``(sMod $(q(p)) 0 = 4)
---   let v ← to_expr ``((by norm_num [LucasLehmer.sMod] : sMod $(q(p)) 0 = 4))
---   let h ← assertv `h t v
---   -- step case, repeated p-2 times
---       iterate_exactly
---       (p - 2) sorry
---   let h
---     ←-- now close the goal
---         get_local
---         `h
---   exact h
--- #align lucas_lehmer.run_test lucas_lehmer.run_test
 
-end LucasLehmer
+lemma testFalseHelper (p : ℕ) (hp : Nat.blt 1 p = true)
+    (h : Nat.ble 1 (sMod' (2 ^ p - 1) (p - 2))) : ¬ LucasLehmerTest p := by
+  rw [Nat.blt_eq] at hp
+  rw [Nat.ble_eq, Nat.succ_le, Nat.pos_iff_ne_zero] at h
+  rw [LucasLehmerTest, LucasLehmer.residue_eq_zero_iff_sMod_eq_zero p hp, ← sMod'_eq_sMod p _ hp]
+  simpa using h
+
+theorem isNat_lucasLehmerTest : {p np : ℕ} →
+    IsNat p np → LucasLehmerTest np → LucasLehmerTest p
+  | _, _, ⟨rfl⟩, h => h
+
+theorem isNat_not_lucasLehmerTest : {p np : ℕ} →
+    IsNat p np → ¬ LucasLehmerTest np → ¬ LucasLehmerTest p
+  | _, _, ⟨rfl⟩, h => h
+
+/-- Calculate `LucasLehmer.LucasLehmerTest p` for `2 ≤ p` by using kernel reduction for the
+`sMod'` function. -/
+@[norm_num LucasLehmer.LucasLehmerTest (_ : ℕ)]
+def evalLucasLehmerTest : NormNumExt where eval {u α} e := do
+  let .app _ (p : Q(ℕ)) ← Meta.whnfR e | failure
+  let sℕ : Q(AddMonoidWithOne ℕ) := q(instAddMonoidWithOneNat)
+  let ⟨ep, hp⟩ ← deriveNat p
+  let np := ep.natLit!
+  unless 1 < np do
+    failure
+  have h1ltp : Q(Nat.blt 1 $ep) := (q(Eq.refl true) : Expr)
+  if sMod' (2 ^ np - 1) (np - 2) = 0 then
+    have hs : Q(sMod' (2 ^ $ep - 1) ($ep - 2) = 0) := (q(Eq.refl 0) : Expr)
+    have pf : Q(LucasLehmerTest $ep) := q(testTrueHelper $ep $h1ltp $hs)
+    have pf' : Q(LucasLehmerTest $p) := q(isNat_lucasLehmerTest $hp $pf)
+    return .isTrue pf'
+  else
+    have hs : Q(Nat.ble 1 (sMod' (2 ^ $ep - 1) ($ep - 2)) = true) := (q(Eq.refl true) : Expr)
+    have pf : Q(¬ LucasLehmerTest $ep) := q(testFalseHelper $ep $h1ltp $hs)
+    have pf' : Q(¬ LucasLehmerTest $p) := q(isNat_not_lucasLehmerTest $hp $pf)
+    return .isFalse pf'
+
+end norm_num_ext
 
--- /-- We verify that the tactic works to prove `127.prime`. -/
--- example : (mersenne 7).Prime :=
---   lucas_lehmer_sufficiency _ (by norm_num)
---     (by
---       run_tac
---         lucas_lehmer.run_test)
+end LucasLehmer
 
 /-!
-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].
+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.
 
-`(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.
+(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