Primes in progressions, Euclid's way

import Mathlib.NumberTheory.LegendreSymbol.QuadraticReciprocity

namespace Lemma.P54

theorem z5 : (5 : ZMod 5) = 0 := by decide

theorem nsq2 : ∀ r : ZMod 5, ((2 : ℕ) : ZMod 5) ≠ r * r := by decide

theorem nsq3 : ∀ r : ZMod 5, ((3 : ℕ) : ZMod 5) ≠ r * r := by decide

theorem key (y : ℕ) (hy : 2 ∣ y) (hy0 : 0 < y) (q : ℕ) (hq : q.Prime)
    (hd : q ∣ 5 * y ^ 2 - 1) : (q % 5 = 1 ∨ q % 5 = 4) ∧ ¬ q ∣ y := by
  haveI := Fact.mk hq
  haveI : Fact (Nat.Prime 5) := Fact.mk Nat.prime_five
  have h1 : ((5 * y ^ 2 - 1 : ℕ) : ZMod q) = 0 :=
    (CharP.cast_eq_zero_iff (ZMod q) q _).2 hd
  have h5 : 1 ≤ 5 * y ^ 2 := by
    have : 1 ≤ y ^ 2 := Nat.one_le_pow _ _ hy0
    omega
  have h2 : (5 : ZMod q) * (y : ZMod q) ^ 2 = 1 := by
    rw [Nat.cast_sub h5] at h1
    push_cast at h1
    linear_combination h1
  have hndvd : ¬ q ∣ y := by
    intro hqy
    have : (y : ZMod q) = 0 := (CharP.cast_eq_zero_iff (ZMod q) q _).2 hqy
    rw [this] at h2
    simp at h2
  refine ⟨?_, hndvd⟩
  have hq2 : q ≠ 2 := by
    rintro rfl
    exact hndvd hy
  have hq5 : q ≠ 5 := by
    rintro rfl
    have : (5 : ZMod 5) = 0 := z5
    rw [this] at h2
    simp at h2
  have hsq : IsSquare ((5 : ℕ) : ZMod q) := ⟨5 * y, by push_cast; linear_combination (-5) * h2⟩
  obtain ⟨r, hr⟩ :=
    (ZMod.exists_sq_eq_prime_iff_of_mod_four_eq_one (p := 5) (q := q) (by norm_num) hq2).2 hsq
  rw [← ZMod.natCast_mod q 5] at hr
  have hlt := Nat.mod_lt q (by norm_num : 5 > 0)
  generalize hm : q % 5 = m at hr hlt
  obtain h | h | h | h | h : m = 0 ∨ m = 1 ∨ m = 2 ∨ m = 3 ∨ m = 4 := by omega
  · subst h
    exact absurd ((Nat.prime_dvd_prime_iff_eq Nat.prime_five hq).1
      (Nat.dvd_of_mod_eq_zero hm)).symm hq5
  · exact Or.inl h
  · subst h
    exact absurd hr (nsq2 r)
  · subst h
    exact absurd hr (nsq3 r)
  · exact Or.inr h

theorem exists4 (N : ℕ) (hN : N % 5 = 4)
    (hf : ∀ q, q.Prime → q ∣ N → q % 5 = 1 ∨ q % 5 = 4) :
    ∃ q, q.Prime ∧ q ∣ N ∧ q % 5 = 4 := by
  induction N using Nat.strong_induction_on with
  | _ N ih =>
    have hN1 : N ≠ 1 := by omega
    have hq := Nat.minFac_prime hN1
    have hqd := Nat.minFac_dvd N
    rcases hf _ hq hqd with h | h
    · obtain ⟨k, hk⟩ := hqd
      have hk5 : k % 5 = 4 := by
        rw [hk, Nat.mul_mod, h] at hN
        omega
      have hklt : k < N := by
        have := hq.two_le
        have hk0 : 0 < k := by omega
        rw [hk]
        nlinarith
      have hdiv : k ∣ N := Dvd.intro_left _ hk.symm
      obtain ⟨p, hp, hpd, hp4⟩ :=
        ih k hklt hk5 (fun p hp hpd => hf p hp (hpd.trans hdiv))
      exact ⟨p, hp, hpd.trans hdiv, hp4⟩
    · exact ⟨_, hq, hqd, h⟩

end Lemma.P54

theorem Lemma.primes_mod_5_eq_4 :
    Set.Infinite {p : ℕ | p.Prime ∧ p % 5 = 4} := by
  apply Set.infinite_of_forall_exists_gt
  intro n
  have hy2 : 2 ∣ (n + 5).factorial := Nat.dvd_factorial (by norm_num) (by omega)
  have hy0 : 0 < (n + 5).factorial := Nat.factorial_pos _
  have hN : (5 * (n + 5).factorial ^ 2 - 1) % 5 = 4 := by
    have : 1 ≤ (n + 5).factorial ^ 2 := Nat.one_le_pow _ _ hy0
    omega
  obtain ⟨q, hq, hqd, hq4⟩ := Lemma.P54.exists4 _ hN
    (fun q hq hd => (Lemma.P54.key _ hy2 hy0 q hq hd).1)
  refine ⟨q, ⟨hq, hq4⟩, ?_⟩
  by_contra hle
  push_neg at hle
  exact (Lemma.P54.key _ hy2 hy0 q hq hqd).2 (Nat.dvd_factorial hq.pos (by omega))

import Mathlib.NumberTheory.LegendreSymbol.QuadraticReciprocity

namespace Lemma.P61

theorem key (x : ℕ) (hx : 6 ∣ x) (q : ℕ) (hq : q.Prime)
    (hd : q ∣ x ^ 2 + x + 1) : q % 6 = 1 ∧ ¬ q ∣ x := by
  haveI := Fact.mk hq
  haveI : Fact (Nat.Prime 3) := Fact.mk Nat.prime_three
  have h1 : ((x ^ 2 + x + 1 : ℕ) : ZMod q) = 0 :=
    (CharP.cast_eq_zero_iff (ZMod q) q _).2 hd
  push_cast at h1
  have hndvd : ¬ q ∣ x := by
    intro hqx
    have : (x : ZMod q) = 0 := (CharP.cast_eq_zero_iff (ZMod q) q _).2 hqx
    rw [this] at h1
    simp at h1
  refine ⟨?_, hndvd⟩
  have hx0 : (x : ZMod q) ≠ 0 := by
    intro h
    exact hndvd ((CharP.cast_eq_zero_iff (ZMod q) q _).1 h)
  have hx1 : (x : ZMod q) ≠ 1 := by
    intro h
    rw [h] at h1
    have h3 : ((3 : ℕ) : ZMod q) = 0 := by push_cast; linear_combination h1
    have hq3 : q ∣ 3 := (CharP.cast_eq_zero_iff (ZMod q) q _).1 h3
    have : q = 3 := (Nat.prime_dvd_prime_iff_eq hq Nat.prime_three).1 hq3
    exact hndvd (this ▸ (Dvd.dvd.trans (by norm_num : 3 ∣ 6) hx))
  have hx3 : (x : ZMod q) ^ 3 = 1 := by
    linear_combination ((x : ZMod q) - 1) * h1
  have hord : orderOf (x : ZMod q) = 3 := orderOf_eq_prime hx3 hx1
  have hdvd : 3 ∣ q - 1 := hord ▸ ZMod.orderOf_dvd_card_sub_one hx0
  have hq2 : q ≠ 2 := by
    rintro rfl
    obtain ⟨k, rfl⟩ := hx
    apply absurd hd
    rw [Nat.dvd_iff_mod_eq_zero]
    have : (6 * k) ^ 2 + 6 * k + 1 = 2 * (18 * k ^ 2 + 3 * k) + 1 := by ring
    rw [this]
    omega
  have h2 := hq.two_le
  rcases Nat.even_or_odd q with he | ho
  · exact absurd ((Nat.prime_dvd_prime_iff_eq Nat.prime_two hq).1 (even_iff_two_dvd.1 he)).symm hq2
  · obtain ⟨k, hk⟩ := ho
    omega

end Lemma.P61

theorem Lemma.primes_mod_6_eq_1 :
    Set.Infinite {p : ℕ | p.Prime ∧ p % 6 = 1} := by
  apply Set.infinite_of_forall_exists_gt
  intro n
  set x := (n + 3).factorial with hxdef
  have hx6 : 6 ∣ x := by
    have : (3 : ℕ).factorial ∣ x := Nat.factorial_dvd_factorial (by omega)
    simpa [Nat.factorial] using this
  have hN : x ^ 2 + x + 1 ≠ 1 := by
    have : 0 < x := Nat.factorial_pos _
    omega
  have hq := Nat.minFac_prime hN
  have hqd := Nat.minFac_dvd (x ^ 2 + x + 1)
  obtain ⟨h6, hnd⟩ := Lemma.P61.key x hx6 _ hq hqd
  refine ⟨_, ⟨hq, h6⟩, ?_⟩
  by_contra hle
  push_neg at hle
  exact hnd (Nat.dvd_factorial hq.pos (by omega))

import Mathlib.NumberTheory.LegendreSymbol.QuadraticReciprocity

namespace Lemma.P65

theorem exists_factor (N : ℕ) (hN : N % 6 = 5) :
    ∃ q, q.Prime ∧ q % 6 = 5 ∧ q ∣ N := by
  induction N using Nat.strong_induction_on with
  | _ N ih =>
    have hN1 : N ≠ 1 := by omega
    have hq := Nat.minFac_prime hN1
    have hqd := Nat.minFac_dvd N
    set q := N.minFac with hqdef
    have hq2 : q % 2 = 1 := by
      by_contra h
      have h2 : 2 ∣ q := by omega
      have : q = 2 := ((Nat.prime_dvd_prime_iff_eq Nat.prime_two hq).1 h2).symm
      rw [this] at hqd
      omega
    have hq3 : q % 3 ≠ 0 := by
      intro h
      have h3 : 3 ∣ q := by omega
      have : q = 3 := ((Nat.prime_dvd_prime_iff_eq Nat.prime_three hq).1 h3).symm
      rw [this] at hqd
      omega
    by_cases h5 : q % 6 = 5
    · exact ⟨q, hq, h5, hqd⟩
    · have h1 : q % 6 = 1 := by omega
      obtain ⟨m, hm⟩ := hqd
      have hmod := Nat.mul_mod q m 6
      rw [h1, ← hm, hN] at hmod
      have hm5 : m % 6 = 5 := by omega
      have hq2' := hq.two_le
      have hmlt : m < N := by
        have hm0 : 0 < m := by
          apply Nat.pos_of_ne_zero
          rintro rfl
          rw [mul_zero] at hm
          omega
        rw [hm]
        calc m = 1 * m := (one_mul m).symm
          _ < q * m := Nat.mul_lt_mul_of_pos_right hq.one_lt hm0
      obtain ⟨r, hr, hr5, hrd⟩ := ih m hmlt hm5
      exact ⟨r, hr, hr5, hrd.trans ⟨q, by rw [hm, mul_comm]⟩⟩

end Lemma.P65

theorem Lemma.primes_mod_6_eq_5 :
    Set.Infinite {p : ℕ | p.Prime ∧ p % 6 = 5} := by
  apply Set.infinite_of_forall_exists_gt
  intro n
  set x := (n + 3).factorial with hxdef
  have hx0 : 0 < x := Nat.factorial_pos _
  obtain ⟨q, hq, h5, hd⟩ := Lemma.P65.exists_factor (6 * x - 1) (by omega)
  refine ⟨q, ⟨hq, h5⟩, ?_⟩
  by_contra hle
  push_neg at hle
  have hqx : q ∣ 6 * x := Dvd.dvd.mul_left (Nat.dvd_factorial hq.pos (by omega)) 6
  have h1 : q ∣ 6 * x - (6 * x - 1) := Nat.dvd_sub hqx hd
  have : 6 * x - (6 * x - 1) = 1 := by omega
  rw [this] at h1
  exact hq.one_lt.ne' (Nat.dvd_one.1 h1)