Classic theorems without their usual tools

import Mathlib.Data.Nat.Choose.Sum
import Mathlib.Data.Nat.Choose.Dvd

theorem Lemma.fermat_little_key (p a : ℕ) (hp : p.Prime) :
    ∃ t, (a + 1) ^ p = a ^ p + 1 + p * t := by
  obtain ⟨m, rfl⟩ : ∃ m, p = m + 1 := ⟨p - 1, by have := hp.two_le; omega⟩
  rw [add_pow, Finset.sum_range_succ, Finset.sum_range_succ']
  simp only [Nat.cast_id]
  have hd : (m + 1) ∣ ∑ k ∈ Finset.range m,
      a ^ (k + 1) * 1 ^ (m + 1 - (k + 1)) * (m + 1).choose (k + 1) := by
    apply Finset.dvd_sum
    intro k hk
    rw [Finset.mem_range] at hk
    exact Dvd.dvd.mul_left (hp.dvd_choose_self (by omega) (by omega)) _
  obtain ⟨t, ht⟩ := hd
  refine ⟨t, ?_⟩
  rw [ht]
  simp only [pow_zero, one_pow, mul_one, one_mul, Nat.choose_zero_right, Nat.choose_self,
    Nat.sub_self]
  generalize a ^ (m + 1) = x
  generalize (m + 1) * t = y
  omega

theorem Lemma.fermat_little :
    ∀ p a : ℕ, p.Prime → a ^ p % p = a % p := by
  intro p a hp
  induction a with
  | zero =>
    rw [zero_pow hp.ne_zero]
  | succ n ih =>
    obtain ⟨t, ht⟩ := Lemma.fermat_little_key p n hp
    rw [ht, Nat.add_mul_mod_self_left, Nat.add_mod, ih, ← Nat.add_mod]

import Mathlib.FieldTheory.Finite.Basic

theorem Lemma.two_squares :
    ∀ p : ℕ, p.Prime → p % 4 = 1 → ∃ a b : ℕ, a ^ 2 + b ^ 2 = p := by
  intro p hp h4
  haveI := Fact.mk hp
  obtain ⟨x, hx⟩ : IsSquare (-1 : ZMod p) := by
    have hne : (-1 : ZMod p) ≠ 0 := neg_ne_zero.mpr one_ne_zero
    have hc : ringChar (ZMod p) ≠ 2 := by
      rw [ZMod.ringChar_zmod_n]; omega
    rw [FiniteField.isSquare_iff hc hne, ZMod.card]
    obtain ⟨t, ht⟩ : ∃ t, p / 2 = 2 * t := ⟨p / 4, by omega⟩
    rw [ht, pow_mul]; simp
  obtain ⟨m, hm⟩ : ∃ m, m = Nat.sqrt p := ⟨_, rfl⟩
  have hm1 : m * m ≤ p := hm ▸ Nat.sqrt_le p
  have hm2 : p < (m + 1) * (m + 1) := hm ▸ Nat.lt_succ_sqrt p
  have hmp : m * m ≠ p := by
    intro h
    rcases hp.eq_one_or_self_of_dvd m (Dvd.intro m h) with h1 | h1
    · subst h1; exact hp.one_lt.ne (by omega)
    · subst h1
      have := Nat.eq_of_mul_eq_mul_left hp.pos (h.trans (mul_one m).symm)
      exact hp.one_lt.ne this.symm
  have hcard : Fintype.card (ZMod p) < Fintype.card (Fin (m + 1) × Fin (m + 1)) := by
    simp only [ZMod.card, Fintype.card_prod, Fintype.card_fin]; exact hm2
  obtain ⟨⟨u1, v1⟩, ⟨u2, v2⟩, hne, heq⟩ := Fintype.exists_ne_map_eq_of_card_lt
    (fun uv : Fin (m + 1) × Fin (m + 1) =>
      ((uv.1 : ℕ) : ZMod p) + x * ((uv.2 : ℕ) : ZMod p)) hcard
  simp only at heq
  have hu1 := u1.isLt
  have hu2 := u2.isLt
  have hv1 := v1.isLt
  have hv2 := v2.isLt
  obtain ⟨a, ha⟩ : ∃ a : ℤ, a = ((u1 : ℕ) : ℤ) - ((u2 : ℕ) : ℤ) := ⟨_, rfl⟩
  obtain ⟨b, hb⟩ : ∃ b : ℤ, b = ((v1 : ℕ) : ℤ) - ((v2 : ℕ) : ℤ) := ⟨_, rfl⟩
  have hdiv : (p : ℤ) ∣ a ^ 2 + b ^ 2 := by
    rw [← ZMod.intCast_zmod_eq_zero_iff_dvd, ha, hb]
    push_cast
    have h' : ((u1 : ℕ) : ZMod p) - ((u2 : ℕ) : ZMod p)
        = -(x * (((v1 : ℕ) : ZMod p) - ((v2 : ℕ) : ZMod p))) := by
      linear_combination heq
    rw [h']
    linear_combination (-((((v1 : ℕ) : ZMod p) - ((v2 : ℕ) : ZMod p)) ^ 2)) * hx
  have hpos : 0 < a ^ 2 + b ^ 2 := by
    by_contra hc
    have ha0 : a = 0 := by nlinarith [sq_nonneg a, sq_nonneg b]
    have hb0 : b = 0 := by nlinarith [sq_nonneg a, sq_nonneg b]
    apply hne
    have e1 : u1 = u2 := Fin.ext (by omega)
    have e2 : v1 = v2 := Fin.ext (by omega)
    rw [e1, e2]
  have hmZ : (m : ℤ) * m < p := by
    have := lt_of_le_of_ne hm1 hmp
    exact_mod_cast this
  have hsa : a ^ 2 ≤ (m : ℤ) * m := by
    have h1 : a ≤ m := by omega
    have h2 : -(m : ℤ) ≤ a := by omega
    nlinarith
  have hsb : b ^ 2 ≤ (m : ℤ) * m := by
    have h1 : b ≤ m := by omega
    have h2 : -(m : ℤ) ≤ b := by omega
    nlinarith
  obtain ⟨k, hk⟩ := hdiv
  have hpZ : (0 : ℤ) < p := by exact_mod_cast hp.pos
  have hk1 : k = 1 := by
    have : 0 < k := by
      by_contra hc
      nlinarith
    have : k < 2 := by
      by_contra hc
      nlinarith
    omega
  subst hk1
  refine ⟨a.natAbs, b.natAbs, ?_⟩
  have : ((a.natAbs ^ 2 + b.natAbs ^ 2 : ℕ) : ℤ) = p := by
    push_cast
    rw [sq_abs, sq_abs]
    linarith
  exact_mod_cast this

import Mathlib.Data.Real.Basic
import Mathlib.Algebra.BigOperators.Fin

lemma amgm_bern (c d : ℝ) (hc : 0 ≤ c) (hcd : 0 ≤ c + d) :
    ∀ m : ℕ, c ^ m * (c + m * d) ≤ c * (c + d) ^ m := by
  intro m
  induction m with
  | zero => simp
  | succ m ih =>
    have h1 : c ^ m * (c + m * d) * (c + d) ≤ c * (c + d) ^ m * (c + d) :=
      mul_le_mul_of_nonneg_right ih hcd
    have h2 : 0 ≤ c ^ m * ((m : ℝ) * (d * d)) :=
      mul_nonneg (pow_nonneg hc m) (mul_nonneg (Nat.cast_nonneg m) (mul_self_nonneg d))
    have e1 : c ^ m * (c + m * d) * (c + d)
        = c ^ (m + 1) * (c + ((m + 1 : ℕ) : ℝ) * d) + c ^ m * ((m : ℝ) * (d * d)) := by
      push_cast; grind
    have e2 : c * (c + d) ^ m * (c + d) = c * (c + d) ^ (m + 1) := by grind
    grind

lemma amgm_step (A x : ℝ) (hA : 0 ≤ A) (hx : 0 ≤ x) (n : ℕ) :
    A ^ n * x ≤ ((n * A + x) / (n + 1)) ^ (n + 1) := by
  have hn : (0 : ℝ) < n + 1 := by
    have := Nat.cast_nonneg (α := ℝ) n
    grind
  have hB : 0 ≤ (n * A + x) / (n + 1) :=
    div_nonneg (add_nonneg (mul_nonneg (Nat.cast_nonneg n) hA) hx) hn.le
  have hmul : ((n : ℝ) + 1) * ((n * A + x) / (n + 1)) = n * A + x := by
    rw [mul_div_assoc']
    exact mul_div_cancel_left₀ _ hn.ne'
  rcases hA.lt_or_eq with hA' | hA'
  · have h := amgm_bern A ((n * A + x) / (n + 1) - A) hA (by grind) (n + 1)
    have e : A + ((n + 1 : ℕ) : ℝ) * ((n * A + x) / (n + 1) - A) = x := by
      push_cast; grind
    rw [e, show A + ((n * A + x) / (n + 1) - A) = (n * A + x) / (n + 1) by grind] at h
    have : A * (A ^ n * x) ≤ A * ((n * A + x) / (n + 1)) ^ (n + 1) := by
      calc A * (A ^ n * x) = A ^ (n + 1) * x := by grind
        _ ≤ _ := h
    exact le_of_mul_le_mul_left this hA'
  · subst hA'
    rcases n with _ | n
    · simp
    · simp only [ne_eq, Nat.add_one_ne_zero, not_false_eq_true, zero_pow, zero_mul]
      exact pow_nonneg hB _

lemma amgm_sum_nonneg : ∀ (n : ℕ) (x : Fin n → ℝ), (∀ i, 0 ≤ x i) → 0 ≤ ∑ i, x i := by
  intro n
  induction n with
  | zero => intro x _; simp
  | succ n ih =>
    intro x hx
    rw [Fin.sum_univ_castSucc]
    exact add_nonneg (ih _ (fun i => hx _)) (hx _)

theorem Lemma.am_gm : ∀ (n : ℕ) (x : Fin n → ℝ), (∀ i, 0 ≤ x i) → ∏ i, x i ≤ ((∑ i, x i) / n) ^ n := by
  intro n
  induction n with
  | zero => intro x _; simp
  | succ n ih =>
    intro x hx
    rw [Fin.prod_univ_castSucc, Fin.sum_univ_castSucc]
    have hS0 : 0 ≤ ∑ i : Fin n, x i.castSucc := amgm_sum_nonneg n _ (fun i => hx _)
    have hih := ih (fun i => x i.castSucc) (fun i => hx _)
    have hS00 : n = 0 → ∑ i : Fin n, x i.castSucc = 0 := by rintro rfl; simp
    generalize ∑ i : Fin n, x i.castSucc = S at hS0 hih hS00 ⊢
    have hSA : S = n * (S / n) := by
      rcases Nat.eq_zero_or_pos n with h | h
      · simp [hS00 h, h]
      · have : (n : ℝ) ≠ 0 := by exact_mod_cast h.ne'
        rw [mul_div_assoc']
        exact (mul_div_cancel_left₀ _ this).symm
    have hA : 0 ≤ S / n := div_nonneg hS0 (Nat.cast_nonneg n)
    have key := amgm_step (S / n) (x (Fin.last n)) hA (hx _) n
    rw [← hSA] at key
    push_cast
    calc (∏ i : Fin n, x i.castSucc) * x (Fin.last n)
        ≤ (S / n) ^ n * x (Fin.last n) := mul_le_mul_of_nonneg_right hih (hx _)
      _ ≤ _ := key

import Mathlib.FieldTheory.Finite.Basic

open Finset

namespace QR

/-- the "absolute least residue" of a*i mod p -/
def m (p a i : ℕ) : ℕ := min ((a * i) % p) (p - (a * i) % p)

lemma r_pos (p a i : ℕ) (hp : p.Prime) (ha : ¬ p ∣ a) (h1 : 1 ≤ i) (h2 : i < p) :
    0 < (a * i) % p := by
  rcases Nat.eq_zero_or_pos ((a * i) % p) with h | h
  · exfalso
    have hd : p ∣ a * i := Nat.dvd_of_mod_eq_zero h
    rcases (Nat.Prime.dvd_mul hp).mp hd with h' | h'
    · exact ha h'
    · have := Nat.le_of_dvd (by omega) h'
      omega
  · exact h

lemma odd_p (p : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) : p % 2 = 1 := by
  rcases hp.eq_two_or_odd with h | h
  · exact absurd h hp2
  · exact h

-- L0 : bijection
theorem m_mem (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) :
    ∀ i ∈ Ico 1 (p / 2 + 1), m p a i ∈ Ico 1 (p / 2 + 1) := by
  intro i hi
  rw [Finset.mem_Ico] at hi
  have hodd := odd_p p hp hp2
  have hr := r_pos p a i hp ha hi.1 (by omega)
  have hlt := Nat.mod_lt (a * i) hp.pos
  rw [Finset.mem_Ico]
  unfold m
  constructor
  · apply le_min <;> omega
  · have := min_le_left ((a * i) % p) (p - (a * i) % p)
    have := min_le_right ((a * i) % p) (p - (a * i) % p)
    rcases le_total ((a * i) % p) (p - (a * i) % p) with h | h
    · rw [min_eq_left h]; omega
    · rw [min_eq_right h]; omega

theorem m_inj (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) :
    Set.InjOn (m p a) (Ico 1 (p / 2 + 1) : Set ℕ) := by
  intro i hi j hj h
  simp only [Finset.coe_Ico, Set.mem_Ico, Finset.mem_coe, Finset.mem_Ico] at hi hj
  have hodd := odd_p p hp hp2
  have hri := r_pos p a i hp ha hi.1 (by omega)
  have hrj := r_pos p a j hp ha hj.1 (by omega)
  have hlti := Nat.mod_lt (a * i) hp.pos
  have hltj := Nat.mod_lt (a * j) hp.pos
  unfold m at h
  have key : (a * i) % p = (a * j) % p ∨ (a * i) % p + (a * j) % p = p := by
    rcases le_total ((a * i) % p) (p - (a * i) % p) with h1 | h1 <;>
    rcases le_total ((a * j) % p) (p - (a * j) % p) with h2 | h2 <;>
    simp only [min_eq_left h1, min_eq_right h1, min_eq_left h2, min_eq_right h2] at h <;>
    omega
  haveI := Fact.mk hp
  rcases key with k | k
  · have e := (ZMod.natCast_eq_natCast_iff' (a * i) (a * j) p).mpr k
    push_cast at e
    have ha' : (a : ZMod p) ≠ 0 := by
      rwa [Ne, CharP.cast_eq_zero_iff (ZMod p) p]
    have e2 := mul_left_cancel₀ ha' e
    have e3 := (ZMod.natCast_eq_natCast_iff' i j p).mp e2
    rwa [Nat.mod_eq_of_lt (by omega), Nat.mod_eq_of_lt (by omega)] at e3
  · exfalso
    have hz : (a * (i + j)) % p = 0 := by
      rw [Nat.mul_add, Nat.add_mod, k, Nat.mod_self]
    have hd : p ∣ a * (i + j) := Nat.dvd_of_mod_eq_zero hz
    rcases (Nat.Prime.dvd_mul hp).mp hd with h' | h'
    · exact ha h'
    · have := Nat.le_of_dvd (by omega) h'
      omega

theorem m_image (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) :
    (Ico 1 (p / 2 + 1)).image (m p a) = Ico 1 (p / 2 + 1) := by
  apply Finset.eq_of_subset_of_card_le
  · intro x hx
    rw [Finset.mem_image] at hx
    obtain ⟨i, hi, rfl⟩ := hx
    exact m_mem p a hp hp2 ha i hi
  · rw [Finset.card_image_of_injOn (m_inj p a hp hp2 ha)]

lemma term (p a i : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) (hi : i ∈ Ico 1 (p / 2 + 1)) :
    ((a : ZMod p) * (i : ZMod p)) =
      (if p / 2 < (a * i) % p then (-1 : ZMod p) else 1) * ((m p a i : ℕ) : ZMod p) := by
  rw [Finset.mem_Ico] at hi
  have hodd := odd_p p hp hp2
  have hr := r_pos p a i hp ha hi.1 (by omega)
  have hlt := Nat.mod_lt (a * i) hp.pos
  have e : ((a : ZMod p) * (i : ZMod p)) = (((a * i) % p : ℕ) : ZMod p) := by
    rw [ZMod.natCast_mod]; push_cast; ring
  rw [e]
  unfold m
  split_ifs with h
  · rw [min_eq_right (by omega), Nat.cast_sub (by omega), ZMod.natCast_self]
    ring
  · rw [min_eq_left (by omega)]
    ring

-- L2a : Gauss lemma
theorem gauss (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) :
    (a : ZMod p) ^ (p / 2) =
      (-1) ^ ((Ico 1 (p / 2 + 1)).filter (fun i => p / 2 < (a * i) % p)).card := by
  haveI := Fact.mk hp
  set S := Ico 1 (p / 2 + 1) with hS
  have hcard : S.card = p / 2 := by simp [hS]
  have h1 : ∏ i ∈ S, ((a : ZMod p) * (i : ZMod p)) = (a : ZMod p) ^ (p / 2) * ∏ i ∈ S, (i : ZMod p) := by
    rw [Finset.prod_mul_distrib, Finset.prod_const, hcard]
  have h2 : ∏ i ∈ S, ((a : ZMod p) * (i : ZMod p)) =
      (-1) ^ (S.filter (fun i => p / 2 < (a * i) % p)).card * ∏ i ∈ S, (i : ZMod p) := by
    rw [Finset.prod_congr rfl (fun i hi => term p a i hp hp2 ha hi), Finset.prod_mul_distrib,
      Finset.prod_ite, Finset.prod_const, Finset.prod_const_one, mul_one]
    congr 1
    have := Finset.prod_image (s := S) (g := m p a) (f := fun j : ℕ => ((j : ℕ) : ZMod p))
      (m_inj p a hp hp2 ha)
    rw [← this, hS, m_image p a hp hp2 ha]
  have hne : ∏ i ∈ S, (i : ZMod p) ≠ 0 := by
    rw [Finset.prod_ne_zero_iff]
    intro i hi
    rw [hS, Finset.mem_Ico] at hi
    rw [Ne, CharP.cast_eq_zero_iff (ZMod p) p]
    intro hd
    have := Nat.le_of_dvd (by omega) hd
    have := odd_p p hp hp2
    omega
  exact mul_right_cancel₀ hne (h1.symm.trans h2)

-- L2b : Eisenstein parity
theorem eisenstein_parity (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a)
    (hodd : a % 2 = 1) :
    ((Ico 1 (p / 2 + 1)).filter (fun i => p / 2 < (a * i) % p)).card % 2 =
      (∑ i ∈ Ico 1 (p / 2 + 1), (a * i) / p) % 2 := by
  have hp1 : p % 2 = 1 := by
    rcases Nat.mod_two_eq_zero_or_one p with h | h
    · rcases hp.eq_one_or_self_of_dvd 2 (Nat.dvd_of_mod_eq_zero h) with h' | h' <;> omega
    · exact h
  have hp0 : 0 < p := hp.pos
  set s := Ico 1 (p / 2 + 1) with hs
  set L := s.filter (fun i => p / 2 < (a * i) % p) with hL
  -- sum of m equals sum of i
  have hM : ∑ i ∈ s, m p a i = ∑ i ∈ s, i := by
    have h1 := m_image p a hp hp2 ha
    have h2 := m_inj p a hp hp2 ha
    rw [hs]
    conv_rhs => rw [← h1]
    rw [Finset.sum_image (fun x hx y hy hxy => h2 (by simpa using hx) (by simpa using hy) hxy)]
  have hmL : ∀ i ∈ L, m p a i + (a * i) % p = p := by
    intro i hi
    have hi' := (Finset.mem_filter.mp hi).2
    have hr : (a * i) % p < p := Nat.mod_lt _ hp0
    unfold m
    rw [Nat.min_def]; split_ifs <;> omega
  have hmN : ∀ i ∈ s.filter (fun i => ¬ p / 2 < (a * i) % p), m p a i = (a * i) % p := by
    intro i hi
    have hi' := (Finset.mem_filter.mp hi).2
    have hr : (a * i) % p < p := Nat.mod_lt _ hp0
    unfold m
    rw [Nat.min_def]; split_ifs <;> omega
  have hsplitM := Finset.sum_filter_add_sum_filter_not s (fun i => p / 2 < (a * i) % p) (m p a)
  have hsplitR := Finset.sum_filter_add_sum_filter_not s (fun i => p / 2 < (a * i) % p)
    (fun i => (a * i) % p)
  have hNsum : ∑ i ∈ s.filter (fun i => ¬ p / 2 < (a * i) % p), m p a i
      = ∑ i ∈ s.filter (fun i => ¬ p / 2 < (a * i) % p), (a * i) % p :=
    Finset.sum_congr rfl hmN
  have hLsum : ∑ i ∈ L, m p a i + ∑ i ∈ L, (a * i) % p = L.card * p := by
    rw [← Finset.sum_add_distrib, Finset.sum_congr rfl hmL, Finset.sum_const, smul_eq_mul]
  have hdiv : a * ∑ i ∈ s, i = p * ∑ i ∈ s, (a * i) / p + ∑ i ∈ s, (a * i) % p := by
    rw [Finset.mul_sum, Finset.mul_sum, ← Finset.sum_add_distrib]
    exact Finset.sum_congr rfl (fun i _ => (Nat.div_add_mod (a * i) p).symm)
  rw [← hL] at hsplitM hsplitR
  generalize ∑ i ∈ s, i = S at *
  generalize ∑ i ∈ s, (a * i) / p = F at *
  generalize ∑ i ∈ s, (a * i) % p = R at *
  generalize ∑ i ∈ s, m p a i = M at *
  generalize ∑ i ∈ L, m p a i = ML at *
  generalize ∑ i ∈ L, (a * i) % p = RL at *
  generalize ∑ i ∈ s.filter (fun i => ¬ p / 2 < (a * i) % p), m p a i = MN at *
  generalize ∑ i ∈ s.filter (fun i => ¬ p / 2 < (a * i) % p), (a * i) % p = RN at *
  generalize L.card = c at *
  obtain ⟨a', rfl⟩ : ∃ a', a = 2 * a' + 1 := ⟨a / 2, by omega⟩
  obtain ⟨p', rfl⟩ : ∃ p', p = 2 * p' + 1 := ⟨p / 2, by omega⟩
  have e1 : (2 * a' + 1) * S = 2 * (a' * S) + S := by ring
  have e2 : c * (2 * p' + 1) = 2 * (c * p') + c := by ring
  have e3 : (2 * p' + 1) * F = 2 * (p' * F) + F := by ring
  rw [e1, e3] at hdiv
  rw [e2] at hLsum
  generalize a' * S = X at *
  generalize c * p' = Y at *
  generalize p' * F = Z at *
  omega

-- L1 : lattice count
theorem lattice_count (a b n : ℕ) (hb : 0 < b) (hn : a / b ≤ n) (hnd : ¬ b ∣ a) :
    ((Ico 1 (n + 1)).filter (fun j => b * j < a)).card = a / b := by
  have : (Ico 1 (n + 1)).filter (fun j => b * j < a) = Ico 1 (a / b + 1) := by
    ext j
    simp only [mem_filter, mem_Ico]
    have h1 : j ≤ a / b ↔ j * b ≤ a := Nat.le_div_iff_mul_le hb
    have h2 : b * j ≠ a := fun h => hnd ⟨j, h.symm⟩
    have h3 : j * b = b * j := Nat.mul_comm j b
    constructor
    · rintro ⟨⟨h4, h5⟩, h6⟩
      refine ⟨h4, ?_⟩
      have : j ≤ a / b := h1.mpr (by omega)
      omega
    · rintro ⟨h4, h5⟩
      have : j * b ≤ a := h1.mp (by omega)
      refine ⟨⟨h4, by omega⟩, by omega⟩
  rw [this]
  simp

theorem lattice_bound (p q i : ℕ) (hp : p % 2 = 1) (hq : q % 2 = 1) (hi : i ≤ p / 2) :
    q * i / p ≤ q / 2 := by
  have hp0 : 0 < p := by omega
  have := (Nat.div_lt_iff_lt_mul hp0 (x := q * i) (y := q / 2 + 1)).mpr
  suffices h : q * i < (q / 2 + 1) * p by have := this h; omega
  obtain ⟨a, rfl⟩ : ∃ a, p = 2 * a + 1 := ⟨p / 2, by omega⟩
  obtain ⟨b, rfl⟩ : ∃ b, q = 2 * b + 1 := ⟨q / 2, by omega⟩
  have ha : (2 * a + 1) / 2 = a := by omega
  have hb : (2 * b + 1) / 2 = b := by omega
  rw [ha] at hi
  rw [hb]
  have : (2 * b + 1) * i ≤ (2 * b + 1) * a := Nat.mul_le_mul_left _ hi
  nlinarith

theorem lattice_nd (p q i : ℕ) (hp : p.Prime) (hq : q.Prime) (hpq : p ≠ q) (hi1 : 1 ≤ i)
    (hi : i ≤ p / 2) : ¬ p ∣ q * i := by
  intro h
  rcases (Nat.Prime.dvd_mul hp).mp h with h | h
  · exact hpq ((Nat.prime_dvd_prime_iff_eq hp hq).mp h)
  · have := Nat.le_of_dvd (by omega) h
    omega

theorem lattice (p q : ℕ) (hp : p.Prime) (hq : q.Prime) (hpq : p ≠ q) (hp2 : p ≠ 2) (hq2 : q ≠ 2) :
    (∑ i ∈ Ico 1 (p / 2 + 1), (q * i) / p) + (∑ j ∈ Ico 1 (q / 2 + 1), (p * j) / q)
      = (p / 2) * (q / 2) := by
  have hpo : p % 2 = 1 := Nat.odd_iff.mp (hp.odd_of_ne_two hp2)
  have hqo : q % 2 = 1 := Nat.odd_iff.mp (hq.odd_of_ne_two hq2)
  have key : ∀ i ∈ Ico 1 (p / 2 + 1), ∀ j ∈ Ico 1 (q / 2 + 1), q * i ≠ p * j := by
    intro i hi j _ h
    rw [mem_Ico] at hi
    exact lattice_nd p q i hp hq hpq (by omega) (by omega) ⟨j, h⟩
  symm
  calc (p / 2) * (q / 2) = ∑ i ∈ Ico 1 (p / 2 + 1), ∑ j ∈ Ico 1 (q / 2 + 1), (1 : ℕ) := by simp
    _ = ∑ i ∈ Ico 1 (p / 2 + 1), ∑ j ∈ Ico 1 (q / 2 + 1),
          ((if p * j < q * i then 1 else 0) + (if q * i < p * j then 1 else 0)) := by
        refine sum_congr rfl fun i hi => sum_congr rfl fun j hj => ?_
        have := key i hi j hj
        split_ifs <;> omega
    _ = ∑ i ∈ Ico 1 (p / 2 + 1), ∑ j ∈ Ico 1 (q / 2 + 1), (if p * j < q * i then 1 else 0)
        + ∑ j ∈ Ico 1 (q / 2 + 1), ∑ i ∈ Ico 1 (p / 2 + 1), (if q * i < p * j then 1 else 0) := by
        rw [sum_comm (s := Ico 1 (q / 2 + 1))]
        simp only [sum_add_distrib]
    _ = _ := by
        congr 1
        · refine sum_congr rfl fun i hi => ?_
          rw [mem_Ico] at hi
          rw [← card_filter]
          exact lattice_count (q * i) p (q / 2) hp.pos (lattice_bound p q i hpo hqo (by omega))
            (lattice_nd p q i hp hq hpq (by omega) (by omega))
        · refine sum_congr rfl fun j hj => ?_
          rw [mem_Ico] at hj
          rw [← card_filter]
          exact lattice_count (p * j) q (p / 2) hq.pos (lattice_bound q p j hqo hpo (by omega))
            (lattice_nd q p j hq hp (Ne.symm hpq) (by omega) (by omega))

theorem euler (p a : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) (ha : ¬ p ∣ a) :
    (∃ x : ℕ, x ^ 2 % p = a % p) ↔ (a : ZMod p) ^ (p / 2) = 1 := by
  haveI := Fact.mk hp
  have hne : (a : ZMod p) ≠ 0 := by
    rwa [Ne, CharP.cast_eq_zero_iff (ZMod p) p a]
  have hc : ringChar (ZMod p) ≠ 2 := by
    rw [ZMod.ringChar_zmod_n]; exact hp2
  have hsq := FiniteField.isSquare_iff hc hne
  rw [ZMod.card] at hsq
  rw [← hsq]
  constructor
  · rintro ⟨x, hx⟩
    refine ⟨(x : ZMod p), ?_⟩
    have : ((x ^ 2 : ℕ) : ZMod p) = (a : ZMod p) := (ZMod.natCast_eq_natCast_iff' _ _ _).mpr hx
    rw [← this]; push_cast; ring
  · rintro ⟨r, hr⟩
    refine ⟨r.val, ?_⟩
    rw [← ZMod.natCast_eq_natCast_iff']
    push_cast
    rw [ZMod.natCast_zmod_val, hr]; ring

theorem neg_one_pow_iff (p n : ℕ) (hp : p.Prime) (hp2 : p ≠ 2) :
    ((-1 : ZMod p) ^ n = 1) ↔ n % 2 = 0 := by
  haveI : Fact (2 < p) := ⟨by have := hp.two_le; omega⟩
  rcases Nat.even_or_odd n with h | h
  · rw [h.neg_one_pow]; simp only [true_iff]; exact Nat.even_iff.mp h
  · rw [h.neg_one_pow]
    constructor
    · intro h1; exact absurd h1 (ZMod.neg_one_ne_one)
    · intro h2; have := Nat.odd_iff.mp h; omega

end QR

theorem Lemma.quadratic_reciprocity :
    ∀ p q : ℕ, p.Prime → q.Prime → p ≠ q → p ≠ 2 → q ≠ 2 → (((∃ x : ℕ, x ^ 2 % q = p % q) ↔ (∃ y : ℕ, y ^ 2 % p = q % p)) ↔ (p % 4 = 1 ∨ q % 4 = 1)) := by
  intro p q hp hq hpq hp2 hq2
  have hqp : ¬ q ∣ p := fun h => hpq ((Nat.prime_dvd_prime_iff_eq hq hp).mp h).symm
  have hpq' : ¬ p ∣ q := fun h => hpq ((Nat.prime_dvd_prime_iff_eq hp hq).mp h)
  have hp1 : p % 2 = 1 := by
    rcases hp.eq_two_or_odd with h | h
    · exact absurd h hp2
    · exact h
  have hq1 : q % 2 = 1 := by
    rcases hq.eq_two_or_odd with h | h
    · exact absurd h hq2
    · exact h
  rw [QR.euler q p hq hq2 hqp, QR.euler p q hp hp2 hpq', QR.gauss q p hq hq2 hqp,
    QR.gauss p q hp hp2 hpq', QR.neg_one_pow_iff q _ hq hq2, QR.neg_one_pow_iff p _ hp hp2,
    QR.eisenstein_parity q p hq hq2 hqp hp1, QR.eisenstein_parity p q hp hp2 hpq' hq1]
  have hl := QR.lattice p q hp hq hpq hp2 hq2
  have hm := Nat.mul_mod (p / 2) (q / 2) 2
  have hs := Nat.add_mod (∑ i ∈ Finset.Ico 1 (p / 2 + 1), q * i / p)
    (∑ j ∈ Finset.Ico 1 (q / 2 + 1), p * j / q) 2
  rw [hl] at hs
  rcases Nat.mod_two_eq_zero_or_one (p / 2) with h1 | h1 <;>
  rcases Nat.mod_two_eq_zero_or_one (q / 2) with h2 | h2 <;>
  simp only [h1, h2] at hm <;> omega