Research-Stack/0-Core-Formalism/lean/Semantics/Semantics/SidonSets.lean

import Mathlib.Data.Set.Basic
import Mathlib.Data.Finset.Basic
import Mathlib.Data.Finset.Sort
import Mathlib.Analysis.SpecialFunctions.Pow.Real
import Mathlib.Algebra.Order.Chebyshev
import Mathlib.Tactic.Zify
import Mathlib.FieldTheory.Finite.GaloisField
import Mathlib.FieldTheory.Finite.Trace
import Mathlib.FieldTheory.Finite.Basic
import Mathlib.FieldTheory.Minpoly.Field
import Mathlib.FieldTheory.IntermediateField.Basic
import Mathlib.FieldTheory.IntermediateField.Adjoin.Basic
import Mathlib.RingTheory.Trace.Basic
import Mathlib.RingTheory.PowerBasis
import Mathlib.LinearAlgebra.FiniteDimensional.Lemmas
import Mathlib.LinearAlgebra.Dimension.RankNullity
import Mathlib.LinearAlgebra.LinearIndependent.Defs
import Mathlib.LinearAlgebra.LinearIndependent.Lemmas
import Mathlib.LinearAlgebra.Span.Basic
import Mathlib.GroupTheory.SpecificGroups.Cyclic
import Mathlib.GroupTheory.QuotientGroup.Basic
import Mathlib.GroupTheory.Index
import Mathlib.GroupTheory.OrderOfElement
import Mathlib.Tactic.LinearCombination
import Mathlib.Tactic.FieldSimp
import Mathlib.NumberTheory.Bertrand
import Semantics.SidonSet

/-!
# Sidon Sets — Singer Construction Infrastructure

Port of the reusable Sidon-set infrastructure from Hulak–Ramos–de Queiroz (2026),
"Formalizing Singer Sidon Constructions and Sidon Set Infrastructure in Lean 4"
(arXiv: 2605.03274).

Original Lean 4 source: https://github.com/d0d1/singer-theorem-lean
Commit: 0c890589afc58e8955a5d7c3a609daff6447da31
License: GPL-3.0-only

This module ports the key reusable definitions and theorem statements from the
Erdos30 development into the Semantics namespace. The heavy algebraic proofs
(Singer construction, Lindström inequality, unconditional bounds) are left as
`sorry` with `TODO(lean-port)` markers, since the original code targets
Mathlib v4.29.0 while this project uses v4.30.0-rc2.

## Reusable components ported

1. **IsSidon** — Finset ℤ Sidon predicate (compatible with paper's Erdos30.Sidon)
2. **IsSidonMod** — Modular Sidon predicate
3. **IsIntervalSidon** — Interval Sidon predicate with containment
4. **IsSidonMaximum / sidonMaximum** — Extremal function h(N)
5. **Singer construction** — sidon set mod p²+p+1 of size p+1
6. **Lindström's cross-difference inequality** — (m-k)·k ≤ N-1
7. **h(N) = Θ(√N) bounds** — unconditional two-sided bounds
8. **Erdos30Statement** — formal Erdős Problem 30

## Relationship to existing Semantics.SidonSet

The existing `Semantics.SidonSet` uses a greedy `List Nat` generator with a
computable `isSidon : List Nat → Prop` check. This module provides the
mathematically rigorous `Finset ℤ` version used in the paper's proofs.
Both coexist: the List Nat version for computation, the Finset ℤ version
for formal combinatorics.

## References

- Singer, J. (1938). A theorem in finite projective geometry and some applications.
  *Trans. Amer. Math. Soc.*, 43, 377–385.
- Lindström, B. (1969). An inequality for B₂-sequences.
  *J. Combin. Theory*, 6(2), 211–212.
- Erdős, P. (1976). Problems and results in combinatorial number theory.
  *Astérisque*, 24–25, 295–310. (Problem 30)
-/

namespace Semantics.SidonSets

open Finset

/-! ## Core Sidon Definitions (Finset ℤ) -/

/-- The Sidon property for a finite set of integers: all pairwise sums a + b
    (with a, b ∈ A) are distinct up to reordering. This is the standard
    combinatorial definition used in the Erdős Problem 30 literature. -/
def IsSidon (A : Finset ℤ) : Prop :=
  ∀ ⦃a b c d : ℤ⦄,
    a ∈ A → b ∈ A → c ∈ A → d ∈ A →
    a + b = c + d →
      (a = c ∧ b = d) ∨ (a = d ∧ b = c)

/-- The Sidon property for a list of natural numbers (computable version).
    Compatible with `Semantics.SidonSet.isSidon`. -/
def IsSidonNat (s : List Nat) : Prop :=
  Semantics.SidonSet.isSidon s

/-! ## Modular Sidon Sets -/

/-- `IsSidonMod M A` means A is Sidon modulo M: for any a, b, c, d ∈ A,
    M ∣ ((a + b) - (c + d)) implies {a, b} = {c, d} as unordered pairs.
    This is the form needed for Singer's construction, which produces
    Sidon sets in Z/(q²+q+1)Z. -/
def IsSidonMod (M : ℤ) (A : Finset ℤ) : Prop :=
  ∀ ⦃a b c d : ℤ⦄,
    a ∈ A → b ∈ A → c ∈ A → d ∈ A →
    (M ∣ ((a + b) - (c + d))) →
    (a = c ∧ b = d) ∨ (a = d ∧ b = c)

/-- Modular Sidon implies integer Sidon. -/
theorem IsSidonMod.toIsSidon {M : ℤ} {A : Finset ℤ} (h : IsSidonMod M A) :
    IsSidon A := by
  intro a b c d ha hb hc hd hsum
  exact h ha hb hc hd (by rw [hsum, sub_self]; exact dvd_zero M)

/-! ## Interval Sidon Sets -/

/-- The interval {1, ..., N} as a Finset ℤ. Empty when N < 1. -/
noncomputable def interval (N : ℤ) : Finset ℤ := Finset.Icc 1 N

/-- A Sidon subset of {1, ..., N}. -/
structure IsIntervalSidon (N : ℤ) (A : Finset ℤ) : Prop where
  subset : ∀ x ∈ A, 1 ≤ x ∧ x ≤ N
  sidon : IsSidon A

/-- Enlarging the ambient interval preserves IsIntervalSidon. -/
theorem IsIntervalSidon.mono {A : Finset ℤ} {N M : ℤ}
    (h : IsIntervalSidon N A) (hle : N ≤ M) : IsIntervalSidon M A where
  subset x hx := ⟨(h.subset x hx).1, le_trans (h.subset x hx).2 hle⟩
  sidon := h.sidon

/-! ## Translation -/

/-- Translate a finset by t. -/
def translate (A : Finset ℤ) (t : ℤ) : Finset ℤ :=
  A.map (⟨fun x => x + t, fun _ _ h => add_right_cancel h⟩ : ℤ ↪ ℤ)

@[simp] theorem card_translate (A : Finset ℤ) (t : ℤ) :
    (translate A t).card = A.card := by
  simp [translate]

/-- Translation preserves the Sidon property. -/
theorem IsSidon.translate {A : Finset ℤ} (hA : IsSidon A) (t : ℤ) :
    IsSidon (translate A t) := by
  intro a b c d ha hb hc hd hsum
  rcases Finset.mem_map.1 ha with ⟨a', ha', ha_eq⟩
  rcases Finset.mem_map.1 hb with ⟨b', hb', hb_eq⟩
  rcases Finset.mem_map.1 hc with ⟨c', hc', hc_eq⟩
  rcases Finset.mem_map.1 hd with ⟨d', hd', hd_eq⟩
  have ha_val : a' + t = a := by simpa using ha_eq
  have hb_val : b' + t = b := by simpa using hb_eq
  have hc_val : c' + t = c := by simpa using hc_eq
  have hd_val : d' + t = d := by simpa using hd_eq
  have hsum' : a' + b' = c' + d' := by
    calc
      a' + b' = (a' + t) + (b' + t) - (t + t) := by ring
      _ = a + b - (t + t) := by simp [ha_val, hb_val]
      _ = c + d - (t + t) := by rw [hsum]
      _ = (c' + t) + (d' + t) - (t + t) := by simp [hc_val, hd_val]
      _ = c' + d' := by ring
  rcases hA ha' hb' hc' hd' hsum' with (⟨hac, hbd⟩ | ⟨had, hbc⟩)
  · left
    constructor
    · calc
        a = a' + t := ha_val.symm
        _ = c' + t := by rw [hac]
        _ = c := hc_val
    · calc
        b = b' + t := hb_val.symm
        _ = d' + t := by rw [hbd]
        _ = d := hd_val
  · right
    constructor
    · calc
        a = a' + t := ha_val.symm
        _ = d' + t := by rw [had]
        _ = d := hd_val
    · calc
        b = b' + t := hb_val.symm
        _ = c' + t := by rw [hbc]
        _ = c := hc_val

/-! ## Extremal Function h(N) -/

/-- `IsSidonMaximum N h` states that h is the maximum cardinality of an
    interval Sidon subset of {1, ..., N}. -/
def IsSidonMaximum (N h : ℕ) : Prop :=
  (∃ A : Finset ℤ, IsIntervalSidon (N : ℤ) A ∧ A.card = h) ∧
    ∀ {A : Finset ℤ}, IsIntervalSidon (N : ℤ) A → A.card ≤ h

/-- Helper: the maximum Sidon cardinality exists for every N. -/
private theorem sidonMaximum_exists (N : ℕ) : ∃ h, IsSidonMaximum N h := by
  let Z := Finset.Icc 1 (N : ℤ)
  let cards : Set ℕ := {k | ∃ (A : Finset ℤ), A ⊆ Z ∧ IsSidon A ∧ A.card = k}
  have h_nonempty : cards.Nonempty := by
    refine ⟨0, ∅, Finset.empty_subset Z, ?_, Finset.card_empty⟩
    intro a b c d ha hb hc hd hsum
    simp at ha
  have h_fin : Set.Finite cards := by
    have h_fin_img : Set.Finite ((Z.powerset.image Finset.card : Finset ℕ) : Set ℕ) :=
      Finset.finite_toSet _
    apply Set.Finite.subset h_fin_img
    intro k hk
    rcases hk with ⟨A, hA_sub, hA_sidon, hcard⟩
    refine Finset.mem_image.mpr ⟨A, ?_, hcard⟩
    exact Finset.mem_powerset.mpr hA_sub
  have h_finset_nonempty : h_fin.toFinset.Nonempty := by
    rcases h_nonempty with ⟨k, hk⟩
    refine ⟨k, h_fin.mem_toFinset.mpr hk⟩
  let m := h_fin.toFinset.max' h_finset_nonempty
  have hm_cards : m ∈ cards :=
    h_fin.mem_toFinset.mp (Finset.max'_mem _ h_finset_nonempty)
  rcases hm_cards with ⟨A, hA_sub, hA_sidon, hcard⟩
  refine ⟨m, ?_⟩
  constructor
  · refine ⟨A, ?_, hcard⟩
    refine { subset := λ x hx => ?_, sidon := hA_sidon }
    have hx_mem_icc : x ∈ Z := hA_sub hx
    rcases Finset.mem_Icc.1 hx_mem_icc with ⟨hx1, hx2⟩
    exact ⟨hx1, hx2⟩
  · intro B hB
    have hB_sub : B ⊆ Z := by
      intro x hx; rcases hB.subset x hx with ⟨hx1, hx2⟩
      exact Finset.mem_Icc.mpr ⟨hx1, hx2⟩
    have hB_card : B.card ∈ cards := ⟨B, hB_sub, hB.sidon, rfl⟩
    have hB_fin : B.card ∈ h_fin.toFinset := h_fin.mem_toFinset.mpr hB_card
    exact Finset.le_max' h_fin.toFinset (B.card) hB_fin

/-- The extremal Sidon function h(N) = max{|A| : A ⊆ {1,...,N} is Sidon}. -/
noncomputable def sidonMaximum (N : ℕ) : ℕ :=
  Classical.choose (sidonMaximum_exists N)

/-- The maximum exists for every N. -/
theorem sidonMaximum_isSidonMaximum (N : ℕ) :
    IsSidonMaximum N (sidonMaximum N) :=
  Classical.choose_spec (sidonMaximum_exists N)

/-- The maximum cardinality is unique. -/
theorem isSidonMaximum_unique {N h k : ℕ}
    (hh : IsSidonMaximum N h) (hk : IsSidonMaximum N k) :
    h = k := by
  rcases hh.1 with ⟨A, hA, hAcard⟩
  rcases hk.1 with ⟨B, hB, hBcard⟩
  have hle : h ≤ k := by rw [← hAcard]; exact hk.2 hA
  have hge : k ≤ h := by rw [← hBcard]; exact hh.2 hB
  omega

/-! ## Difference-Counting Upper Bound -/

/-- First upper bound: for any interval Sidon set A ⊆ {1,...,N},
    |A| ≤ √(2N) + 1. This follows from pair-difference counting. -/
theorem IsIntervalSidon.card_le {A : Finset ℤ} {N : ℕ}
    (h : IsIntervalSidon (N : ℤ) A) (hN : 1 ≤ N) :
    A.card ≤ Nat.sqrt (2 * N) + 1 := by
  set m := A.card with hm
  have hA_sidon : IsSidon A := h.sidon
  have hbound : ∀ x ∈ A, 1 ≤ x ∧ x ≤ N := h.subset
  -- All positive differences a-b (a,b ∈ A, a > b) are distinct by the Sidon property
  have h_diff_unique : ∀ a b c d : ℤ, a ∈ A → b ∈ A → c ∈ A → d ∈ A → a > b → c > d →
      a - b = c - d → a = c ∧ b = d := by
    intro a b c d ha hb hc hd ha_gt hc_gt h_eq
    have hsum : a + d = b + c := by omega
    rcases hA_sidon ha hd hb hc hsum with (⟨h1, h2⟩ | ⟨h1, h2⟩)
    · omega
    · exact ⟨h1, h2.symm⟩
  -- P = ordered pairs (a,b) ∈ A×A with a > b
  let P := (A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 > ab.2)
  have hP_injOn : Set.InjOn (λ (ab : ℤ × ℤ) => ab.1 - ab.2) (P : Set (ℤ × ℤ)) := by
    intro x hx y hy h
    rcases x with ⟨a,b⟩; rcases y with ⟨c,d⟩
    have haA : a ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hx).1)).1
    have hbA : b ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hx).1)).2
    have hcA : c ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hy).1)).1
    have hdA : d ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hy).1)).2
    have ha_gt_b : a > b := (Finset.mem_filter.1 hx).2
    have hc_gt_d : c > d := (Finset.mem_filter.1 hy).2
    rcases h_diff_unique a b c d haA hbA hcA hdA ha_gt_b hc_gt_d h with ⟨hac, hbd⟩
    ext <;> assumption
  have hP_card_diffs : (Finset.image (λ (ab : ℤ × ℤ) => ab.1 - ab.2) P).card = P.card :=
    Finset.card_image_of_injOn hP_injOn
  have hP_card : P.card = m * (m - 1) / 2 := by
    have h_total : (A.product A).card = m * m := by simp [hm, Finset.card_product]
    have h_swap_card : (Finset.image Prod.swap P).card = P.card :=
      Finset.card_image_of_injective _ Prod.swap_injective
    have h_swap_eq : Finset.image Prod.swap P = ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)) := by
      ext ⟨a, b⟩
      constructor
      · intro hmem
        rcases Finset.mem_image.1 hmem with ⟨⟨x, y⟩, hxy, hswap⟩
        rcases Finset.mem_filter.1 hxy with ⟨hprod, hgt⟩
        rcases Finset.mem_product.1 hprod with ⟨hx, hy⟩
        have h1 : y = a := congrArg Prod.fst hswap
        have h2 : x = b := congrArg Prod.snd hswap
        subst h1; subst h2
        exact Finset.mem_filter.2 ⟨Finset.mem_product.2 ⟨hy, hx⟩, hgt⟩
      · intro hmem
        rcases Finset.mem_filter.1 hmem with ⟨hprod, hlt⟩
        rcases Finset.mem_product.1 hprod with ⟨ha, hb⟩
        exact Finset.mem_image.2 ⟨(b, a),
          Finset.mem_filter.2 ⟨Finset.mem_product.2 ⟨hb, ha⟩, hlt⟩, rfl⟩
    have h_diag_card : ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)).card = m := by
      have himg : ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)) =
          A.image (λ (x : ℤ) => (x, x)) := by
        ext ⟨a, b⟩
        constructor
        · intro hmem
          rcases Finset.mem_filter.1 hmem with ⟨hprod, heq⟩
          have ha : a ∈ A := (Finset.mem_product.1 hprod).1
          have heq' : a = b := heq
          subst heq'
          exact Finset.mem_image.2 ⟨a, ha, rfl⟩
        · intro hmem
          rcases Finset.mem_image.1 hmem with ⟨x, hx, hxy⟩
          have h1 : x = a := congrArg Prod.fst hxy
          have h2 : x = b := congrArg Prod.snd hxy
          subst h1; subst h2
          exact Finset.mem_filter.2 ⟨Finset.mem_product.2 ⟨hx, hx⟩, rfl⟩
      rw [himg, Finset.card_image_of_injective _ (fun x y hxy => congrArg Prod.fst hxy)]
    -- Partition product into >, <, =
    have h_partition : (A.product A) = P ∪ ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)) ∪
        ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)) := by
      ext ⟨a, b⟩
      constructor
      · intro hmem
        rcases lt_trichotomy a b with hlt | heq | hgt
        · exact Finset.mem_union.2 (Or.inl (Finset.mem_union.2 (Or.inr
            (Finset.mem_filter.2 ⟨hmem, hlt⟩))))
        · exact Finset.mem_union.2 (Or.inr (Finset.mem_filter.2 ⟨hmem, heq⟩))
        · exact Finset.mem_union.2 (Or.inl (Finset.mem_union.2 (Or.inl
            (Finset.mem_filter.2 ⟨hmem, hgt⟩))))
      · intro hmem
        rcases Finset.mem_union.1 hmem with h | h
        · rcases Finset.mem_union.1 h with h' | h'
          · exact (Finset.mem_filter.1 h').1
          · exact (Finset.mem_filter.1 h').1
        · exact (Finset.mem_filter.1 h).1
    have h_disjoint_gt_lt : Disjoint P ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)) := by
      rw [Finset.disjoint_left]
      rintro ⟨a, b⟩ hP hlt
      have h1 : a > b := (Finset.mem_filter.1 hP).2
      have h2 : a < b := (Finset.mem_filter.1 hlt).2
      omega
    have h_disjoint_union_eq : Disjoint (P ∪ ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)))
        ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)) := by
      rw [Finset.disjoint_left]
      rintro ⟨a, b⟩ hmem heq
      have h2 : a = b := (Finset.mem_filter.1 heq).2
      rcases Finset.mem_union.1 hmem with h | h
      · have h1 : a > b := (Finset.mem_filter.1 h).2
        omega
      · have h1 : a < b := (Finset.mem_filter.1 h).2
        omega
    -- Count: |P| + |<| + |=| = m*m, and |P| = |<|
    have h_lt_card : ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)).card = P.card := by
      calc
        ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)).card =
            (Finset.image Prod.swap P).card := by rw [h_swap_eq]
        _ = P.card := h_swap_card
    have h_total_card : P.card + ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)).card +
        ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)).card = m * m := by
      calc
        P.card + ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)).card +
          ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2)).card =
          (P ∪ ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 < ab.2)) ∪
            ((A.product A).filter (λ (ab : ℤ × ℤ) => ab.1 = ab.2))).card := by
          rw [Finset.card_union_of_disjoint h_disjoint_union_eq,
            Finset.card_union_of_disjoint h_disjoint_gt_lt]
        _ = (A.product A).card := by conv_lhs => rw [← h_partition]
        _ = m * m := h_total
    rw [h_lt_card, h_diag_card] at h_total_card
    have hmm : m * (m - 1) + m = m * m := by
      rcases Nat.eq_zero_or_pos m with hm0 | hm0
      · simp [hm0]
      · calc m * (m - 1) + m = m * (m - 1 + 1) := by ring
          _ = m * m := by rw [Nat.sub_add_cancel hm0]
    omega
  have hD_bound : Finset.image (λ (ab : ℤ × ℤ) => ab.1 - ab.2) P ⊆ Finset.Icc 1 (N - 1 : ℤ) := by
    intro d hd
    rcases Finset.mem_image.1 hd with ⟨⟨a, b⟩, hP, hd_eq⟩
    have haA : a ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hP).1)).1
    have hbA : b ∈ A := (Finset.mem_product.1 ((Finset.mem_filter.1 hP).1)).2
    have ha_gt_b : a > b := (Finset.mem_filter.1 hP).2
    rcases hbound a haA with ⟨ha1, haN⟩
    rcases hbound b hbA with ⟨hb1, hbN⟩
    rw [← hd_eq]
    have h_pos : 1 ≤ a - b := by omega
    have h_le : a - b ≤ (N : ℤ) - 1 := by omega
    exact Finset.mem_Icc.mpr ⟨h_pos, h_le⟩
  have h_icc_card : (Finset.Icc 1 (N - 1 : ℤ) : Finset ℤ).card = N - 1 := by
    simp
  have hP_card_le : m * (m - 1) / 2 ≤ N - 1 := by
    calc
      m * (m - 1) / 2 = P.card := hP_card.symm
      _ = (Finset.image (λ (ab : ℤ × ℤ) => ab.1 - ab.2) P).card := hP_card_diffs.symm
      _ ≤ (Finset.Icc 1 (N - 1 : ℤ) : Finset ℤ).card := Finset.card_le_card hD_bound
      _ = N - 1 := h_icc_card
  -- From m*(m-1)/2 ≤ N-1, prove m ≤ √(2N) + 1 by contradiction
  have hpar : 2 ∣ m * (m - 1) := by
    rcases Nat.even_or_odd m with he | ho
    · exact he.two_dvd.mul_right _
    · exact ((Nat.Odd.sub_odd ho odd_one).two_dvd).mul_left _
  have hm_sq_sub_m_le : m * (m - 1) ≤ 2 * (N - 1) := by omega
  by_contra! H
  have hm_gt : m > Nat.sqrt (2 * N) + 1 := H
  set s := Nat.sqrt (2 * N) with hs
  have hm_ge : m ≥ s + 2 := by omega
  have hm_sq_sub_m_gt : m * (m - 1) > 2 * (N - 1) := by
    have h_sq_lt : 2 * N < (s + 1) * (s + 1) := Nat.lt_succ_sqrt (2 * N)
    have hm_m_hm1_ge : m * (m - 1) ≥ (s + 2) * (s + 1) :=
      Nat.mul_le_mul hm_ge (by omega)
    have h_gt : (s + 2) * (s + 1) > 2 * (N - 1) := by
      have hexp : (s + 2) * (s + 1) = (s + 1) * (s + 1) + (s + 1) := by ring
      omega
    omega
  omega

/-- The quadratic upper bound on sidonMaximum: h(N) ≤ √(2N) + 1. -/
theorem sidonMaximum_le_sqrt_two (N : ℕ) (hN : 1 ≤ N) :
    sidonMaximum N ≤ Nat.sqrt (2 * N) + 1 := by
  have hmax := sidonMaximum_isSidonMaximum N
  rcases hmax.1 with ⟨A, hA, hAcard⟩
  have hcard := hA.card_le hN
  rw [hAcard] at hcard
  exact hcard

/-! ## Lindström's Cross-Difference Inequality -/

/-- In a strictly sorted list, elements from `take k` are strictly less than
    elements from `drop k`. -/
private lemma sortedLT_take_lt_drop (l : List ℤ) (hs : l.SortedLT)
    {k : ℕ} (hk_lt : k < l.length) :
    ∀ a ∈ l.take k, ∀ b ∈ l.drop k, a < b := by
  intro a ha b hb
  rw [List.mem_take_iff_getElem] at ha
  rw [List.mem_drop_iff_getElem] at hb
  obtain ⟨i, hi, rfl⟩ := ha
  obtain ⟨j, hj, rfl⟩ := hb
  have hi' : i < l.length := by omega
  have hkj : k + j < l.length := by omega
  exact hs (show (⟨i, hi'⟩ : Fin l.length) < ⟨k + j, hkj⟩ from by
    simp [Fin.lt_def]; omega)

/-- In a Sidon set, cross-differences between disjoint subsets are distinct.
    If `L, R ⊆ A` are disjoint and `b - a = d - c` with `a, c ∈ L`, `b, d ∈ R`,
    then `a = c` and `b = d`. -/
theorem IsSidon.cross_diff_eq {A : Finset ℤ}
    (hA : IsSidon A) {L R : Finset ℤ}
    (hL : L ⊆ A) (hR : R ⊆ A) (hLR : Disjoint L R)
    {a b c d : ℤ} (ha : a ∈ L) (hb : b ∈ R)
    (hc : c ∈ L) (hd : d ∈ R)
    (heq : b - a = d - c) :
    a = c ∧ b = d := by
  have hsum : b + c = d + a := by linarith
  rcases hA (hR hb) (hL hc) (hR hd) (hL ha) hsum with h | h
  · exact ⟨h.2.symm, h.1⟩
  · exfalso
    have hba : b = a := h.1
    rw [hba] at hb
    exact Finset.disjoint_left.mp hLR ha hb

/-- If `L` and `R` are disjoint subsets of an interval Sidon set in `{1,...,N}`,
    and every element of `L` is strictly less than every element of `R`, then
    the cross-difference map `(a, b) ↦ b - a` is injective from `L × R` into
    `{1, ..., N-1}`, giving `|L| · |R| ≤ N - 1`. -/
theorem IsIntervalSidon.ordered_cross_diff_le {A : Finset ℤ} {N : ℤ}
    (hIS : IsIntervalSidon N A) {L R : Finset ℤ}
    (hL : L ⊆ A) (hR : R ⊆ A) (hLR : Disjoint L R)
    (hord : ∀ a ∈ L, ∀ b ∈ R, a < b)
    (hN : 1 ≤ N) (_hLne : L.Nonempty) (_hRne : R.Nonempty) :
    (L.card : ℤ) * R.card ≤ N - 1 := by
  let f : L × R → ℤ := fun ⟨⟨a, _⟩, ⟨b, _⟩⟩ => b - a
  have hf_inj : Function.Injective f := by
    intro ⟨⟨a, ha⟩, ⟨b, hb⟩⟩ ⟨⟨c, hc⟩, ⟨d, hd⟩⟩ heq
    simp only [f] at heq
    have := hIS.sidon.cross_diff_eq hL hR hLR ha hb hc hd heq
    simp [this.1, this.2]
  have hf_pos : ∀ x : L × R, 1 ≤ f x := by
    intro ⟨⟨a, ha⟩, ⟨b, hb⟩⟩
    simp only [f]
    linarith [hord a ha b hb]
  have hf_le : ∀ x : L × R, f x ≤ N - 1 := by
    intro ⟨⟨a, ha⟩, ⟨b, hb⟩⟩
    simp only [f]
    have ha_bound := hIS.subset a (hL ha)
    have hb_bound := hIS.subset b (hR hb)
    linarith
  have himage_sub : Finset.univ.image f ⊆ Finset.Icc 1 (N - 1) := by
    intro z hz
    rcases Finset.mem_image.mp hz with ⟨x, _, rfl⟩
    exact Finset.mem_Icc.mpr ⟨hf_pos x, hf_le x⟩
  have hcard : Fintype.card (L × R) = L.card * R.card := by
    simp [Fintype.card_prod, Fintype.card_coe]
  have himage_card : (Finset.univ.image f).card = L.card * R.card := by
    rw [Finset.card_image_of_injective _ hf_inj]
    simp [Fintype.card_prod, Fintype.card_coe]
  have hprod_le : L.card * R.card ≤ (N - 1).toNat := by
    calc L.card * R.card
        = (Finset.univ.image f).card := himage_card.symm
      _ ≤ (Finset.Icc 1 (N - 1)).card := Finset.card_le_card himage_sub
      _ ≤ (N - 1).toNat := by simp
  have hN1 : (0 : ℤ) ≤ N - 1 := by linarith
  calc (L.card : ℤ) * R.card
      = ↑(L.card * R.card) := by push_cast; ring
    _ ≤ ↑(N - 1).toNat := by exact_mod_cast hprod_le
    _ = N - 1 := Int.toNat_of_nonneg hN1

/-- **Lindström's cross-difference inequality.** For a Sidon set in {1,...,N}
    of cardinality m, and any k with 1 ≤ k ≤ m, we have (m - k) * k ≤ N - 1.

    This follows from splitting the sorted set into bottom-`k` and top-`(m-k)`
    elements and applying `ordered_cross_diff_le`.

    Reference: Lindström, B. (1969). An inequality for B₂-sequences.
    *J. Combin. Theory*, 6(2), 211–212. -/
theorem IsIntervalSidon.lindstrom_cross_ineq {A : Finset ℤ} {N : ℕ}
    (hIS : IsIntervalSidon (N : ℤ) A) (hN : 1 ≤ N)
    {k : ℕ} (hk : 1 ≤ k) (hkm : k ≤ A.card) :
    (A.card - k) * k ≤ N - 1 := by
  by_cases hkm_eq : k = A.card
  · simp [hkm_eq]
  have hk_lt : k < A.card := lt_of_le_of_ne hkm hkm_eq
  set sorted := A.sort (· ≤ ·)
  have hnd : sorted.Nodup := A.sort_nodup (· ≤ ·)
  have hlen : sorted.length = A.card := Finset.length_sort _
  have hst : sorted.SortedLT := Finset.sortedLT_sort A
  set L := (sorted.take k).toFinset
  set R := (sorted.drop k).toFinset
  have hLA : L ⊆ A := by
    intro x hx; rw [List.mem_toFinset] at hx
    have := List.mem_of_mem_take hx; rwa [Finset.mem_sort] at this
  have hRA : R ⊆ A := by
    intro x hx; rw [List.mem_toFinset] at hx
    have := List.mem_of_mem_drop hx; rwa [Finset.mem_sort] at this
  have hLR : Disjoint L R := by
    rw [Finset.disjoint_left]; intro x hxL hxR
    rw [List.mem_toFinset] at hxL hxR
    exact absurd hxR ((List.disjoint_take_drop hnd (le_refl k)) hxL)
  have hLcard : L.card = k := by
    rw [List.toFinset_card_of_nodup (hnd.sublist (List.take_sublist k sorted))]
    rw [List.length_take, hlen]; exact Nat.min_eq_left (le_of_lt hk_lt)
  have hRcard : R.card = A.card - k := by
    rw [List.toFinset_card_of_nodup (hnd.sublist (List.drop_sublist k sorted))]
    rw [List.length_drop, hlen]
  have hord : ∀ a ∈ L, ∀ b ∈ R, a < b := by
    intro a ha b hb; rw [List.mem_toFinset] at ha hb
    exact sortedLT_take_lt_drop sorted hst (by rw [hlen]; omega) a ha b hb
  have hLne : L.Nonempty := by
    rw [Finset.nonempty_iff_ne_empty]; intro h; simp [h] at hLcard; omega
  have hRne : R.Nonempty := by
    rw [Finset.nonempty_iff_ne_empty]; intro h; simp [h] at hRcard; omega
  have hN_int : (1 : ℤ) ≤ (N : ℤ) := by exact_mod_cast hN
  have hint : (L.card : ℤ) * R.card ≤ (N : ℤ) - 1 :=
    hIS.ordered_cross_diff_le hLA hRA hLR hord hN_int hLne hRne
  rw [hLcard, hRcard] at hint
  zify [hN, show k ≤ A.card from le_of_lt hk_lt]
  have hconv : (↑(A.card - k) : ℤ) = (↑A.card : ℤ) - ↑k :=
    Nat.cast_sub (le_of_lt hk_lt)
  rw [hconv] at hint
  linarith

/-! ## Lindström Improved Bound — Johnson/Cauchy-Schwarz machinery -/

/-- For a Sidon set A, the shifted copies A+h₁ and A+h₂ intersect in at most
    one element when h₁ ≠ h₂. -/
theorem IsSidon.shift_inter_le_one {A : Finset ℤ} (hA : IsSidon A)
    {h₁ h₂ : ℤ} (hne : h₁ ≠ h₂) :
    ((A.image (· + h₁)) ∩ (A.image (· + h₂))).card ≤ 1 := by
  by_contra hgt
  push_neg at hgt
  obtain ⟨x, hx, y, hy, hxy⟩ := Finset.one_lt_card.mp (by omega : 1 < ((A.image (· + h₁)) ∩ (A.image (· + h₂))).card)
  rw [Finset.mem_inter, Finset.mem_image, Finset.mem_image] at hx hy
  obtain ⟨⟨a₁, ha₁, rfl⟩, ⟨b₁, hb₁, hx_eq⟩⟩ := hx
  obtain ⟨⟨a₂, ha₂, rfl⟩, ⟨b₂, hb₂, hy_eq⟩⟩ := hy
  have heq1 : a₁ - b₁ = h₂ - h₁ := by linarith [hx_eq]
  have heq2 : a₂ - b₂ = h₂ - h₁ := by linarith [hy_eq]
  have hsum : a₁ + b₂ = a₂ + b₁ := by linarith
  rcases hA ha₁ hb₂ ha₂ hb₁ hsum with ⟨h1, h2⟩ | ⟨h1, h2⟩
  · have : a₁ + h₁ = a₂ + h₁ := by rw [h1]
    exact hxy this
  · have : h₁ = h₂ := by linarith [heq1, h1]
    exact hne this

/-- **Johnson's bound (numerical form).** If (km)² ≤ v·m·(m+k-1), then
    k²·m ≤ v·(m+k-1). -/
theorem johnson_numerical {k m v : ℕ} (hm : 0 < m)
    (hcs_moment : (k * m) ^ 2 ≤ v * (m * (m + k - 1))) :
    k ^ 2 * m ≤ v * (m + k - 1) := by
  have hrw1 : (k * m) ^ 2 = k ^ 2 * m * m := by ring
  have hrw2 : v * (m * (m + k - 1)) = v * (m + k - 1) * m := by ring
  rw [hrw1, hrw2] at hcs_moment
  exact Nat.le_of_mul_le_mul_right hcs_moment hm

/-- **Incidence inequality.** For any family of finsets S₁,...,Sₘ all contained
    in a universe U, (∑ᵢ |Sᵢ|)² ≤ |U| · ∑ᵢ ∑ⱼ |Sᵢ ∩ Sⱼ|.
    Uses Cauchy-Schwarz via the degree function d(x) = #{i : x ∈ Sᵢ}. -/
theorem incidence_inequality {α : Type*} [DecidableEq α] (m : ℕ)
    (shifts : Fin m → Finset α) (U : Finset α)
    (hsub : ∀ i, shifts i ⊆ U) :
    (∑ i : Fin m, (shifts i).card) ^ 2 ≤
      U.card * ∑ i : Fin m, ∑ j : Fin m, ((shifts i) ∩ (shifts j)).card := by
  set deg : α → ℕ := fun x => (Finset.univ.filter (fun i : Fin m => x ∈ shifts i)).card
  have hfilt_eq : ∀ i : Fin m, shifts i = U.filter (· ∈ shifts i) := by
    intro i; ext x; simp only [Finset.mem_filter]
    exact ⟨fun h => ⟨hsub i h, h⟩, fun h => h.2⟩
  have h_dc : ∑ i : Fin m, (shifts i).card = ∑ x ∈ U, deg x := by
    simp only [deg]
    trans ∑ i : Fin m, ∑ x ∈ U, if x ∈ shifts i then 1 else 0
    · congr 1; ext i
      conv_lhs => rw [hfilt_eq i, Finset.card_eq_sum_ones, Finset.sum_filter]
    · rw [Finset.sum_comm]
      congr 1; ext x; rw [Finset.card_eq_sum_ones, Finset.sum_filter]
  have hinter_eq : ∀ i j : Fin m, (shifts i) ∩ (shifts j) =
      U.filter (fun x => x ∈ shifts i ∧ x ∈ shifts j) := by
    intro i j; ext x; simp only [Finset.mem_inter, Finset.mem_filter]
    exact ⟨fun ⟨hi, hj⟩ => ⟨hsub i hi, hi, hj⟩, fun ⟨_, hi, hj⟩ => ⟨hi, hj⟩⟩
  have h_sm : ∑ i : Fin m, ∑ j : Fin m, ((shifts i) ∩ (shifts j)).card =
      ∑ x ∈ U, deg x ^ 2 := by
    simp only [deg, sq]
    trans ∑ x ∈ U, ∑ i : Fin m, ∑ j : Fin m,
        if x ∈ shifts i ∧ x ∈ shifts j then (1 : ℕ) else 0
    · trans ∑ i : Fin m, ∑ j : Fin m, ∑ x ∈ U,
          if x ∈ shifts i ∧ x ∈ shifts j then (1 : ℕ) else 0
      · congr 1; ext i; congr 1; ext j
        conv_lhs => rw [hinter_eq i j, Finset.card_eq_sum_ones, Finset.sum_filter]
      · conv_lhs => arg 2; ext i; rw [Finset.sum_comm]
        exact Finset.sum_comm
    · congr 1; ext x
      trans (∑ i : Fin m, if x ∈ shifts i then (1 : ℕ) else 0) *
            (∑ j : Fin m, if x ∈ shifts j then 1 else 0)
      · rw [Finset.sum_mul]; congr 1; ext i; rw [Finset.mul_sum]
        congr 1; ext j
        by_cases h1 : x ∈ shifts i <;> by_cases h2 : x ∈ shifts j <;> simp [h1, h2]
      · congr 1 <;> rw [Finset.card_eq_sum_ones, Finset.sum_filter]
  have h_cs : (∑ x ∈ U, deg x) ^ 2 ≤ U.card * ∑ x ∈ U, deg x ^ 2 := by
    suffices h : (∑ x ∈ U, (deg x : ℤ)) ^ 2 ≤ ↑U.card * ∑ x ∈ U, (deg x : ℤ) ^ 2 by
      exact_mod_cast h
    exact sq_sum_le_card_mul_sum_sq
  calc (∑ i : Fin m, (shifts i).card) ^ 2
      = (∑ x ∈ U, deg x) ^ 2 := by rw [h_dc]
    _ ≤ U.card * ∑ x ∈ U, deg x ^ 2 := h_cs
    _ = U.card * ∑ i : Fin m, ∑ j : Fin m, ((shifts i) ∩ (shifts j)).card := by rw [h_sm]

/-- The intersection matrix row-sum bound for shifted Sidon copies.
    Diagonal terms contribute k each, off-diagonal ≤ 1 each,
    total ≤ mk + m(m-1) = m(m+k-1). -/
theorem sidon_intersection_sum_bound {A : Finset ℤ} (hA : IsSidon A) (m : ℕ) :
    ∑ i : Fin m, ∑ j : Fin m,
      ((A.image (· + (↑(i : ℕ) : ℤ))) ∩ (A.image (· + (↑(j : ℕ) : ℤ)))).card
      ≤ m * (m + A.card - 1) := by
  set k := A.card
  have hrow : ∀ i : Fin m,
      ∑ j : Fin m,
        ((A.image (· + (↑(i : ℕ) : ℤ))) ∩ (A.image (· + (↑(j : ℕ) : ℤ)))).card
        ≤ m + k - 1 := by
    intro i
    have hi_mem : i ∈ (Finset.univ : Finset (Fin m)) := Finset.mem_univ i
    rw [← Finset.sum_erase_add _ _ hi_mem]
    have hdiag : ((A.image (· + (↑(i : ℕ) : ℤ))) ∩ (A.image (· + (↑(i : ℕ) : ℤ)))).card = k := by
      rw [Finset.inter_self]
      exact Finset.card_image_of_injective _ (fun a b h => by linarith)
    have hoff : ∑ j ∈ Finset.univ.erase i,
        ((A.image (· + (↑(i : ℕ) : ℤ))) ∩ (A.image (· + (↑(j : ℕ) : ℤ)))).card ≤ m - 1 := by
      calc ∑ j ∈ Finset.univ.erase i, _
          ≤ ∑ j ∈ Finset.univ.erase i, 1 := Finset.sum_le_sum (fun j hj => by
            have hjmem := Finset.mem_erase.mp hj
            have hne : (↑(i : ℕ) : ℤ) ≠ ↑(j : ℕ) := by
              exact_mod_cast Fin.val_ne_of_ne (Ne.symm hjmem.1)
            exact hA.shift_inter_le_one hne)
        _ = (Finset.univ.erase i).card := by simp
        _ = m - 1 := by simp [Finset.card_erase_of_mem hi_mem, Fintype.card_fin]
    have hm_pos : 0 < m := Fin.pos i
    omega
  calc ∑ i : Fin m, ∑ j : Fin m, _
      ≤ ∑ i : Fin m, (m + k - 1) := Finset.sum_le_sum (fun i _ => hrow i)
    _ = m * (m + k - 1) := by simp [Finset.sum_const, Finset.card_univ, Fintype.card_fin]

/-- **Johnson bound for shifted Sidon sets.** For a Sidon set A ⊆ {1,...,N}
    with |A| = k, the m shifted copies A, A+1, ..., A+(m-1) satisfy
    k²·m ≤ (N+m-1)·(m+k-1). -/
theorem IsIntervalSidon.sidon_johnson_bound {A : Finset ℤ} {N : ℕ}
    (hIS : IsIntervalSidon (N : ℤ) A) (hN : 1 ≤ N)
    (m : ℕ) (hm : 0 < m) :
    A.card ^ 2 * m ≤ (N + m - 1) * (m + A.card - 1) := by
  set k := A.card
  have hA := hIS.sidon
  set shifts : Fin m → Finset ℤ := fun i => A.image (· + (↑(i : ℕ) : ℤ))
  set U := Finset.Icc (1 : ℤ) (↑N + ↑m - 1)
  have hshift_card : ∀ i : Fin m, (shifts i).card = k := fun i =>
    Finset.card_image_of_injective _ (fun a b h => by linarith)
  have hsub : ∀ i : Fin m, shifts i ⊆ U := by
    intro i x hx
    simp only [shifts, Finset.mem_image] at hx
    obtain ⟨a, ha, rfl⟩ := hx
    have hAint := hIS.subset a ha
    simp only [U, Finset.mem_Icc]
    constructor
    · linarith [hAint.1, (i : ℕ).zero_le]
    · have hi : (↑(i : ℕ) : ℤ) ≤ ↑m - 1 := by
        have := i.isLt
        omega
      linarith [hAint.2]
  have hU_card : U.card = N + m - 1 := by
    simp only [U, Int.card_Icc]
    omega
  have hinc := incidence_inequality m shifts U hsub
  have hsum_card : ∑ i : Fin m, (shifts i).card = m * k := by
    simp [hshift_card, Finset.sum_const, Finset.card_univ, Fintype.card_fin]
  have hint : ∑ i : Fin m, ∑ j : Fin m, ((shifts i) ∩ (shifts j)).card
      ≤ m * (m + k - 1) := sidon_intersection_sum_bound hA m
  have hkey : (k * m) ^ 2 ≤ (N + m - 1) * (m * (m + k - 1)) :=
    calc (k * m) ^ 2 = (m * k) ^ 2 := by ring
      _ = (∑ i : Fin m, (shifts i).card) ^ 2 := by rw [hsum_card]
      _ ≤ U.card * ∑ i : Fin m, ∑ j : Fin m, ((shifts i) ∩ (shifts j)).card := hinc
      _ ≤ U.card * (m * (m + k - 1)) := Nat.mul_le_mul_left _ hint
      _ = (N + m - 1) * (m * (m + k - 1)) := by rw [hU_card]
  exact johnson_numerical hm hkey

/-- Johnson bound with Nat subtraction implies the cleaner relaxed bound. -/
theorem lindstrom_monotone {k m n : ℕ}
    (hjohnson : k ^ 2 * m ≤ (n + m - 1) * (m + k - 1)) :
    k ^ 2 * m ≤ (n + m) * (m + k) :=
  hjohnson.trans (Nat.mul_le_mul (Nat.sub_le _ _) (Nat.sub_le _ _))

set_option maxHeartbeats 1600000 in
/-- **Lindström's upper bound.** For a Sidon set A ⊆ {1,...,N} with N ≥ 16,
    |A| ≤ √N + ⁴√N + 2.
    Uses the Johnson bound with optimal choice m = √N · ⁴√N ≈ N^{3/4}. -/
theorem IsIntervalSidon.lindstrom_bound {A : Finset ℤ} {N : ℕ}
    (hIS : IsIntervalSidon (N : ℤ) A) (hN : 16 ≤ N) :
    A.card ≤ Nat.sqrt N + Nat.sqrt (Nat.sqrt N) + 2 := by
  by_contra h_bad
  push_neg at h_bad
  set k := A.card
  set s := Nat.sqrt N
  set t := Nat.sqrt s
  set m := s * t
  have hs_ge : 4 ≤ s := Nat.le_sqrt.mpr (by omega : 4 ^ 2 ≤ N)
  have ht_ge : 2 ≤ t := Nat.le_sqrt.mpr (by omega : 2 ^ 2 ≤ s)
  have hm_pos : 0 < m := by positivity
  have hk_ge : s + t + 3 ≤ k := by omega
  have hN_lt : N < (s + 1) ^ 2 := Nat.lt_succ_sqrt' N
  have hs_lt : s < (t + 1) ^ 2 := Nat.lt_succ_sqrt' s
  have hN_le : N ≤ s ^ 2 + 2 * s := by nlinarith [hN_lt]
  have hs_le : s ≤ t ^ 2 + 2 * t := by nlinarith [hs_lt]
  have hJ := hIS.sidon_johnson_bound (by omega : 1 ≤ N) m hm_pos
  have hJz : (k : ℤ) ^ 2 * ((s : ℤ) * t) ≤
      ((N : ℤ) + s * t - 1) * (s * t + k - 1) := by
    have h : k ^ 2 * (s * t) ≤ (N + s * t - 1) * (s * t + k - 1) := hJ
    have hge1 : 1 ≤ N + s * t := by omega
    have hge2 : 1 ≤ s * t + k := by omega
    zify [hge1, hge2] at h
    linarith
  have hs_z : (s : ℤ) ≤ (t : ℤ) ^ 2 + 2 * t := by exact_mod_cast hs_le
  have hN_z : (N : ℤ) ≤ (s : ℤ) ^ 2 + 2 * s := by exact_mod_cast hN_le
  have hk_z : (s : ℤ) + t + 3 ≤ (k : ℤ) := by exact_mod_cast hk_ge
  have hs_pos : (0 : ℤ) < (s : ℤ) := by
    linarith [show (4 : ℤ) ≤ (s : ℤ) from by exact_mod_cast hs_ge]
  have ht_pos : (0 : ℤ) < (t : ℤ) := by
    linarith [show (2 : ℤ) ≤ (t : ℤ) from by exact_mod_cast ht_ge]
  have hD0 : ((s : ℤ) + t + 3) ^ 2 * (s * t) >
      ((N : ℤ) + s * t - 1) * (s * t + s + t + 2) := by
    have h_id : ((s : ℤ) + t + 3) ^ 2 * (s * t) -
        ((s : ℤ) ^ 2 + 2 * s + s * t - 1) * (s * t + s + t + 2)
        = ((s : ℤ) ^ 2 + 4 * s) * ((t : ℤ) ^ 2 + 2 * t - s)
          + s * ((t : ℤ) ^ 3 + t ^ 2 - 2 * t - 3) + t + 2 := by ring
    have h1 : (0 : ℤ) ≤ ((s : ℤ) ^ 2 + 4 * s) * ((t : ℤ) ^ 2 + 2 * t - s) := by
      apply mul_nonneg <;> nlinarith
    have h2 : (0 : ℤ) < (s : ℤ) * ((t : ℤ) ^ 3 + t ^ 2 - 2 * t - 3) + t + 2 := by
      have ht_cube : (0 : ℤ) < (t : ℤ) ^ 3 + t ^ 2 - 2 * t - 3 := by
        nlinarith [mul_nonneg (show (0 : ℤ) ≤ t from by linarith)
                              (sq_nonneg ((t : ℤ) - 2))]
      nlinarith
    have h_stpos : (0 : ℤ) ≤ (s : ℤ) * t + s + t + 2 := by positivity
    have h_mono : ((N : ℤ) + s * t - 1) * (s * t + s + t + 2) ≤
        ((s : ℤ) ^ 2 + 2 * s + s * t - 1) * (s * t + s + t + 2) := by
      apply mul_le_mul_of_nonneg_right <;> linarith
    nlinarith
  have hDeriv : (0 : ℤ) ≤ (s : ℤ) * t * (k + s + t + 3) - ((N : ℤ) + s * t - 1) := by
    have ht1 : (1 : ℤ) ≤ t := ht_pos
    have h1 : (s : ℤ) * t * (2 * s + 2 * t + 6) ≤ s * t * (k + s + t + 3) :=
      mul_le_mul_of_nonneg_left (by linarith) (mul_nonneg hs_pos.le ht_pos.le)
    have h2 : (s : ℤ) ^ 2 ≤ s ^ 2 * t := le_mul_of_one_le_right (sq_nonneg _) ht1
    have h3 : (s : ℤ) * t ≤ s * t * t :=
      le_mul_of_one_le_right (mul_nonneg hs_pos.le ht_pos.le) ht1
    have h4 : (s : ℤ) ≤ s * t := le_mul_of_one_le_right hs_pos.le ht1
    nlinarith [h1, h2, h3, h4, hN_z]
  have hFact : (k : ℤ) ^ 2 * (s * t) - ((N : ℤ) + s * t - 1) * (s * t + k - 1)
      = ((s : ℤ) + t + 3) ^ 2 * (s * t) - ((N : ℤ) + s * t - 1) * (s * t + s + t + 2)
        + ((k : ℤ) - s - t - 3) *
          ((s : ℤ) * t * (k + s + t + 3) - ((N : ℤ) + s * t - 1)) := by ring
  have hknn : (0 : ℤ) ≤ (k : ℤ) - s - t - 3 := by linarith
  have hprod : (0 : ℤ) ≤ ((k : ℤ) - s - t - 3) *
      ((s : ℤ) * t * (k + s + t + 3) - ((N : ℤ) + s * t - 1)) :=
    mul_nonneg hknn hDeriv
  linarith

/-- The Lindström upper bound: h(N) ≤ √N + √(√N) + 2 for N ≥ 16. -/
theorem sidonMaximum_le_lindstrom (N : ℕ) (hN : 16 ≤ N) :
    sidonMaximum N ≤ Nat.sqrt N + Nat.sqrt (Nat.sqrt N) + 2 := by
  rcases (sidonMaximum_isSidonMaximum N).1 with ⟨A, hA, hAcard⟩
  have h := hA.lindstrom_bound hN
  omega

/-! ## Singer's Construction -/

/-!
Port of the Singer construction from Erdos30/{Singer, SingerBridge, SingerSidon,
SingerTheorem}.lean (Hulak–Ramos–de Queiroz, arXiv:2605.03274), adapted from
Mathlib v4.29.0 to v4.30.0-rc2 and renamespaced into `Semantics.SidonSets.Singer`.
-/

namespace Singer

set_option maxHeartbeats 8000000
set_option linter.unusedSimpArgs false

open Module Submodule Polynomial LinearMap

variable (p : ℕ) [hp : Fact (Nat.Prime p)]

/-! ### Dimension facts (Erdos30/Singer.lean) -/

theorem finrank_ext : finrank (ZMod p) (GaloisField p 3) = 3 :=
  @GaloisField.finrank p hp 3 (by norm_num)

theorem trace_surjective :
    Function.Surjective (Algebra.trace (ZMod p) (GaloisField p 3)) :=
  Algebra.trace_surjective (ZMod p) (GaloisField p 3)

/-- ker(Tr) has dimension 2 (rank-nullity). -/
theorem finrank_ker_trace :
    finrank (ZMod p) (Algebra.trace (ZMod p) (GaloisField p 3)).ker = 2 := by
  have h_rn := LinearMap.finrank_range_add_finrank_ker
    (Algebra.trace (ZMod p) (GaloisField p 3))
  rw [@GaloisField.finrank p hp 3 (by norm_num)] at h_rn
  have htop : (Algebra.trace (ZMod p) (GaloisField p 3)).range = ⊤ :=
    LinearMap.range_eq_top_of_surjective _ (trace_surjective p)
  rw [show finrank (ZMod p) ↥(Algebra.trace (ZMod p) (GaloisField p 3)).range =
      finrank (ZMod p) (ZMod p) from by rw [htop]; exact finrank_top (ZMod p) (ZMod p),
    finrank_self] at h_rn
  omega

/-! ### Minimal polynomial and linear independence -/

theorem minpoly_degree_eq_three (α : GaloisField p 3)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3))) :
    (minpoly (ZMod p) α).natDegree = 3 := by
  have hint : IsIntegral (ZMod p) α := Algebra.IsIntegral.isIntegral α
  have hdvd : (minpoly (ZMod p) α).natDegree ∣ 3 := by
    have h := minpoly.degree_dvd hint
    rwa [@GaloisField.finrank p hp 3 (by norm_num)] at h
  have hne1 : (minpoly (ZMod p) α).natDegree ≠ 1 := by
    intro h1
    exact hα (IntermediateField.mem_bot.mp (by
      rw [← IntermediateField.finrank_eq_one_iff.mp
        (by rw [IntermediateField.adjoin.finrank hint]; exact h1)]
      exact IntermediateField.subset_adjoin (ZMod p) {α} (Set.mem_singleton α)))
  exact (Nat.Prime.eq_one_or_self_of_dvd (by decide) _ hdvd).resolve_left hne1

/-- {α⁰·v, α¹·v, α²·v} are GF(p)-linearly independent when α ∉ GF(p) and v ≠ 0. -/
theorem linIndep_smul_v (α v : GaloisField p 3)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3)))
    (hv : v ≠ 0) :
    LinearIndependent (ZMod p) (fun i : Fin 3 => α ^ (i : ℕ) * v) := by
  rw [Fintype.linearIndependent_iff]
  intro g hg
  have hfactor : (∑ i : Fin 3, g i • α ^ (i : ℕ)) * v = 0 := by
    have heq : ∑ i : Fin 3, g i • (α ^ (i : ℕ) * v) =
               (∑ i : Fin 3, g i • α ^ (i : ℕ)) * v := by
      simp [Finset.sum_mul, Algebra.smul_mul_assoc]
    rw [← heq]; exact hg
  have hsum : ∑ i : Fin 3, g i • α ^ (i : ℕ) = 0 :=
    (mul_eq_zero.mp hfactor).resolve_right hv
  have hdeg := minpoly_degree_eq_three p α hα
  have hli := @linearIndependent_pow _ _ (ZMod p) _ _ α
  rw [Fintype.linearIndependent_iff] at hli
  intro i
  have hsum_t : ∑ j : Fin (minpoly (ZMod p) α).natDegree,
      (g ∘ Fin.cast hdeg) j • α ^ (j : ℕ) = 0 := by
    convert hsum using 1
    exact Fintype.sum_equiv (Fin.castOrderIso hdeg).toEquiv _ _
      (fun j => by simp [Function.comp, Fin.castOrderIso, Fin.cast])
  have h := hli _ hsum_t (Fin.cast hdeg.symm i)
  simp [Function.comp, Fin.cast] at h; exact h

/-! ### No proper invariant subspace -/

/-- Multiplication by α ∉ GF(p) has no proper invariant subspace in GF(p³)/GF(p). -/
theorem no_proper_invariant_subspace (α : GaloisField p 3)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3)))
    (V : Submodule (ZMod p) (GaloisField p 3)) (hVbot : V ≠ ⊥) (hVtop : V ≠ ⊤)
    (hinv : ∀ v : GaloisField p 3, v ∈ V → α • v ∈ V) : False := by
  have hinv_pow : ∀ (n : ℕ) (v : GaloisField p 3), v ∈ V → α ^ n • v ∈ V := by
    intro n; induction n with
    | zero => intro v hv; simpa using hv
    | succ n ih => intro v hv; have h := ih _ (hinv v hv); rwa [← mul_smul, ← pow_succ] at h
  obtain ⟨v, hv_mem, hv_ne⟩ := Submodule.exists_mem_ne_zero_of_ne_bot hVbot
  have hVlt : finrank (ZMod p) V < 3 := by
    have := finrank_lt_finrank_of_lt (lt_top_iff_ne_top.mpr hVtop)
    rw [finrank_top, @GaloisField.finrank p hp 3 (by norm_num)] at this; exact this
  have hmem : ∀ i : Fin 3, α ^ (i : ℕ) * v ∈ V := fun i => by
    have h := hinv_pow i v hv_mem; rwa [Algebra.smul_def] at h
  have hli : LinearIndependent (ZMod p) (fun i : Fin 3 => α ^ (i : ℕ) * v) :=
    linIndep_smul_v p α v hα hv_ne
  have hli_V : LinearIndependent (ZMod p) (fun i : Fin 3 =>
      (⟨α ^ (i : ℕ) * v, hmem i⟩ : V)) := by
    rw [Fintype.linearIndependent_iff]; intro g hg
    rw [Fintype.linearIndependent_iff] at hli
    apply hli g
    have := congr_arg Subtype.val hg
    simpa using this
  exact absurd hVlt (not_lt.mpr
    (le_of_eq (Fintype.card_fin 3).symm |>.trans hli_V.fintype_card_le_finrank))

/-! ### Intersection dimension -/

/-- Two distinct 2-dim subspaces of GF(p³)/GF(p) intersect in dimension 1. -/
theorem finrank_inf_of_distinct_twodim
    (V W : Submodule (ZMod p) (GaloisField p 3))
    (hV : finrank (ZMod p) V = 2) (hW : finrank (ZMod p) W = 2) (hne : V ≠ W) :
    finrank (ZMod p) ↥(V ⊓ W) = 1 := by
  have hgrass := Submodule.finrank_sup_add_finrank_inf_eq V W
  have h_sup_le : finrank (ZMod p) ↥(V ⊔ W) ≤ 3 := by
    have := Submodule.finrank_le (V ⊔ W)
    rw [@GaloisField.finrank p hp 3 (by norm_num)] at this; exact this
  have hV_lt_sup : V < V ⊔ W := lt_of_le_of_ne le_sup_left (fun heq =>
    hne (eq_of_le_of_finrank_le (heq ▸ le_sup_right) (by omega)).symm)
  have h_sup_gt : 2 < finrank (ZMod p) ↥(V ⊔ W) := by
    have := finrank_lt_finrank_of_lt hV_lt_sup; rw [hV] at this; exact this
  rw [hV, hW] at hgrass; omega

/-! ### Multiplication linear equivalence (Erdos30/SingerBridge.lean) -/

noncomputable instance instFintypeGF3 : Fintype (GaloisField p 3) := Fintype.ofFinite _
noncomputable instance instFintypeGF3units : Fintype (GaloisField p 3)ˣ := Fintype.ofFinite _
noncomputable instance instFintypeZModUnits : Fintype (ZMod p)ˣ := Fintype.ofFinite _

/-- Multiplication by a nonzero element is a GF(p)-linear automorphism. -/
noncomputable def mulLinearEquiv (α : GaloisField p 3) (hα : α ≠ 0) :
    GaloisField p 3 ≃ₗ[ZMod p] GaloisField p 3 where
  toFun := fun x => α * x
  map_add' := mul_add α
  map_smul' := fun r x => by simp [Algebra.smul_def, mul_left_comm]
  invFun := fun x => α⁻¹ * x
  left_inv := fun x => by
    show α⁻¹ * (α * x) = x; rw [← mul_assoc, inv_mul_cancel₀ hα, one_mul]
  right_inv := fun x => by
    show α * (α⁻¹ * x) = x; rw [← mul_assoc, mul_inv_cancel₀ hα, one_mul]

/-- The scaled submodule αV = {αv : v ∈ V}. -/
noncomputable def scaledSubmodule (α : GaloisField p 3) (hα : α ≠ 0)
    (V : Submodule (ZMod p) (GaloisField p 3)) :
    Submodule (ZMod p) (GaloisField p 3) :=
  V.map (mulLinearEquiv p α hα).toLinearMap

lemma mem_scaledSubmodule_iff (α : GaloisField p 3) (hα : α ≠ 0)
    (V : Submodule (ZMod p) (GaloisField p 3)) (x : GaloisField p 3) :
    x ∈ scaledSubmodule p α hα V ↔ ∃ v ∈ V, α * v = x := by
  simp [scaledSubmodule, Submodule.mem_map, mulLinearEquiv]

lemma finrank_scaledSubmodule (α : GaloisField p 3) (hα : α ≠ 0)
    (V : Submodule (ZMod p) (GaloisField p 3)) :
    finrank (ZMod p) (scaledSubmodule p α hα V) = finrank (ZMod p) V :=
  LinearEquiv.finrank_eq ((mulLinearEquiv p α hα).submoduleMap V |>.symm)

/-! ### Trace kernel basic properties -/

private lemma ker_trace_ne_bot :
    (Algebra.trace (ZMod p) (GaloisField p 3)).ker ≠ ⊥ :=
  fun h => by have := finrank_ker_trace p; rw [h, finrank_bot (R := ZMod p) (M := GaloisField p 3)] at this; omega

private lemma ker_trace_ne_top :
    (Algebra.trace (ZMod p) (GaloisField p 3)).ker ≠ ⊤ :=
  fun h => by
    have h2 := finrank_ker_trace p
    have h3 := @GaloisField.finrank p hp 3 (by norm_num)
    rw [h, finrank_top, h3] at h2; omega

/-! ### Non-invariance of trace kernel under non-base multiplication -/

/-- The scaled trace kernel αV ≠ V when α ∉ GF(p). -/
lemma scaledSubmodule_ne_ker_trace (α : GaloisField p 3) (hα_ne : α ≠ 0)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3))) :
    scaledSubmodule p α hα_ne (Algebra.trace (ZMod p) (GaloisField p 3)).ker ≠
    (Algebra.trace (ZMod p) (GaloisField p 3)).ker := by
  set V := (Algebra.trace (ZMod p) (GaloisField p 3)).ker
  intro heq
  have hinv : ∀ v : GaloisField p 3, v ∈ V → α • v ∈ V := by
    intro v hv
    have hmem : α * v ∈ scaledSubmodule p α hα_ne V :=
      (mem_scaledSubmodule_iff p α hα_ne V (α * v)).mpr ⟨v, hv, rfl⟩
    rw [heq] at hmem
    rwa [Algebra.smul_def]
  exact no_proper_invariant_subspace p α hα V (ker_trace_ne_bot p) (ker_trace_ne_top p) hinv

/-- α⁻¹ ∉ GF(p) when α ∉ GF(p). -/
lemma inv_not_in_range (α : GaloisField p 3) (hα_ne : α ≠ 0)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3))) :
    α⁻¹ ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3)) := by
  intro ⟨a, ha⟩
  apply hα
  exact ⟨a⁻¹, by rw [map_inv₀, ha, inv_inv]⟩

/-! ### Intersection dimension (geometric core of Sidon proof) -/

/-- When α ∉ GF(p), V ∩ α⁻¹V has dimension 1.
This is the key geometric fact for Singer's Sidon argument. -/
theorem finrank_inf_scaled_ker_trace (α : GaloisField p 3) (hα_ne : α ≠ 0)
    (hα : α ∉ Set.range (algebraMap (ZMod p) (GaloisField p 3))) :
    finrank (ZMod p) ↥((Algebra.trace (ZMod p) (GaloisField p 3)).ker ⊓
      scaledSubmodule p α⁻¹ (inv_ne_zero hα_ne)
        (Algebra.trace (ZMod p) (GaloisField p 3)).ker) = 1 := by
  set V := (Algebra.trace (ZMod p) (GaloisField p 3)).ker
  have hV2 : finrank (ZMod p) V = 2 := finrank_ker_trace p
  have hαinv_ne := inv_ne_zero hα_ne
  have hαinv_not_base := inv_not_in_range p α hα_ne hα
  have hW2 : finrank (ZMod p) (scaledSubmodule p α⁻¹ hαinv_ne V) = 2 := by
    rw [finrank_scaledSubmodule, hV2]
  have hne : V ≠ scaledSubmodule p α⁻¹ hαinv_ne V :=
    fun h => scaledSubmodule_ne_ker_trace p α⁻¹ hαinv_ne hαinv_not_base h.symm
  exact finrank_inf_of_distinct_twodim p V _ hV2 hW2 hne

/-! ### Base units subgroup and index -/

/-- GF(p)× embedded in GF(p³)× via algebraMap. -/
noncomputable def baseUnitsSubgroup : Subgroup (GaloisField p 3)ˣ :=
  (Units.map (algebraMap (ZMod p) (GaloisField p 3)).toMonoidHom).range

instance : (baseUnitsSubgroup p).Normal := inferInstance

private lemma units_map_injective :
    Function.Injective (Units.map (algebraMap (ZMod p) (GaloisField p 3)).toMonoidHom) := by
  apply Units.map_injective
  exact (algebraMap (ZMod p) (GaloisField p 3)).injective

lemma baseUnitsSubgroup_card : Nat.card (baseUnitsSubgroup p) = p - 1 := by
  exact (Nat.card_congr
    ((Units.map (algebraMap (ZMod p) (GaloisField p 3)).toMonoidHom).ofInjective
      (units_map_injective p)).toEquiv.symm).trans
    (by rw [Nat.card_units, Nat.card_zmod])

lemma gf3_units_card : Nat.card (GaloisField p 3)ˣ = p ^ 3 - 1 := by
  rw [Nat.card_units, GaloisField.card p 3 (by norm_num)]

/-- The index [GF(p³)× : GF(p)×] = p²+p+1. -/
lemma baseUnitsSubgroup_index (hp' : Nat.Prime p) :
    (baseUnitsSubgroup p).index = p ^ 2 + p + 1 := by
  have hmul := Subgroup.card_mul_index (baseUnitsSubgroup p)
  rw [baseUnitsSubgroup_card, gf3_units_card] at hmul
  have hp1 : 0 < p - 1 := Nat.sub_pos_of_lt hp'.one_lt
  have hfact : p ^ 3 - 1 = (p - 1) * (p ^ 2 + p + 1) := by
    zify [Nat.one_le_pow 3 p hp'.pos, hp'.pos]; ring
  rw [hfact] at hmul
  exact Nat.eq_of_mul_eq_mul_left hp1 hmul

/-- Membership characterization: u ∈ baseUnitsSubgroup iff ↑u ∈ range(algebraMap). -/
lemma mem_baseUnitsSubgroup_iff (u : (GaloisField p 3)ˣ) :
    u ∈ baseUnitsSubgroup p ↔
    (↑u : GaloisField p 3) ∈ Set.range (algebraMap (ZMod p) (GaloisField p 3)) := by
  simp only [baseUnitsSubgroup, MonoidHom.mem_range]
  constructor
  · rintro ⟨v, rfl⟩; exact ⟨v.val, by simp [Units.coe_map]⟩
  · rintro ⟨a, ha⟩
    have ha_ne : a ≠ 0 := by
      intro h; simp [h] at ha; exact Units.ne_zero u ha.symm
    exact ⟨Units.mk0 a ha_ne, Units.ext (by simp [Units.coe_map, ha])⟩

/-! ### Quotient Sidon property (Erdos30/SingerSidon.lean) -/

/-- Map a nonzero element of GF(p³) to its class in the quotient group
    (GaloisField p 3)ˣ / baseUnitsSubgroup p. -/
noncomputable def singerMk (x : GaloisField p 3) (hx : x ≠ 0) :
    (GaloisField p 3)ˣ ⧸ baseUnitsSubgroup p :=
  QuotientGroup.mk (Units.mk0 x hx)

lemma singerMk_eq_iff (a b : GaloisField p 3) (ha : a ≠ 0) (hb : b ≠ 0) :
    singerMk p a ha = singerMk p b hb ↔
    a⁻¹ * b ∈ Set.range (algebraMap (ZMod p) (GaloisField p 3)) := by
  unfold singerMk; rw [QuotientGroup.eq, mem_baseUnitsSubgroup_iff]
  simp [Units.val_inv_eq_inv_val, Units.val_mul, Units.val_mk0]

private lemma proportional_of_finrank_one
    (W : Submodule (ZMod p) (GaloisField p 3))
    (hW : finrank (ZMod p) W = 1)
    (w₁ w₂ : GaloisField p 3) (hw₁ : w₁ ∈ W) (hw₂ : w₂ ∈ W)
    (hw₁_ne : w₁ ≠ 0) (hw₂_ne : w₂ ≠ 0) :
    ∃ c : ZMod p, c ≠ 0 ∧ w₂ = c • w₁ := by
  have hsub : span (ZMod p) {w₁} ≤ W :=
    span_le.mpr (Set.singleton_subset_iff.mpr hw₁)
  have heq' : W = span (ZMod p) ({w₁} : Set (GaloisField p 3)) :=
    (eq_of_le_of_finrank_le hsub (by rw [finrank_span_singleton hw₁_ne]; omega)).symm
  have hmem : w₂ ∈ span (ZMod p) ({w₁} : Set (GaloisField p 3)) := heq' ▸ hw₂
  rw [mem_span_singleton] at hmem
  obtain ⟨c, hc⟩ := hmem
  exact ⟨c, fun h => hw₂_ne (by simp [h] at hc; exact hc.symm), hc.symm⟩

/-- **Singer Sidon property in the quotient group.**
If u*v = α*(w*x) with u,v,w,x nonzero elements of ker(Tr) and α ∈ (ZMod p)×,
then in the quotient (GaloisField p 3)ˣ / (ZMod p)ˣ we have either
mk u = mk w ∧ mk v = mk x, or mk u = mk x ∧ mk v = mk w. -/
theorem singer_quotient_sidon
    (u v w x : GaloisField p 3)
    (hu : u ∈ (Algebra.trace (ZMod p) (GaloisField p 3)).ker)
    (hv : v ∈ (Algebra.trace (ZMod p) (GaloisField p 3)).ker)
    (hw : w ∈ (Algebra.trace (ZMod p) (GaloisField p 3)).ker)
    (hx : x ∈ (Algebra.trace (ZMod p) (GaloisField p 3)).ker)
    (hu0 : u ≠ 0) (hv0 : v ≠ 0) (hw0 : w ≠ 0) (hx0 : x ≠ 0)
    (α : ZMod p) (hα : α ≠ 0)
    (hmul : u * v = (algebraMap (ZMod p) (GaloisField p 3) α) * (w * x)) :
    (singerMk p u hu0 = singerMk p w hw0 ∧ singerMk p v hv0 = singerMk p x hx0) ∨
    (singerMk p u hu0 = singerMk p x hx0 ∧ singerMk p v hv0 = singerMk p w hw0) := by
  set V := (Algebra.trace (ZMod p) (GaloisField p 3)).ker
  have hαF_ne : (algebraMap (ZMod p) (GaloisField p 3) α) ≠ 0 := by simp [hα]
  have hβ_ne : u * w⁻¹ ≠ 0 := mul_ne_zero hu0 (inv_ne_zero hw0)
  have h_key : (u * w⁻¹) * v = (algebraMap (ZMod p) (GaloisField p 3) α) * x := by
    have h := hmul; field_simp at h ⊢; linear_combination h
  by_cases hβ_base : (u * w⁻¹) ∈ Set.range (algebraMap (ZMod p) (GaloisField p 3))
  case pos =>
    left; constructor
    · rw [singerMk_eq_iff]; obtain ⟨c, hc⟩ := hβ_base
      refine ⟨c⁻¹, ?_⟩
      simp only [map_inv₀]; rw [hc]; field_simp
    · rw [singerMk_eq_iff]; obtain ⟨c, hc⟩ := hβ_base
      have hcF_ne : (algebraMap (ZMod p) (GaloisField p 3) c) ≠ 0 :=
        fun heq => hβ_ne (by rw [← hc, heq])
      refine ⟨α⁻¹ * c, ?_⟩
      simp only [map_mul, map_inv₀]
      have h := h_key; rw [← hc] at h
      generalize algebraMap (ZMod p) (GaloisField p 3) α = αF at h hαF_ne ⊢
      generalize algebraMap (ZMod p) (GaloisField p 3) c = cF at h hcF_ne ⊢
      field_simp; linear_combination h
  case neg =>
    right
    have h1dim : finrank (ZMod p) ↥(V ⊓ scaledSubmodule p (u * w⁻¹)⁻¹ (inv_ne_zero hβ_ne) V) = 1 :=
      finrank_inf_scaled_ker_trace p (u * w⁻¹) hβ_ne hβ_base
    have hw_inf : w ∈ V ⊓ scaledSubmodule p (u * w⁻¹)⁻¹ (inv_ne_zero hβ_ne) V := by
      refine Submodule.mem_inf.mpr ⟨hw, ?_⟩
      rw [mem_scaledSubmodule_iff]
      exact ⟨u, hu, by field_simp⟩
    have hv_inf : v ∈ V ⊓ scaledSubmodule p (u * w⁻¹)⁻¹ (inv_ne_zero hβ_ne) V := by
      refine Submodule.mem_inf.mpr ⟨hv, ?_⟩
      rw [mem_scaledSubmodule_iff]
      refine ⟨(u * w⁻¹) * v, ?_, by field_simp⟩
      rw [h_key, show (algebraMap (ZMod p) (GaloisField p 3) α) * x = α • x
        from (Algebra.smul_def α x).symm]
      exact V.smul_mem α hx
    obtain ⟨c, hc_ne, hvc⟩ := proportional_of_finrank_one p _ h1dim w v hw_inf hv_inf hw0 hv0
    have hcF_ne : (algebraMap (ZMod p) (GaloisField p 3) c) ≠ 0 := by simp [hc_ne]
    have hvc' : v = (algebraMap (ZMod p) (GaloisField p 3) c) * w := by
      rw [hvc, Algebra.smul_def]
    have hcu : (algebraMap (ZMod p) (GaloisField p 3) c) * u =
        (algebraMap (ZMod p) (GaloisField p 3) α) * x := by
      have h := h_key; rw [hvc'] at h
      generalize algebraMap (ZMod p) (GaloisField p 3) α = αF at h hαF_ne ⊢
      generalize algebraMap (ZMod p) (GaloisField p 3) c = cF at h hcF_ne hvc' ⊢
      field_simp at h ⊢; linear_combination h
    constructor
    · rw [singerMk_eq_iff]; refine ⟨c * α⁻¹, ?_⟩
      simp only [map_mul, map_inv₀]
      generalize algebraMap (ZMod p) (GaloisField p 3) α = αF at hαF_ne hcu ⊢
      generalize algebraMap (ZMod p) (GaloisField p 3) c = cF at hcF_ne hcu ⊢
      field_simp; linear_combination hcu
    · rw [singerMk_eq_iff]; refine ⟨c⁻¹, ?_⟩
      simp only [map_inv₀]
      rw [hvc']
      generalize algebraMap (ZMod p) (GaloisField p 3) c = cF at hcF_ne ⊢
      field_simp

/-! ### Combinatorial bridge (Erdos30/SingerTheorem.lean) -/

noncomputable section

private abbrev V' := (Algebra.trace (ZMod p) (GaloisField p 3)).ker
private abbrev Q' := (GaloisField p 3)ˣ ⧸ baseUnitsSubgroup p

private def kerBasis' :=
  (Module.finBasis (ZMod p) ↥(V' p)).reindex
    ((Fin.castOrderIso (finrank_ker_trace p)).toEquiv)

private def repV (i : Option (ZMod p)) : V' p :=
  match i with
  | none => (kerBasis' p) ⟨0, by omega⟩
  | some c => (kerBasis' p) ⟨1, by omega⟩ + c • (kerBasis' p) ⟨0, by omega⟩

private def rep (i : Option (ZMod p)) : GaloisField p 3 := (repV p i).val

private lemma rep_mem (i : Option (ZMod p)) :
    rep p i ∈ (Algebra.trace (ZMod p) (GaloisField p 3)).ker := (repV p i).2

private lemma rep_ne_zero (i : Option (ZMod p)) : rep p i ≠ 0 := by
  intro h; have hV : repV p i = 0 := Subtype.ext h
  cases i with
  | none => exact (kerBasis' p).ne_zero ⟨0, by omega⟩ hV
  | some c =>
    have heq := congr_arg (kerBasis' p).repr hV
    simp only [repV, map_add, map_smul, (kerBasis' p).repr_self, map_zero] at heq
    have h1 := DFunLike.congr_fun heq ⟨1, by omega⟩
    simp only [Finsupp.add_apply, Finsupp.smul_apply, Finsupp.single_apply, smul_eq_mul,
      Finsupp.zero_apply,
      show (⟨1, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2) from rfl,
      show ((⟨0, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2)) = False from by simp [Fin.ext_iff],
      ite_true, ite_false, mul_zero, add_zero] at h1
    exact one_ne_zero h1

private lemma rep_proportional_imp_eq (i j : Option (ZMod p)) (α : ZMod p)
    (hprop : repV p j = α • repV p i) : i = j := by
  have heq := congr_arg (kerBasis' p).repr hprop
  cases i with
  | none =>
    cases j with
    | none => rfl
    | some c =>
      simp only [repV, map_add, map_smul, (kerBasis' p).repr_self] at heq
      have h1 := DFunLike.congr_fun heq ⟨1, by omega⟩
      simp only [Finsupp.add_apply, Finsupp.smul_apply, Finsupp.single_apply, smul_eq_mul,
        show (⟨1, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2) from rfl,
        show ((⟨0, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2)) = False from by simp [Fin.ext_iff],
        ite_true, ite_false, mul_zero, add_zero, mul_one] at h1
      exact absurd h1 one_ne_zero
  | some a =>
    cases j with
    | none =>
      simp only [repV, map_add, map_smul, (kerBasis' p).repr_self] at heq
      have h1 := DFunLike.congr_fun heq ⟨1, by omega⟩
      simp only [Finsupp.add_apply, Finsupp.smul_apply, Finsupp.single_apply, smul_eq_mul,
        show (⟨1, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2) from rfl,
        show ((⟨0, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2)) = False from by simp [Fin.ext_iff],
        ite_true, ite_false, mul_zero, add_zero, mul_one] at h1
      exfalso; exact rep_ne_zero p none (show rep p none = 0 from by
        show (repV p none).val = 0
        have := hprop; rw [show α = 0 from h1.symm, zero_smul] at this
        exact congrArg Subtype.val this)
    | some b =>
      simp only [repV, map_add, map_smul, (kerBasis' p).repr_self] at heq
      have h1 := DFunLike.congr_fun heq ⟨1, by omega⟩
      have h0 := DFunLike.congr_fun heq ⟨0, by omega⟩
      simp only [Finsupp.add_apply, Finsupp.smul_apply, Finsupp.single_apply, smul_eq_mul,
        show (⟨1, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2) from rfl,
        show ((⟨0, by omega⟩ : Fin 2) = (⟨1, by omega⟩ : Fin 2)) = False from by simp [Fin.ext_iff],
        show (⟨0, by omega⟩ : Fin 2) = (⟨0, by omega⟩ : Fin 2) from rfl,
        show ((⟨1, by omega⟩ : Fin 2) = (⟨0, by omega⟩ : Fin 2)) = False from by simp [Fin.ext_iff],
        ite_true, ite_false, mul_zero, add_zero, mul_one, zero_add] at h1 h0
      congr 1; rw [← h1, one_mul] at h0; exact h0.symm

private lemma singerMk_rep_injective :
    Function.Injective (fun i => singerMk p (rep p i) (rep_ne_zero p i)) := by
  intro i j hij
  rw [singerMk_eq_iff] at hij
  obtain ⟨α, hα⟩ := hij
  have hri := rep_ne_zero p i
  have hprop_field : rep p j = (algebraMap (ZMod p) (GaloisField p 3) α) * rep p i := by
    have h1 : (rep p i) * ((rep p i)⁻¹ * rep p j) = rep p j := by
      rw [← mul_assoc, mul_inv_cancel₀ hri, one_mul]
    rw [← h1, hα]; ring
  have hV : repV p j = α • repV p i := by
    apply Subtype.ext
    change rep p j = (α • repV p i).val
    rw [show (α • repV p i).val = α • (repV p i).val from rfl]
    rw [show α • (repV p i).val = (algebraMap (ZMod p) (GaloisField p 3) α) * (repV p i).val
      from Algebra.smul_def α _]
    exact hprop_field
  exact rep_proportional_imp_eq p i j α hV

/-! ### Cyclic group isomorphism -/

private lemma Q_card_eq (hp' : Nat.Prime p) : Nat.card (Q' p) = p ^ 2 + p + 1 := by
  rw [show Nat.card (Q' p) = (baseUnitsSubgroup p).index from
    (Subgroup.index_eq_card _).symm]
  exact baseUnitsSubgroup_index p hp'

private noncomputable def mulEquivQ (hp' : Nat.Prime p) :
    Multiplicative (ZMod (p ^ 2 + p + 1)) ≃* Q' p := by
  haveI : NeZero (p ^ 2 + p + 1) := ⟨by omega⟩
  haveI : IsCyclic (Q' p) :=
    isCyclic_of_surjective (QuotientGroup.mk' (baseUnitsSubgroup p))
      (QuotientGroup.mk'_surjective _)
  let g := Classical.choose (IsCyclic.exists_generator (α := Q' p))
  have hg : ∀ x : Q' p, x ∈ Subgroup.zpowers g :=
    Classical.choose_spec (IsCyclic.exists_generator (α := Q' p))
  have htop : Subgroup.zpowers g = ⊤ := by ext x; exact ⟨fun _ => trivial, fun _ => hg x⟩
  have hord : orderOf g = p ^ 2 + p + 1 := by
    have h1 : Nat.card ↥(Subgroup.zpowers g) = orderOf g := Nat.card_zpowers g
    have h2 : Nat.card ↥(Subgroup.zpowers g) = Nat.card (Q' p) := by
      rw [htop]; exact Nat.card_congr Subgroup.topEquiv.toEquiv
    have := Q_card_eq p hp'; omega
  let φ : Multiplicative (ZMod (p ^ 2 + p + 1)) →* Q' p :=
    MonoidHom.mk' (fun k => g ^ (ZMod.val (Multiplicative.toAdd k))) (fun a b => by
      show g ^ ZMod.val (Multiplicative.toAdd a + Multiplicative.toAdd b) =
        g ^ ZMod.val (Multiplicative.toAdd a) * g ^ ZMod.val (Multiplicative.toAdd b)
      rw [← pow_add, pow_eq_pow_iff_modEq, hord]
      unfold Nat.ModEq; rw [ZMod.val_add]
      exact Nat.mod_mod_of_dvd _ (dvd_refl _))
  have hφ_apply : ∀ k, φ k = g ^ (ZMod.val (Multiplicative.toAdd k)) := fun _ => rfl
  exact MulEquiv.ofBijective φ ⟨by
    intro a b hab
    have hab' : g ^ (ZMod.val (Multiplicative.toAdd a)) =
        g ^ (ZMod.val (Multiplicative.toAdd b)) := by rw [← hφ_apply, ← hφ_apply]; exact hab
    have hinj := @pow_injOn_Iio_orderOf (Q' p) _ g
    have ha : ZMod.val (Multiplicative.toAdd a) ∈ Set.Iio (orderOf g) := by
      rw [Set.mem_Iio, hord]; exact ZMod.val_lt _
    have hb : ZMod.val (Multiplicative.toAdd b) ∈ Set.Iio (orderOf g) := by
      rw [Set.mem_Iio, hord]; exact ZMod.val_lt _
    cases a; cases b; exact congrArg _ (ZMod.val_injective _ (hinj ha hb hab')),
  by
    intro x
    obtain ⟨m, hm⟩ := (hg x : ∃ m : ℤ, g ^ m = x)
    have hpos : (0 : ℤ) < ((p ^ 2 + p + 1 : ℕ) : ℤ) := by positivity
    have hmod_nn : 0 ≤ m % ((p ^ 2 + p + 1 : ℕ) : ℤ) := Int.emod_nonneg _ (by linarith)
    have hmod_lt : (m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat < p ^ 2 + p + 1 := by
      have := Int.emod_lt_of_pos m hpos; omega
    let k : Multiplicative (ZMod (p ^ 2 + p + 1)) :=
      Multiplicative.ofAdd ((m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat : ZMod (p ^ 2 + p + 1))
    refine ⟨k, ?_⟩
    have hφk : φ k = g ^ (m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat := by
      rw [hφ_apply]; congr 1; exact ZMod.val_natCast_of_lt hmod_lt
    have key : g ^ ((m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat : ℤ) = g ^ m := by
      rw [zpow_eq_zpow_iff_modEq, hord]
      change ((m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat : ℤ) % ((p ^ 2 + p + 1 : ℕ) : ℤ) =
        m % ((p ^ 2 + p + 1 : ℕ) : ℤ)
      rw [Int.toNat_of_nonneg hmod_nn, Int.emod_emod_of_dvd _ dvd_rfl]
    calc φ k = g ^ (m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat := hφk
      _ = g ^ ((m % ((p ^ 2 + p + 1 : ℕ) : ℤ)).toNat : ℤ) := (zpow_natCast g _).symm
      _ = g ^ m := key
      _ = x := hm⟩

/-! ### Finset construction and IsSidonMod proof -/

set_option maxHeartbeats 200000000 in
/-- The Singer Sidon set: a Finset ℤ of size p+1 that is IsSidonMod (p²+p+1). -/
theorem singer_sidon_set_of (hp' : Nat.Prime p) :
    ∃ S : Finset ℤ, IsSidonMod (↑p * ↑p + ↑p + 1 : ℤ) S ∧ S.card = p + 1 := by
  haveI : NeZero (p ^ 2 + p + 1) := ⟨by omega⟩
  -- Cyclic isomorphism
  let φ := mulEquivQ p hp'
  -- Map each representative to its ZMod coordinate via φ⁻¹
  let f : Option (ZMod p) → ℤ := fun i =>
    ↑(ZMod.val (Multiplicative.toAdd (φ.symm (singerMk p (rep p i) (rep_ne_zero p i)))))
  let S : Finset ℤ := Finset.univ.image f
  refine ⟨S, ?_, ?_⟩
  · -- IsSidonMod
    intro a b c d ha hb hc hd hdvd
    -- a, b, c, d ∈ S = image of f
    rw [Finset.mem_image] at ha hb hc hd
    obtain ⟨ia, _, rfl⟩ := ha; obtain ⟨ib, _, rfl⟩ := hb
    obtain ⟨ic, _, rfl⟩ := hc; obtain ⟨id, _, rfl⟩ := hd
    -- Abbreviations for the four quotient elements
    set qa := singerMk p (rep p ia) (rep_ne_zero p ia)
    set qb := singerMk p (rep p ib) (rep_ne_zero p ib)
    set qc := singerMk p (rep p ic) (rep_ne_zero p ic)
    set qd := singerMk p (rep p id) (rep_ne_zero p id)
    -- Step 1: divisibility → ZMod equality
    have hzmod : Multiplicative.toAdd (φ.symm qa) + Multiplicative.toAdd (φ.symm qb) =
        Multiplicative.toAdd (φ.symm qc) + Multiplicative.toAdd (φ.symm qd) := by
      have h0 : ((((ZMod.val (Multiplicative.toAdd (φ.symm qa)) : ℤ) +
          (ZMod.val (Multiplicative.toAdd (φ.symm qb)) : ℤ)) -
          ((ZMod.val (Multiplicative.toAdd (φ.symm qc)) : ℤ) +
          (ZMod.val (Multiplicative.toAdd (φ.symm qd)) : ℤ)) : ℤ) : ZMod (p ^ 2 + p + 1)) = 0 := by
        rw [ZMod.intCast_zmod_eq_zero_iff_dvd]
        convert hdvd using 1
        push_cast; ring
      simp only [Int.cast_sub, Int.cast_add, Int.cast_natCast, ZMod.natCast_zmod_val] at h0
      exact sub_eq_zero.mp h0
    -- Step 2: ZMod equality → Multiplicative equality → Q' product equality
    have hQ : qa * qb = qc * qd := by
      have hmult : φ.symm qa * φ.symm qb = φ.symm qc * φ.symm qd := by
        show Multiplicative.ofAdd (Multiplicative.toAdd (φ.symm qa) +
            Multiplicative.toAdd (φ.symm qb)) =
          Multiplicative.ofAdd (Multiplicative.toAdd (φ.symm qc) +
            Multiplicative.toAdd (φ.symm qd))
        exact congrArg _ hzmod
      have hφ := congrArg φ hmult
      simp only [map_mul, MulEquiv.apply_symm_apply] at hφ
      exact hφ
    -- Step 3: Products become singerMk of field products
    have hmul_l : qa * qb = singerMk p (rep p ia * rep p ib)
        (mul_ne_zero (rep_ne_zero p ia) (rep_ne_zero p ib)) := by
      show QuotientGroup.mk (Units.mk0 _ _) * QuotientGroup.mk (Units.mk0 _ _) =
        QuotientGroup.mk (Units.mk0 _ _)
      rw [← QuotientGroup.mk_mul]; congr 1; ext; rfl
    have hmul_r : qc * qd = singerMk p (rep p ic * rep p id)
        (mul_ne_zero (rep_ne_zero p ic) (rep_ne_zero p id)) := by
      show QuotientGroup.mk (Units.mk0 _ _) * QuotientGroup.mk (Units.mk0 _ _) =
        QuotientGroup.mk (Units.mk0 _ _)
      rw [← QuotientGroup.mk_mul]; congr 1; ext; rfl
    -- Step 4: singerMk equality → algebraMap factor via singerMk_eq_iff
    rw [hmul_l, hmul_r] at hQ
    rw [singerMk_eq_iff] at hQ
    obtain ⟨α, hα⟩ := hQ
    have hab_ne : rep p ia * rep p ib ≠ 0 := mul_ne_zero (rep_ne_zero p ia) (rep_ne_zero p ib)
    have hcd_ne : rep p ic * rep p id ≠ 0 := mul_ne_zero (rep_ne_zero p ic) (rep_ne_zero p id)
    have hα_ne : α ≠ 0 := by
      intro h0; rw [h0, map_zero] at hα
      exact mul_ne_zero (inv_ne_zero hab_ne) hcd_ne hα.symm
    have hcd_eq : rep p ic * rep p id =
        (algebraMap (ZMod p) (GaloisField p 3)) α * (rep p ia * rep p ib) := by
      calc rep p ic * rep p id
          = (rep p ia * rep p ib) * ((rep p ia * rep p ib)⁻¹ * (rep p ic * rep p id)) := by
            rw [← mul_assoc, mul_inv_cancel₀ hab_ne, one_mul]
        _ = (rep p ia * rep p ib) * (algebraMap (ZMod p) (GaloisField p 3)) α := by rw [hα]
        _ = (algebraMap (ZMod p) (GaloisField p 3)) α * (rep p ia * rep p ib) := mul_comm _ _
    -- Step 5: Apply singer_quotient_sidon
    have hsq := singer_quotient_sidon p (rep p ic) (rep p id) (rep p ia) (rep p ib)
      (rep_mem p ic) (rep_mem p id) (rep_mem p ia) (rep_mem p ib)
      (rep_ne_zero p ic) (rep_ne_zero p id) (rep_ne_zero p ia) (rep_ne_zero p ib)
      α hα_ne hcd_eq
    -- Step 6: From singerMk equality to index equality via injectivity
    cases hsq with
    | inl h =>
      left; constructor
      · exact congrArg f (singerMk_rep_injective p h.1.symm)
      · exact congrArg f (singerMk_rep_injective p h.2.symm)
    | inr h =>
      right; constructor
      · exact congrArg f (singerMk_rep_injective p h.2.symm)
      · exact congrArg f (singerMk_rep_injective p h.1.symm)
  · -- card S = p + 1
    rw [Finset.card_image_of_injective _ (by
      intro i j hij
      -- f(i) = f(j) means val(toAdd(φ⁻¹(singerMk(rep i)))) = val(toAdd(φ⁻¹(singerMk(rep j))))
      -- nat→int cast is injective, val is injective, toAdd is bijective, φ⁻¹ is bijective
      -- So singerMk(rep i) = singerMk(rep j), hence i = j by singerMk_rep_injective
      have h1 : (ZMod.val (Multiplicative.toAdd (φ.symm (singerMk p (rep p i) (rep_ne_zero p i)))) : ℤ) =
        ↑(ZMod.val (Multiplicative.toAdd (φ.symm (singerMk p (rep p j) (rep_ne_zero p j))))) := hij
      have h2 := Nat.cast_injective h1
      have h3 := ZMod.val_injective _ h2
      -- h3 : toAdd(φ⁻¹(singerMk(rep i))) = toAdd(φ⁻¹(singerMk(rep j)))
      have h4 : φ.symm (singerMk p (rep p i) (rep_ne_zero p i)) =
          φ.symm (singerMk p (rep p j) (rep_ne_zero p j)) :=
        Multiplicative.toAdd.injective h3
      have h5 := φ.symm.injective h4
      exact singerMk_rep_injective p h5)]
    simp [Finset.card_univ, Fintype.card_option, ZMod.card]

end
end Singer

/-- **Singer's theorem.** For each prime p, there exists a Sidon set
    modulo p² + p + 1 of cardinality p + 1.

    This is the classical algebraic construction using the trace kernel
    of GF(p³)/GF(p). The proof proceeds through:
    1. Construction of GF(p) and its degree-3 extension GF(p³)
    2. Analysis of ker(Tr) as a 2-dimensional subspace
    3. Geometric argument via subspace intersections
    4. Transfer from quotient multiplication to modular integer addition

    Reference: Singer, J. (1938). A theorem in finite projective geometry
    and some applications. *Trans. Amer. Math. Soc.*, 43, 377–385. -/
theorem singer_sidon_set (p : ℕ) (hp : Nat.Prime p) :
    ∃ S : Finset ℤ, IsSidonMod (↑p * ↑p + ↑p + 1 : ℤ) S ∧ S.card = p + 1 := by
  haveI : Fact (Nat.Prime p) := ⟨hp⟩
  exact Singer.singer_sidon_set_of p hp

/-- The Singer family hypothesis: for every prime p, there exists a
    Sidon set mod (p²+p+1) of size p+1. -/
def SingerFamilyHypothesis : Prop :=
  ∀ p : ℕ, Nat.Prime p →
    ∃ S : Finset ℤ, IsSidonMod (↑p * ↑p + ↑p + 1 : ℤ) S ∧ S.card = p + 1

/-- Singer's theorem establishes the Singer family hypothesis. -/
theorem singerFamilyHypothesis_holds : SingerFamilyHypothesis :=
  fun p hp => singer_sidon_set p hp

/-! ## Unconditional h(N) = Θ(√N) Bounds -/

/-- Bounded-lift lemma: an IsSidonMod M set whose elements all lie in [0, M-1]
    is automatically an IsSidon set in ℤ (no wraparound can occur).
    
    Proof: If a+b = c+d + M, then a+b ≡ c+d (mod M), so IsSidonMod gives
    {a,b} = {c,d}. But then a+b = a+b + M ⇒ M = 0 — contradiction.
    Therefore a+b = c+d in ℤ, which is the Sidon property. -/
theorem IsSidonMod.isSidon_of_bounded {M : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M S) (h_bound : ∀ x ∈ S, 0 ≤ x ∧ x < M) : IsSidon S := by
  intro a b c d ha hb hc hd hsum
  have ha_bound := h_bound a ha; have hb_bound := h_bound b hb
  have hc_bound := h_bound c hc; have hd_bound := h_bound d hd
  have ha_nonneg : 0 ≤ a := ha_bound.1; have ha_lt : a < M := ha_bound.2
  have hb_nonneg : 0 ≤ b := hb_bound.1; have hb_lt : b < M := hb_bound.2
  have hc_nonneg : 0 ≤ c := hc_bound.1; have hc_lt : c < M := hc_bound.2
  have hd_nonneg : 0 ≤ d := hd_bound.1; have hd_lt : d < M := hd_bound.2
  -- From IsSidonMod, a+b ≡ c+d (mod M)
  have hmod : M ∣ (a + b) - (c + d) := by
    -- hsum states a + b = c + d in ℤ, so (a+b) - (c+d) = 0 which is divisible by M
    rw [hsum, sub_self]
    exact dvd_zero M
  -- If a+b = c+d, Sidon property follows directly
  rcases hS ha hb hc hd hmod with (⟨hac, hbd⟩ | ⟨had, hbc⟩)
  · left; exact ⟨hac, hbd⟩
  · right; exact ⟨had, hbc⟩



/-- The Sidon maximum is positive for N ≥ 1. -/
theorem sidonMaximum_pos (N : ℕ) (hN : 1 ≤ N) : 1 ≤ sidonMaximum N := by
  have hmax := sidonMaximum_isSidonMaximum N
  have hSidon : IsSidon ({1} : Finset ℤ) := by
    intro a b c d ha hb hc hd _; simp at ha hb hc hd
    left; exact ⟨ha ▸ hc.symm, hb ▸ hd.symm⟩
  have h1 : IsIntervalSidon (N : ℤ) ({1} : Finset ℤ) := by
    constructor
    · intro x hx; simp at hx; subst hx; exact ⟨le_refl 1, by exact_mod_cast hN⟩
    · exact hSidon
  have hle := hmax.2 h1; simp at hle; exact hle

/-- The Sidon maximum function is monotone non-decreasing. -/
theorem sidonMaximum_mono {N M : ℕ} (hNM : N ≤ M) :
    sidonMaximum N ≤ sidonMaximum M := by
  have hmax_N := sidonMaximum_isSidonMaximum N
  have hmax_M := sidonMaximum_isSidonMaximum M
  rcases hmax_N.1 with ⟨A, hA, hAcard⟩
  have hA_M : IsIntervalSidon (M : ℤ) A := hA.mono (by exact_mod_cast hNM)
  have hle := hmax_M.2 hA_M; omega

/-! ## Erdős Problem 30 Statement -/

/-- The formal Erdős Problem 30 statement: h(N) = √N + O_ε(N^ε) for every ε > 0. -/
def Erdos30Statement : Prop :=
  ∀ ε : ℝ, 0 < ε →
    ∃ C : ℝ, ∃ N0 : ℕ,
      0 ≤ C ∧
        ∀ {N h : ℕ}, N0 ≤ N → IsSidonMaximum N h →
          abs ((h : ℝ) - Real.sqrt (N : ℝ)) ≤ C * Real.rpow (N : ℝ) ε

/-- **Partial discharge for ε ≥ 1/2** (unconditional).
    For all ε ≥ 1/2, |h(N) - √N| ≤ 2·N^ε for all N ≥ 5. -/
theorem erdos30_partial_half :
    ∀ ε : ℝ, (1 : ℝ) / 2 ≤ ε → 0 < ε →
      ∃ C : ℝ, ∃ N0 : ℕ,
        0 ≤ C ∧
          ∀ {N h : ℕ}, N0 ≤ N → IsSidonMaximum N h →
            abs ((h : ℝ) - Real.sqrt (N : ℝ)) ≤ C * Real.rpow (N : ℝ) ε :=
by
  intro ε hε_ge_half hε_pos
  have h_two_nonneg : 0 ≤ (2 : ℝ) := by norm_num
  have hC_nonneg : 0 ≤ Real.sqrt 2 := Real.sqrt_nonneg _
  refine ⟨Real.sqrt 2, 1, hC_nonneg, ?_⟩
  intro N h hN1 hmax
  have hN_pos : 1 ≤ N := hN1
  have hN_pos_real : (1 : ℝ) ≤ (N : ℝ) := by exact_mod_cast hN_pos
  -- Quadratic upper bound: h ≤ √(2N) + 1
  have h_bound_nat : h ≤ Nat.sqrt (2 * N) + 1 := by
    rcases hmax.1 with ⟨A, hA, hAcard⟩
    have hcard := hA.card_le hN_pos
    rw [hAcard] at hcard
    exact hcard
  have h_bound_real : (h : ℝ) ≤ (Nat.sqrt (2 * N) : ℝ) + 1 := by exact_mod_cast h_bound_nat
  -- (Nat.sqrt (2*N) : ℝ) ≤ Real.sqrt (2*(N:ℝ))
  have h_sqrt_nat_sq : (Nat.sqrt (2 * N) : ℝ) * (Nat.sqrt (2 * N) : ℝ) ≤ 2 * (N : ℝ) := by
    have h_sq_nat : (Nat.sqrt (2 * N)) * (Nat.sqrt (2 * N)) ≤ 2 * N :=
      (Nat.le_sqrt.1 (le_refl (Nat.sqrt (2 * N))))
    exact_mod_cast h_sq_nat
  have h_nat_sqrt_nonneg : 0 ≤ (Nat.sqrt (2 * N) : ℝ) := by exact_mod_cast Nat.zero_le _
  have h_nat_sqrt_le_real_sqrt : (Nat.sqrt (2 * N) : ℝ) ≤ Real.sqrt (2 * (N : ℝ)) := by
    calc
      (Nat.sqrt (2 * N) : ℝ) = Real.sqrt (((Nat.sqrt (2 * N) : ℝ)) * ((Nat.sqrt (2 * N) : ℝ))) := by
        rw [Real.sqrt_mul_self h_nat_sqrt_nonneg]
      _ ≤ Real.sqrt (2 * (N : ℝ)) := Real.sqrt_le_sqrt h_sqrt_nat_sq
  have h_upper_real_sqrt : (h : ℝ) ≤ Real.sqrt (2 * (N : ℝ)) + 1 := by
    linarith
  have h_pos_nat : 1 ≤ h := by
    have h_eq : h = sidonMaximum N :=
      isSidonMaximum_unique hmax (sidonMaximum_isSidonMaximum N)
    rw [h_eq]; exact sidonMaximum_pos N hN_pos
  have h_lower_one : (1 : ℝ) ≤ (h : ℝ) := by exact_mod_cast h_pos_nat
  have h_N_ge_one_sqrt : 1 ≤ Real.sqrt (N : ℝ) := by
    calc
      (1 : ℝ) = Real.sqrt ((1 : ℝ)) := by norm_num
      _ ≤ Real.sqrt (N : ℝ) := Real.sqrt_le_sqrt (by exact_mod_cast hN_pos)
  have h_sqrt_eq_rpow : Real.sqrt (N : ℝ) = Real.rpow (N : ℝ) ((1 : ℝ) / 2) :=
    Real.sqrt_eq_rpow _
  have h_sqrt_N_le_N_pow_eps : Real.sqrt (N : ℝ) ≤ Real.rpow (N : ℝ) ε := by
    rw [h_sqrt_eq_rpow]
    exact Real.rpow_le_rpow_of_exponent_le hN_pos_real hε_ge_half
  have h_one_le_N_pow_eps : (1 : ℝ) ≤ Real.rpow (N : ℝ) ε := by
    have : (1 : ℝ) ^ ε = (1 : ℝ) := by simp
    have h_rpow_mono : (1 : ℝ) ^ ε ≤ (N : ℝ) ^ ε :=
      Real.rpow_le_rpow (by norm_num) hN_pos_real (hε_pos.le)
    simpa [this] using h_rpow_mono
  have ha_nonneg : 0 ≤ Real.sqrt 2 - 1 := by
    have : 1 ≤ Real.sqrt 2 := by
      calc
        (1 : ℝ) = Real.sqrt (1 : ℝ) := by norm_num
        _ ≤ Real.sqrt 2 := Real.sqrt_le_sqrt (by norm_num)
    linarith
  -- Case 1: h - √N ≥ 0
  by_cases h_nonneg_diff : (h : ℝ) - Real.sqrt (N : ℝ) ≥ 0
  · rw [abs_of_nonneg h_nonneg_diff]
    have h_diff_upper : (h : ℝ) - Real.sqrt (N : ℝ) ≤ Real.sqrt 2 * Real.rpow (N : ℝ) ε := by
      calc
        (h : ℝ) - Real.sqrt (N : ℝ) ≤ (Real.sqrt (2 * (N : ℝ)) + 1) - Real.sqrt (N : ℝ) := by
          linarith
        _ = Real.sqrt (2 * (N : ℝ)) - Real.sqrt (N : ℝ) + 1 := by ring
        _ = Real.sqrt 2 * Real.sqrt (N : ℝ) - Real.sqrt (N : ℝ) + 1 := by
          rw [Real.sqrt_mul (by norm_num : 0 ≤ (2 : ℝ))]
        _ = (Real.sqrt 2 - 1) * Real.sqrt (N : ℝ) + 1 := by ring
        _ ≤ (Real.sqrt 2 - 1) * Real.rpow (N : ℝ) ε + Real.rpow (N : ℝ) ε := by
          nlinarith
        _ = Real.sqrt 2 * Real.rpow (N : ℝ) ε := by ring
    exact h_diff_upper
  · rw [abs_of_neg (by linarith), neg_sub]
    have h_diff_lower : Real.sqrt (N : ℝ) - (h : ℝ) ≤ Real.sqrt 2 * Real.rpow (N : ℝ) ε := by
      calc
        Real.sqrt (N : ℝ) - (h : ℝ) ≤ Real.sqrt (N : ℝ) - 1 := by nlinarith
        _ ≤ Real.sqrt (N : ℝ) := by nlinarith
        _ ≤ Real.rpow (N : ℝ) ε := h_sqrt_N_le_N_pow_eps
        _ ≤ Real.sqrt 2 * Real.rpow (N : ℝ) ε := by nlinarith
    exact h_diff_lower

/-- **Lindström upper bound for ε ≥ 1/4** (unconditional).
    For all ε ≥ 1/4, h(N) ≤ √N + 2·N^ε for all N ≥ 16. -/
theorem sidonUpperBound_quarter :
    ∀ ε : ℝ, (1 : ℝ) / 4 ≤ ε → 0 < ε →
      ∃ C : ℝ, ∃ N0 : ℕ,
        0 ≤ C ∧
          ∀ {N h : ℕ}, N0 ≤ N → IsSidonMaximum N h →
            (h : ℝ) ≤ Real.sqrt (N : ℝ) + C * Real.rpow (N : ℝ) ε :=
by
  intro ε hε_ge_quarter hε_pos
  by_cases hε_ge_half : (1 : ℝ) / 2 ≤ ε
  · -- For ε ≥ 1/2, use the stronger bilateral bound from erdos30_partial_half
    rcases erdos30_partial_half ε hε_ge_half hε_pos with ⟨C, N0, hC_nonneg, hC⟩
    refine ⟨C, N0, hC_nonneg, ?_⟩
    intro N h hN hmax
    have h_abs := hC hN hmax
    have h_abs_le := abs_le.mp h_abs
    nlinarith
  · -- For 1/4 ≤ ε < 1/2, use the Lindström bound h(N) ≤ √N + ⁴√N + 2
    refine ⟨3, 16, by norm_num, ?_⟩
    intro N h hN hmax
    have hN16 : 16 ≤ N := hN
    have hNpos : 1 ≤ N := by omega
    have hN_real : (1 : ℝ) ≤ (N : ℝ) := by exact_mod_cast hNpos
    have hN_nonneg : 0 ≤ (N : ℝ) := by exact_mod_cast Nat.zero_le _
    have h_bound_nat : h ≤ Nat.sqrt N + Nat.sqrt (Nat.sqrt N) + 2 := by
      have h_eq : h = sidonMaximum N :=
        isSidonMaximum_unique hmax (sidonMaximum_isSidonMaximum N)
      rw [h_eq]; exact sidonMaximum_le_lindstrom N hN16
    have h_bound_real : (h : ℝ) ≤ (Nat.sqrt N : ℝ) + (Nat.sqrt (Nat.sqrt N) : ℝ) + 2 :=
      by exact_mod_cast h_bound_nat
    have h_sq_s : (Nat.sqrt N : ℝ)^2 ≤ (N : ℝ) := by
      have h := Nat.sqrt_le' N
      exact_mod_cast h
    have h_sqrt_nat_le_real : (Nat.sqrt N : ℝ) ≤ Real.sqrt (N : ℝ) := by
      calc
        (Nat.sqrt N : ℝ) = Real.sqrt (((Nat.sqrt N : ℝ))^2) := by
          rw [Real.sqrt_sq (show 0 ≤ (Nat.sqrt N : ℝ) from by exact_mod_cast Nat.zero_le _)]
        _ ≤ Real.sqrt (N : ℝ) := Real.sqrt_le_sqrt h_sq_s
    have h_t_sq_s : ((Nat.sqrt (Nat.sqrt N) : ℝ))^2 ≤ Real.sqrt (N : ℝ) := by
      calc
        ((Nat.sqrt (Nat.sqrt N) : ℝ))^2 ≤ (Nat.sqrt N : ℝ) := by
          have h := Nat.sqrt_le' (Nat.sqrt N)
          exact_mod_cast h
        _ ≤ Real.sqrt (N : ℝ) := h_sqrt_nat_le_real
    have h_sqrt_sqrt_rpow : Real.sqrt (Real.sqrt (N : ℝ)) = (N : ℝ) ^ ((1 : ℝ) / 4) := by
      calc
        Real.sqrt (Real.sqrt (N : ℝ)) = ((N : ℝ) ^ ((1 : ℝ) / 2)) ^ ((1 : ℝ) / 2) := by
          simp [Real.sqrt_eq_rpow]
        _ = (N : ℝ) ^ (((1 : ℝ) / 2) * ((1 : ℝ) / 2)) := by
          rw [Real.rpow_mul hN_nonneg]
        _ = (N : ℝ) ^ ((1 : ℝ) / 4) := by ring
    have h_t_le_rpow : (Nat.sqrt (Nat.sqrt N) : ℝ) ≤ (N : ℝ) ^ ε := by
      calc
        (Nat.sqrt (Nat.sqrt N) : ℝ) = Real.sqrt (((Nat.sqrt (Nat.sqrt N) : ℝ))^2) := by
          rw [Real.sqrt_sq (show 0 ≤ (Nat.sqrt (Nat.sqrt N) : ℝ) from by exact_mod_cast Nat.zero_le _)]
        _ ≤ Real.sqrt (Real.sqrt (N : ℝ)) := Real.sqrt_le_sqrt h_t_sq_s
        _ = (N : ℝ) ^ ((1 : ℝ) / 4) := h_sqrt_sqrt_rpow
        _ ≤ (N : ℝ) ^ ε := Real.rpow_le_rpow_of_exponent_le hN_real hε_ge_quarter
    have h_two_le_rpow : (2 : ℝ) ≤ 2 * (N : ℝ) ^ ε := by
      have h_one_le : (1 : ℝ) ≤ (N : ℝ) ^ ε := by
        have h_one_rpow : (1 : ℝ) ^ ε = (1 : ℝ) := by simp
        have h_rpow_mono : (1 : ℝ) ^ ε ≤ (N : ℝ) ^ ε :=
          Real.rpow_le_rpow (by norm_num) hN_real (hε_pos.le)
        simpa [h_one_rpow] using h_rpow_mono
      nlinarith
    calc
      (h : ℝ) ≤ (Nat.sqrt N : ℝ) + (Nat.sqrt (Nat.sqrt N) : ℝ) + 2 := h_bound_real
      _ ≤ Real.sqrt (N : ℝ) + (Nat.sqrt (Nat.sqrt N) : ℝ) + 2 := by nlinarith
      _ ≤ Real.sqrt (N : ℝ) + (N : ℝ) ^ ε + 2 := by nlinarith
      _ ≤ Real.sqrt (N : ℝ) + (N : ℝ) ^ ε + 2 * (N : ℝ) ^ ε := by nlinarith
      _ = Real.sqrt (N : ℝ) + 3 * (N : ℝ) ^ ε := by ring

/-! ## Conditional Erdős Problem 30 Reduction

`conditional_erdos30` appears later in this file (after `singerIntervalSidon`),
since its lower-bound side is discharged via the Singer interval construction. -/

/-! ## Representation Function -/

/-- For a Sidon set, the representation function is bounded by 2:
    at most 2 ordered pairs (a,b) ∈ A×A satisfy a + b = n.
    Uses Finset.product instead of the ×ˢ notation. -/
theorem IsSidon.repr_le_two {A : Finset ℤ} (hA : IsSidon A) (n : ℤ) :
    ((A.product A).filter (fun ab : ℤ × ℤ => ab.1 + ab.2 = n)).card ≤ 2 := by
  set S := (A.product A).filter (λ ab : ℤ × ℤ => ab.1 + ab.2 = n) with hS
  by_cases h_empty : S.Nonempty
  · rcases h_empty with ⟨⟨a, b⟩, hab⟩
    have ha_mem_filter : (a, b) ∈ (A.product A).filter (λ ab : ℤ × ℤ => ab.1 + ab.2 = n) := by
      simpa [hS] using hab
    have ha_mem_product : (a, b) ∈ A.product A :=
      (Finset.mem_filter.1 ha_mem_filter).1
    have ha_all : a ∈ A ∧ b ∈ A := Finset.mem_product.1 ha_mem_product
    have ha : a ∈ A := ha_all.1
    have hb : b ∈ A := ha_all.2
    have hsum : a + b = n := by
      simpa using (Finset.mem_filter.1 ha_mem_filter).2
    have h_sub : S ⊆ {(a, b), (b, a)} := by
      intro ⟨x, y⟩ hxy
      have hx_mem_filter : (x, y) ∈ (A.product A).filter (λ ab : ℤ × ℤ => ab.1 + ab.2 = n) := by
        simpa [hS] using hxy
      have hx_mem_product : (x, y) ∈ A.product A :=
        (Finset.mem_filter.1 hx_mem_filter).1
      have hx_all : x ∈ A ∧ y ∈ A := Finset.mem_product.1 hx_mem_product
      have hx : x ∈ A := hx_all.1
      have hy : y ∈ A := hx_all.2
      have hsum_xy : x + y = n := by
        simpa using (Finset.mem_filter.1 hx_mem_filter).2
      have hab_eq : a + b = x + y := by
        calc
          a + b = n := hsum
          _ = x + y := hsum_xy.symm
      rcases hA ha hb hx hy hab_eq with (⟨hac, hbd⟩ | ⟨had, hbc⟩)
      · simp [hac, hbd]
      · simp [had, hbc]
    have h_card_sub : S.card ≤ ({(a, b), (b, a)} : Finset (ℤ × ℤ)).card :=
      Finset.card_le_card h_sub
    have h_card_two : ({(a, b), (b, a)} : Finset (ℤ × ℤ)).card ≤ 2 := by
      by_cases h_eq : (a, b) = (b, a)
      · simp [h_eq]
      · simp [h_eq]
    omega
  · have h_card_zero : S.card = 0 := by
      apply Finset.card_eq_zero.mpr
      ext x; simp; intro hx; exact h_empty ⟨x, hx⟩
    omega

/-! ## No-Wraparound Lemma -/

/-- **No-wraparound lemma.** If all elements of A are in {1,...,N} and
    M ≥ 2N - 1, then IsSidon A → IsSidonMod M A. This is the key step
    that lets interval Sidon sets be embedded into a cyclic ambient group
    without creating new sum collisions. -/
theorem IsSidon.isSidonMod_of_interval {A : Finset ℤ} {N M : ℤ}
    (hA : IsSidon A)
    (hbound : ∀ x ∈ A, 1 ≤ x ∧ x ≤ N)
    (hM : 2 * N - 1 ≤ M) : IsSidonMod M A := by
  intro a b c d ha hb hc hd hdiv
  have ha_bound := hbound a ha
  have hb_bound := hbound b hb
  have hc_bound := hbound c hc
  have hd_bound := hbound d hd
  have hN_pos : 1 ≤ N := le_trans ha_bound.1 ha_bound.2
  have hM_nonneg : 0 ≤ M := by omega
  have hdiff_bound : |(a + b) - (c + d)| ≤ 2 * N - 2 := by
    apply abs_le.mpr
    constructor <;> omega
  have hlt : |(a + b) - (c + d)| < M := by
    have : 2 * N - 2 < 2 * N - 1 := by omega
    omega
  rcases hdiv with ⟨k, hk⟩
  by_cases hk0 : k = 0
  · rw [hk0, mul_zero] at hk
    have hsum_eq : a + b = c + d := by omega
    exact hA ha hb hc hd hsum_eq
  · have hk_abs_ge_one : 1 ≤ |k| := by
      have hk_ne_zero : k ≠ 0 := hk0
      have hk_abs_pos : 0 < |k| := abs_pos.mpr hk_ne_zero
      omega
    have hM_eq_abs : |M| = M := abs_of_nonneg hM_nonneg
    have hdiff_bound_M : M ≤ |(a + b) - (c + d)| := by
      calc
        M = |M| := hM_eq_abs.symm
        _ ≤ |M| * |k| := by
          calc
            |M| = |M| * 1 := by simp
            _ ≤ |M| * |k| :=
              mul_le_mul_of_nonneg_left hk_abs_ge_one (abs_nonneg _)
        _ = |M * k| := by rw [abs_mul]
        _ = |(a + b) - (c + d)| := by rw [hk]
    linarith

/-! ## Singer ↔ Golden Angle Connection -/

/-- The Singer construction modulus for prime p: q² + q + 1 where q = p.
    For p = 2: 2² + 2 + 1 = 7. For p = 3: 3² + 3 + 1 = 13.
    These are the orders of the cyclic difference sets. -/
def singerModulus (p : ℕ) : ℕ := p * p + p + 1

/-- The Singer set cardinality for prime p: p + 1 elements. -/
def singerCardinality (p : ℕ) : ℕ := p + 1

/-- The Singer Sidon density ratio: numerator = p+1, denominator = p²+p+1.
    For large p, this ratio ≈ 1/p → 0, while the golden angle density
    1/φ ≈ 0.618 exceeds all finite Singer densities. -/
def singerDensityNum (p : ℕ) : ℕ := p + 1
def singerDensityDen (p : ℕ) : ℕ := p * p + p + 1

-- Executable witnesses for small primes
#eval singerModulus 2    -- 7
#eval singerCardinality 2  -- 3
#eval singerModulus 3    -- 13
#eval singerCardinality 3  -- 4
#eval singerModulus 5    -- 31
#eval singerCardinality 5  -- 6

/-! ## Cyclic Window Infrastructure (from SidonGap.lean) -/

/-!
Port of the cyclic gap infrastructure from Hulak–Ramos–de Queiroz (2026),
Erdos30/SidonGap.lean. This provides the bridge from modular Sidon sets to
interval Sidon sets via the cyclic gap structure.

Key results:
1. residueImageNat — the natural-number residue image of a modular Sidon set
2. sortedResidues — the sorted residues indexed by Fin S.card
3. sortedResidueWitness — a chosen element realizing each sorted residue
4. cyclicGapAt — the cyclic gap function on sorted residues
5. exists_full_intervalSidon_of_quantitative_gap_bound — the main theorem
-/

-- ============================================================
-- Part 1: Residue image and injectivity
-- ============================================================

/-- The natural-number residue image of `S` modulo `M`. -/
noncomputable def residueImageNat (M : ℤ) (S : Finset ℤ) : Finset ℕ :=
  S.image fun s => Int.toNat (s % M)

theorem mem_residueImageNat {M : ℤ} {S : Finset ℤ} {n : ℕ} :
    n ∈ residueImageNat M S ↔ ∃ s ∈ S, Int.toNat (s % M) = n := by
  simp [residueImageNat]

private theorem exists_mem_modEq_of_mem_residueImageNat
    {M_int : ℤ} {S : Finset ℤ} (hM : 0 < M_int) {n : ℕ}
    (hn : n ∈ residueImageNat M_int S) :
    ∃ s ∈ S, (n : ℤ) = s % M_int := by
  rcases mem_residueImageNat.mp hn with ⟨s, hs, hsmod⟩
  refine ⟨s, hs, ?_⟩
  have hnonneg : 0 ≤ s % M_int := Int.emod_nonneg _ (ne_of_gt hM)
  simpa [Int.toNat_of_nonneg hnonneg] using
    (congrArg (fun m : ℕ => (m : ℤ)) hsmod).symm

/-- The residue map `s ↦ s % M` is injective on any modular Sidon set. -/
theorem IsSidonMod.residue_injOn {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) :
    Set.InjOn (· % M_int) (↑S : Set ℤ) := by
  intro a ha b hb heq
  have ha' : a ∈ S := by exact ha
  have hb' : b ∈ S := by exact hb
  have hdvd : M_int ∣ (a - b) := by
    have ha_mod := Int.emod_def a M_int
    have hb_mod := Int.emod_def b M_int
    have heq' : a % M_int = b % M_int := heq
    rw [heq'] at ha_mod
    exact ⟨a / M_int - b / M_int, by linarith⟩
  have hS' : M_int ∣ ((a + a) - (b + a)) := by
    convert hdvd using 1; ring
  have hresult := hS (a := a) (b := a) (c := b) (d := a) ha' ha' hb' ha' hS'
  rcases hresult with h | h
  · exact h.1
  · exact h.2

/-- The residue map to natural representatives preserves cardinality on a
modular Sidon set. -/
theorem IsSidonMod.residueImageNat_card {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) :
    (residueImageNat M_int S).card = S.card := by
  unfold residueImageNat
  apply Finset.card_image_of_injOn
  intro a ha b hb hab
  have hcast : (Int.toNat (a % M_int) : ℤ) = Int.toNat (b % M_int) := by
    exact congrArg (fun n : ℕ => (n : ℤ)) hab
  have ha_nonneg : 0 ≤ a % M_int := Int.emod_nonneg _ (ne_of_gt hM)
  have hb_nonneg : 0 ≤ b % M_int := Int.emod_nonneg _ (ne_of_gt hM)
  have hmod : a % M_int = b % M_int := by
    simpa [Int.toNat_of_nonneg ha_nonneg, Int.toNat_of_nonneg hb_nonneg] using hcast
  exact hS.residue_injOn hM (by exact ha) (by exact hb) hmod

theorem residueImageNat_lt_modulus {M_int : ℤ} {S : Finset ℤ}
    (hM : 0 < M_int) {n : ℕ} (hn : n ∈ residueImageNat M_int S) :
    n < M_int.toNat := by
  rcases mem_residueImageNat.mp hn with ⟨s, _hs, hsmod⟩
  have hnonneg : 0 ≤ s % M_int := Int.emod_nonneg _ (ne_of_gt hM)
  have hlt : s % M_int < M_int := Int.emod_lt_of_pos _ hM
  have hcast : ((Int.toNat (s % M_int) : ℤ) : ℤ) < (M_int.toNat : ℤ) := by
    rw [Int.toNat_of_nonneg hnonneg, Int.toNat_of_nonneg (le_of_lt hM)]
    exact hlt
  have hnat : Int.toNat (s % M_int) < M_int.toNat := by
    exact_mod_cast hcast
  simpa [hsmod] using hnat

private theorem IsSidonMod.card_le_modulus {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) :
    S.card ≤ M_int.toNat := by
  calc
    S.card = (residueImageNat M_int S).card := by
      symm
      exact hS.residueImageNat_card hM
    _ ≤ (Finset.range M_int.toNat).card := by
      apply Finset.card_le_card
      intro n hn
      exact Finset.mem_range.mpr (residueImageNat_lt_modulus hM hn)
    _ = M_int.toNat := by
      simp

-- ============================================================
-- Part 2: Sorted residues and witnesses
-- ============================================================

/-- The sorted residues of a modular Sidon set, indexed by `Fin S.card`. -/
noncomputable def IsSidonMod.sortedResidues {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) :
    Fin S.card ↪o ℕ :=
  (residueImageNat M_int S).orderEmbOfFin (hS.residueImageNat_card hM)

theorem IsSidonMod.sortedResidues_mem {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) (i : Fin S.card) :
    hS.sortedResidues hM i ∈ residueImageNat M_int S := by
  simpa [IsSidonMod.sortedResidues] using
    (residueImageNat M_int S).orderEmbOfFin_mem
      (hS.residueImageNat_card hM) i

/-- A chosen element of `S` realizing the `i`th sorted residue. -/
noncomputable def IsSidonMod.sortedResidueWitness {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) (i : Fin S.card) : ℤ :=
  Classical.choose <|
    exists_mem_modEq_of_mem_residueImageNat hM <|
      hS.sortedResidues_mem hM i

theorem IsSidonMod.sortedResidueWitness_mem {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) (i : Fin S.card) :
    hS.sortedResidueWitness hM i ∈ S :=
  (Classical.choose_spec <|
    exists_mem_modEq_of_mem_residueImageNat hM <|
      hS.sortedResidues_mem hM i).1

theorem IsSidonMod.sortedResidueWitness_mod {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int) (i : Fin S.card) :
    ((hS.sortedResidues hM i : ℕ) : ℤ) = hS.sortedResidueWitness hM i % M_int :=
  (Classical.choose_spec <|
    exists_mem_modEq_of_mem_residueImageNat hM <|
      hS.sortedResidues_mem hM i).2

-- ============================================================
-- Part 3: Cyclic gap function
-- ============================================================

/-- Cyclic gap at position `i` for sorted residues `r : ℕ → ℕ` of a set
of `k` elements in `{0, …, M−1}`. For `i + 1 < k` this is the forward
step `r(i+1) − r(i)`; for the last index it is the wraparound step
`M − r(k−1) + r(0)`. -/
def cyclicGapAt (M k : ℕ) (r : ℕ → ℕ) (i : ℕ) : ℕ :=
  if i + 1 < k then r (i + 1) - r i
  else M - r (k - 1) + r 0

-- ============================================================
-- Part 4: Gap existence lemmas
-- ============================================================

private theorem no_mem_between_orderEmbOfFin {s : Finset ℕ} {k : ℕ}
    (h : s.card = k) {i : Fin k} (hi : (i : ℕ) + 1 < k) {n : ℕ} (hn : n ∈ s)
    (hleft : s.orderEmbOfFin h i < n)
    (hright : n < s.orderEmbOfFin h ⟨(i : ℕ) + 1, hi⟩) : False := by
  have hn_range : n ∈ Set.range (s.orderEmbOfFin h) := by
    rw [Finset.range_orderEmbOfFin]
    exact hn
  rcases hn_range with ⟨j, rfl⟩
  have hij_left : i < j := by
    by_contra hij_left
    exact not_lt_of_ge ((s.orderEmbOfFin h).monotone (le_of_not_gt hij_left)) hleft
  have hij_right : j < ⟨(i : ℕ) + 1, hi⟩ := by
    by_contra hij_right
    exact not_lt_of_ge ((s.orderEmbOfFin h).monotone (le_of_not_gt hij_right)) hright
  have hij_left' : (i : ℕ) + 1 ≤ (j : ℕ) := Nat.succ_le_of_lt hij_left
  have hij_right' : (j : ℕ) < (i : ℕ) + 1 := hij_right
  omega

private theorem IsSidonMod.no_residue_between_successive_sortedResidues
    {M_int : ℤ} {S : Finset ℤ} (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    {i : Fin S.card} (hi : (i : ℕ) + 1 < S.card) {n : ℕ}
    (hn : n ∈ residueImageNat M_int S)
    (hleft : hS.sortedResidues hM i < n)
    (hright : n < hS.sortedResidues hM ⟨(i : ℕ) + 1, hi⟩) : False := by
  exact no_mem_between_orderEmbOfFin (hS.residueImageNat_card hM) hi hn hleft hright

/-- **Difference distinctness.** In a modular Sidon set, if
`a − b ≡ c − d (mod M)` and `a ≠ b`, then `a = c` and `b = d`. -/
theorem IsSidonMod.diff_eq {M : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M S) {a b c d : ℤ}
    (ha : a ∈ S) (hb : b ∈ S) (hc : c ∈ S) (hd : d ∈ S)
    (hab : a ≠ b)
    (hdiff : M ∣ ((a - b) - (c - d))) :
    a = c ∧ b = d := by
  have hconv : M ∣ ((a + d) - (b + c)) := by convert hdiff using 1; ring
  rcases hS ha hd hb hc hconv with h | h
  · exact absurd h.1 hab
  · exact ⟨h.1, h.2.symm⟩

/-- If `M` divides the difference of residues of two pairs, then `M`
divides the difference of the original pairs. -/
theorem dvd_diff_of_residue_diff {M a b c d : ℤ}
    (h : M ∣ ((b % M - a % M) - (d % M - c % M))) :
    M ∣ ((b - a) - (d - c)) := by
  have key : (b - a) - (d - c) =
      M * (b / M - a / M - (d / M - c / M)) +
      ((b % M - a % M) - (d % M - c % M)) := by
    have hb := Int.emod_def b M
    have ha := Int.emod_def a M
    have hd := Int.emod_def d M
    have hc := Int.emod_def c M
    linarith
  rw [key]
  exact dvd_add (dvd_mul_right M _) h

-- ============================================================
-- Part 5: Window Sidon property (needed for gap lemmas)
-- ============================================================

/-- The cyclic-window relabeling map. For a window starting at `u` in
`ℤ/Mℤ`, sends representative `x` to position `(x − u) % M + 1`.
Elements in a window of length `N` land in `{1, …, N}`. -/
def windowRelabel (M u x : ℤ) : ℤ := (x - u) % M + 1

theorem windowRelabel_bounds {M u x N : ℤ} (hM : 0 < M)
    (hw : (x - u) % M < N) :
    1 ≤ windowRelabel M u x ∧ windowRelabel M u x ≤ N := by
  have hnn : 0 ≤ (x - u) % M := Int.emod_nonneg _ (ne_of_gt hM)
  exact ⟨by unfold windowRelabel; linarith,
         by unfold windowRelabel; linarith⟩

/-- The restriction of `S` to the cyclic window of length `N` starting at
`u`, relabeled into `{1, …, N}`. An element `s ∈ S` is kept iff
`(s − u) % M < N`. -/
def windowImage (M u N : ℤ) (S : Finset ℤ) : Finset ℤ :=
  (S.filter (fun s => (s - u) % M < N)).image (windowRelabel M u)

/-- In a modular Sidon set, `M ∣ (x − y)` forces `x = y`. -/
theorem IsSidonMod.eq_of_dvd {M : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M S) {x y : ℤ} (hx : x ∈ S) (hy : y ∈ S)
    (hdvd : M ∣ (x - y)) : x = y := by
  have hdvd' : M ∣ ((x + x) - (y + x)) := by
    convert hdvd using 1; ring
  rcases hS hx hx hy hx hdvd' with h | h
  · exact h.1
  · exact h.2

/-- Two elements of a modular Sidon set that relabel to the same value
must be equal. -/
theorem IsSidonMod.windowRelabel_injective {M : ℤ} {S : Finset ℤ} {u : ℤ}
    (hS : IsSidonMod M S)
    {x y : ℤ} (hx : x ∈ S) (hy : y ∈ S)
    (heq : windowRelabel M u x = windowRelabel M u y) : x = y := by
  unfold windowRelabel at heq
  have heq' : (x - u) % M = (y - u) % M := by linarith
  have hdvd : M ∣ (x - y) :=
    ⟨(x - u) / M - (y - u) / M, by
      have hx_mod := Int.emod_def (x - u) M
      have hy_mod := Int.emod_def (y - u) M
      linarith⟩
  exact hS.eq_of_dvd hx hy hdvd

/-- **Sum-transfer lemma.** If the relabeled sums of two pairs are equal
as integers, then the original sums are congruent modulo `M`. -/
theorem windowRelabel_sum_dvd {M u x y z w : ℤ}
    (hsum : windowRelabel M u x + windowRelabel M u y =
            windowRelabel M u z + windowRelabel M u w) :
    M ∣ ((x + y) - (z + w)) := by
  unfold windowRelabel at hsum
  have hx_mod := Int.emod_def (x - u) M
  have hy_mod := Int.emod_def (y - u) M
  have hz_mod := Int.emod_def (z - u) M
  have hw_mod := Int.emod_def (w - u) M
  exact ⟨(x - u) / M + (y - u) / M - (z - u) / M - (w - u) / M,
         by linarith⟩

/-- **Cyclic-window Sidon theorem.** The relabeled window restriction of
a modular Sidon set is an interval Sidon set. -/
theorem IsSidonMod.windowSidon {M : ℤ} {S : Finset ℤ} {u N : ℤ}
    (hS : IsSidonMod M S) (hM : 0 < M) (_hN : 0 < N) (_hNM : N ≤ M) :
    IsIntervalSidon N (windowImage M u N S) where
  subset := by
    intro a ha
    simp only [windowImage, mem_image, mem_filter] at ha
    obtain ⟨s, ⟨_, hsw⟩, rfl⟩ := ha
    exact windowRelabel_bounds hM (by simpa using hsw)
  sidon := by
    intro a b c d ha hb hc hd hsum
    simp only [windowImage, mem_image, mem_filter] at ha hb hc hd
    obtain ⟨x, ⟨hxS, _⟩, rfl⟩ := ha
    obtain ⟨y, ⟨hyS, _⟩, rfl⟩ := hb
    obtain ⟨z, ⟨hzS, _⟩, rfl⟩ := hc
    obtain ⟨w, ⟨hwS, _⟩, rfl⟩ := hd
    have hdvd := windowRelabel_sum_dvd hsum
    rcases hS hxS hyS hzS hwS hdvd with h | h
    · left
      exact ⟨congrArg (windowRelabel M u) h.1, congrArg (windowRelabel M u) h.2⟩
    · right
      exact ⟨congrArg (windowRelabel M u) h.1, congrArg (windowRelabel M u) h.2⟩

/-- The relabeled image has the same cardinality as the window
restriction, because `windowRelabel` is injective on any modular
Sidon set. -/
theorem IsSidonMod.windowImage_card {M : ℤ} {S : Finset ℤ} {u N : ℤ}
    (hS : IsSidonMod M S) :
    (windowImage M u N S).card =
      (S.filter (fun s => (s - u) % M < N)).card := by
  simp only [windowImage]
  apply Finset.card_image_of_injOn
  intro x hx y hy heq
  rw [Finset.mem_coe, Finset.mem_filter] at hx hy
  exact hS.windowRelabel_injective hx.1 hy.1 heq

-- ============================================================
-- Part 6: Gap existence from sorted residues
-- ============================================================

/-- **Zero-loss window existence.** If a modular Sidon set of size `k`
in `ℤ/Mℤ` has a cyclic gap of length at least `L`, then the window
of length `N = M − L` starting just past that gap captures all `k`
elements, giving an interval Sidon set of full size `k`. -/
theorem IsSidonMod.exists_full_intervalSidon_of_gap
    {M : ℤ} {S : Finset ℤ} {L : ℤ}
    (hS : IsSidonMod M S) (hM : 0 < M) (hL : 0 ≤ L) (hLM : L < M)
    (hgap : ∃ u : ℤ, ∀ s ∈ S, ¬((s - u) % M < L)) :
    ∃ A : Finset ℤ, IsIntervalSidon (M - L) A ∧ A.card = S.card := by
  obtain ⟨u, hu⟩ := hgap
  have hN_pos : 0 < M - L := by omega
  have hNM : M - L ≤ M := by omega
  refine ⟨windowImage M (u + L) (M - L) S,
         hS.windowSidon hM hN_pos hNM, ?_⟩
  have hcard_filter : (S.filter (fun s => (s - (u + L)) % M < M - L)).card = S.card := by
    have hall : ∀ s ∈ S, (s - (u + L)) % M < M - L := by
      intro s hs
      have hu_gap : ¬((s - u) % M < L) := hu s hs
      have hge : L ≤ (s - u) % M := not_lt.mp hu_gap
      have hmod_bound : (s - u) % M < M := Int.emod_lt_of_pos _ hM
      have hkey : (s - (u + L)) % M = ((s - u) % M - L) % M := by
        have := Int.emod_def (s - u) M
        have h_eq : s - (u + L) = (s - u) % M - L + M * ((s - u) / M) := by linarith
        rw [h_eq, Int.add_mul_emod_self_left]
      rw [hkey]
      have hsub_nonneg : 0 ≤ (s - u) % M - L := by omega
      have hsub_lt : (s - u) % M - L < M - L := by omega
      have hmod_result : ((s - u) % M - L) % M = (s - u) % M - L := by
        exact Int.emod_eq_of_lt hsub_nonneg (by omega : (s - u) % M - L < M)
      rw [hmod_result]
      omega
    have hsub : S.filter (fun s => (s - (u + L)) % M < M - L) = S := by
      ext s
      constructor
      · intro h; exact (Finset.mem_filter.mp h).1
      · intro hs; exact Finset.mem_filter.mpr ⟨hs, hall s hs⟩
    rw [hsub]
  rw [hS.windowImage_card, hcard_filter]

/-- A forward cyclic gap between consecutive sorted residues yields a full-size
interval Sidon set in the complementary interval. -/
theorem IsSidonMod.exists_full_intervalSidon_of_forward_sortedGap
    {M_int : ℤ} {S : Finset ℤ} (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    {i : Fin S.card} (hi : (i : ℕ) + 1 < S.card) :
    ∃ A : Finset ℤ,
      IsIntervalSidon
          (M_int
            - (hS.sortedResidues hM ⟨(i : ℕ) + 1, hi⟩ : ℕ)
            + hS.sortedResidues hM i + 1) A
        ∧ A.card = S.card := by
  classical
  let j : Fin S.card := ⟨(i : ℕ) + 1, hi⟩
  let aNat : ℕ := hS.sortedResidues hM i
  let bNat : ℕ := hS.sortedResidues hM j
  let a : ℤ := aNat
  let b : ℤ := bNat
  let L : ℤ := b - a - 1
  have hab : aNat < bNat := by
    simpa [aNat, bNat, j] using
      ((hS.sortedResidues hM).strictMono (show i < j by
        rw [Fin.lt_def]
        dsimp [j]
        omega))
  have hb_lt_M : bNat < M_int.toNat := by
    exact residueImageNat_lt_modulus hM (hS.sortedResidues_mem hM j)
  have hL_nonneg : 0 ≤ L := by
    dsimp [L, a, b]
    omega
  have hLM : L < M_int := by
    dsimp [L, a, b]
    omega
  have hgap : ∃ u : ℤ, ∀ s ∈ S, ¬((s - u) % M_int < L) := by
    refine ⟨a + 1, ?_⟩
    intro s hs hslt
    set n : ℕ := Int.toNat (s % M_int)
    have hsmod_nonneg : 0 ≤ s % M_int := Int.emod_nonneg _ (ne_of_gt hM)
    have hn_cast : (n : ℤ) = s % M_int := by
      dsimp [n]
      rw [Int.toNat_of_nonneg hsmod_nonneg]
    have hn_mem : n ∈ residueImageNat M_int S := by
      refine mem_residueImageNat.mpr ⟨s, hs, ?_⟩
      simp [n]
    have hn_lt_M : n < M_int.toNat := residueImageNat_lt_modulus hM hn_mem
    have hslt' : (((n : ℤ) - (a + 1)) % M_int) < L := by
      have hshift : (s - (a + 1)) % M_int = (((n : ℤ) - (a + 1)) % M_int) := by
        rw [hn_cast]
        have h_eq : s - (a + 1) = (s % M_int - (a + 1)) + M_int * (s / M_int) := by
          have hsdef := Int.emod_def s M_int
          linarith
        rw [h_eq, Int.add_mul_emod_self_left]
      rwa [hshift] at hslt
    by_cases hna : n < aNat + 1
    · have hmod_nonneg : 0 ≤ M_int + n - (a + 1) := by
        dsimp [a]
        omega
      have hmod_lt : M_int + n - (a + 1) < M_int := by
        dsimp [a]
        omega
      have hrewrite : (((n : ℤ) - (a + 1)) % M_int) = M_int + n - (a + 1) := by
        calc
          (((n : ℤ) - (a + 1)) % M_int)
            = ((((n : ℤ) - (a + 1)) + 1 * M_int) % M_int) := by
                simpa [sub_eq_add_neg, add_assoc, add_left_comm, add_comm] using
                  (Int.add_mul_emod_self_right ((n : ℤ) - (a + 1)) (1 : ℤ) M_int).symm
          _ = M_int + n - (a + 1) := by
                have hsum :
                    (((n : ℤ) - (a + 1)) + 1 * M_int) = M_int + n - (a + 1) := by
                  ring
                rw [hsum, Int.emod_eq_of_lt hmod_nonneg hmod_lt]
      have : M_int + n - (a + 1) < L := by
        simpa [hrewrite] using hslt'
      omega
    · have hna' : aNat + 1 ≤ n := le_of_not_gt hna
      have hmod_nonneg : 0 ≤ (n : ℤ) - (a + 1) := by
        dsimp [a]
        omega
      have hmod_lt : (n : ℤ) - (a + 1) < M_int := by
        dsimp [a]
        omega
      have hrewrite : (((n : ℤ) - (a + 1)) % M_int) = (n : ℤ) - (a + 1) := by
        rw [Int.emod_eq_of_lt hmod_nonneg hmod_lt]
      have hnb : n < bNat := by
        have : (n : ℤ) - (a + 1) < L := by
          simpa [hrewrite] using hslt'
        dsimp [L, a, b] at this
        omega
      have hna_lt : aNat < n := by
        omega
      exact hS.no_residue_between_successive_sortedResidues hM hi hn_mem hna_lt hnb
  have htransfer := hS.exists_full_intervalSidon_of_gap hM hL_nonneg hLM hgap
  have hlen :
      M_int - L =
        M_int - (hS.sortedResidues hM ⟨(i : ℕ) + 1, hi⟩ : ℕ) + hS.sortedResidues hM i + 1 := by
    dsimp [L, a, b, aNat, bNat, j]
    ring
  simpa [hlen] using htransfer

/-- The wraparound cyclic gap between the last and first sorted residues yields
a full-size interval Sidon set in the complementary interval. -/
theorem IsSidonMod.exists_full_intervalSidon_of_wraparound_sortedGap
    {M_int : ℤ} {S : Finset ℤ} (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    (hk : 2 ≤ S.card) :
    let first : Fin S.card := ⟨0, by omega⟩
    let last : Fin S.card := ⟨S.card - 1, Nat.sub_lt (by omega) (Nat.succ_pos 0)⟩
    ∃ A : Finset ℤ,
      IsIntervalSidon ((hS.sortedResidues hM last : ℕ) - hS.sortedResidues hM first + 1) A
        ∧ A.card = S.card := by
  classical
  have hcard_pos : 0 < S.card := by omega
  let firstIdx : Fin S.card := ⟨0, hcard_pos⟩
  let lastIdx : Fin S.card := ⟨S.card - 1, Nat.sub_lt hcard_pos (Nat.succ_pos 0)⟩
  let aNat : ℕ := hS.sortedResidues hM firstIdx
  let bNat : ℕ := hS.sortedResidues hM lastIdx
  let a : ℤ := aNat
  let b : ℤ := bNat
  let L : ℤ := M_int - b + a - 1
  have hb_lt_M : bNat < M_int.toNat := by
    exact residueImageNat_lt_modulus hM (hS.sortedResidues_mem hM lastIdx)
  have hab : aNat < bNat := by
    simpa [aNat, bNat, firstIdx, lastIdx] using
      ((hS.sortedResidues hM).strictMono (show firstIdx < lastIdx by
        rw [Fin.lt_def]
        dsimp [firstIdx, lastIdx]
        omega))
  have hL_nonneg : 0 ≤ L := by
    dsimp [L, a, b]
    omega
  have hLM : L < M_int := by
    dsimp [L, a, b]
    omega
  have hfirst_min :
      aNat = (residueImageNat M_int S).min' (Finset.card_pos.mp <| by
        simpa [hS.residueImageNat_card hM] using hcard_pos) := by
    simpa [aNat, firstIdx, IsSidonMod.sortedResidues] using
      (Finset.orderEmbOfFin_zero (s := residueImageNat M_int S)
        (h := hS.residueImageNat_card hM) hcard_pos)
  have hlast_max :
      bNat = (residueImageNat M_int S).max' (Finset.card_pos.mp <| by
        simpa [hS.residueImageNat_card hM] using hcard_pos) := by
    simpa [bNat, lastIdx, IsSidonMod.sortedResidues] using
      (Finset.orderEmbOfFin_last (s := residueImageNat M_int S)
        (h := hS.residueImageNat_card hM) hcard_pos)
  have hgap : ∃ u : ℤ, ∀ s ∈ S, ¬((s - u) % M_int < L) := by
    refine ⟨b + 1, ?_⟩
    intro s hs hslt
    set n : ℕ := Int.toNat (s % M_int)
    have hsmod_nonneg : 0 ≤ s % M_int := Int.emod_nonneg _ (ne_of_gt hM)
    have hn_cast : (n : ℤ) = s % M_int := by
      dsimp [n]
      rw [Int.toNat_of_nonneg hsmod_nonneg]
    have hn_mem : n ∈ residueImageNat M_int S := by
      refine mem_residueImageNat.mpr ⟨s, hs, ?_⟩
      simp [n]
    have hmin_le : aNat ≤ n := by
      rw [hfirst_min]
      exact (residueImageNat M_int S).min'_le _ hn_mem
    have hmax_ge : n ≤ bNat := by
      rw [hlast_max]
      exact (residueImageNat M_int S).le_max' n hn_mem
    have hslt' : (((n : ℤ) - (b + 1)) % M_int) < L := by
      have hshift : (s - (b + 1)) % M_int = (((n : ℤ) - (b + 1)) % M_int) := by
        rw [hn_cast]
        have h_eq : s - (b + 1) = (s % M_int - (b + 1)) + M_int * (s / M_int) := by
          have hsdef := Int.emod_def s M_int
          linarith
        rw [h_eq, Int.add_mul_emod_self_left]
      rwa [hshift] at hslt
    by_cases hbn : bNat < n
    · exact absurd hbn (not_lt_of_ge hmax_ge)
    · have hnb : n ≤ bNat := le_of_not_gt hbn
      have hmod_nonneg : 0 ≤ M_int + n - (b + 1) := by
        dsimp [b]
        omega
      have hmod_lt : M_int + n - (b + 1) < M_int := by
        dsimp [b]
        omega
      have hrewrite : (((n : ℤ) - (b + 1)) % M_int) = M_int + n - (b + 1) := by
        calc
          (((n : ℤ) - (b + 1)) % M_int)
            = ((((n : ℤ) - (b + 1)) + 1 * M_int) % M_int) := by
                simpa [sub_eq_add_neg, add_assoc, add_left_comm, add_comm] using
                  (Int.add_mul_emod_self_right ((n : ℤ) - (b + 1)) (1 : ℤ) M_int).symm
          _ = M_int + n - (b + 1) := by
                have hsum :
                    (((n : ℤ) - (b + 1)) + 1 * M_int) = M_int + n - (b + 1) := by
                  ring
                rw [hsum, Int.emod_eq_of_lt hmod_nonneg hmod_lt]
      have hna : n < aNat := by
        have : M_int + n - (b + 1) < L := by
          simpa [hrewrite] using hslt'
        dsimp [L, a, b] at this
        omega
      exact absurd hna (not_lt_of_ge hmin_le)
  have htransfer := hS.exists_full_intervalSidon_of_gap hM hL_nonneg hLM hgap
  have hlen :
      M_int - L = b - a + 1 := by
    dsimp [L]
    ring
  simpa [firstIdx, lastIdx, aNat, bNat, a, b, hlen] using htransfer

/-- Any cyclic gap in the sorted residue model yields a full-size interval
Sidon set in the complementary interval of length `M - gap + 1`. -/
theorem IsSidonMod.exists_full_intervalSidon_of_cyclicGapAt
    {M_int : ℤ} {S : Finset ℤ} (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    (hk : 2 ≤ S.card) (i : Fin S.card) :
    let r : ℕ → ℕ := fun n =>
      if hn : n < S.card then hS.sortedResidues hM ⟨n, hn⟩ else 0
    ∃ A : Finset ℤ,
      IsIntervalSidon (M_int - cyclicGapAt M_int.toNat S.card r i + 1) A ∧ A.card = S.card := by
  classical
  let r : ℕ → ℕ := fun n =>
    if hn : n < S.card then hS.sortedResidues hM ⟨n, hn⟩ else 0
  by_cases hi_wrap : (i : ℕ) + 1 < S.card
  · let j : Fin S.card := ⟨(i : ℕ) + 1, hi_wrap⟩
    have hforward := hS.exists_full_intervalSidon_of_forward_sortedGap hM (i := i) hi_wrap
    have hij : i < j := by
      rw [Fin.lt_def]
      dsimp [j]
      omega
    have hri_le : hS.sortedResidues hM i ≤ hS.sortedResidues hM j :=
      le_of_lt ((hS.sortedResidues hM).strictMono hij)
    have hcyc :
        (cyclicGapAt M_int.toNat S.card r i : ℤ) =
          (hS.sortedResidues hM j : ℕ) - hS.sortedResidues hM i := by
      simp [cyclicGapAt, r, hi_wrap, i.2, j, Nat.cast_sub hri_le]
    have hlen :
        M_int - (cyclicGapAt M_int.toNat S.card r i : ℕ) + 1 =
          M_int - (hS.sortedResidues hM j : ℕ) + hS.sortedResidues hM i + 1 := by
      rw [hcyc]
      ring
    simpa [r, j, hlen] using hforward
  · have hcard_pos : 0 < S.card := by omega
    let firstIdx : Fin S.card := ⟨0, hcard_pos⟩
    let lastIdx : Fin S.card := ⟨S.card - 1, Nat.sub_lt hcard_pos (Nat.succ_pos 0)⟩
    have hi_last : i = lastIdx := by
      apply Fin.ext
      dsimp [lastIdx]
      omega
    have hwrap := hS.exists_full_intervalSidon_of_wraparound_sortedGap hM hk
    have hb_lt_M :
        hS.sortedResidues hM lastIdx < M_int.toNat := by
      exact residueImageNat_lt_modulus hM (hS.sortedResidues_mem hM lastIdx)
    have hlast_le : hS.sortedResidues hM lastIdx ≤ M_int.toNat := le_of_lt hb_lt_M
    have hM_cast : (M_int.toNat : ℤ) = M_int := by
      rw [Int.toNat_of_nonneg (le_of_lt hM)]
    have hcyc :
        (cyclicGapAt M_int.toNat S.card r i : ℤ) =
          M_int - hS.sortedResidues hM lastIdx + hS.sortedResidues hM firstIdx := by
      rw [hi_last]
      have hnot_last_succ : ¬((S.card - 1 : ℕ) + 1 < S.card) := by omega
      simp [cyclicGapAt, r, hnot_last_succ, hcard_pos, firstIdx, lastIdx,
        Nat.cast_sub hlast_le, hM_cast]
    have hlen :
        M_int - (cyclicGapAt M_int.toNat S.card r i : ℕ) + 1 =
          (hS.sortedResidues hM lastIdx : ℕ) - hS.sortedResidues hM firstIdx + 1 := by
      rw [hcyc]
      ring
    simpa [r, firstIdx, lastIdx, hlen] using hwrap

-- ============================================================
-- Part 7: Gap-distinctness structure
-- ============================================================

/-- A finset of positive naturals summing to `M` models the cyclic gap
structure of a Sidon set. -/
structure DistinctPosPartsOf (M : ℕ) (k : ℕ) where
  parts : Finset ℕ
  card_eq : parts.card = k
  pos : 0 ∉ parts
  sum_eq : ∑ x ∈ parts, x = M

theorem DistinctPosPartsOf.nonempty {M k : ℕ}
    (g : DistinctPosPartsOf M k) (hk : 0 < k) : g.parts.Nonempty := by
  rw [Finset.nonempty_iff_ne_empty]
  intro h
  have h1 : g.parts.card = 0 := by rw [h, Finset.card_empty]
  have h2 := g.card_eq
  omega

/-- The sum of the elements of any `Finset ℕ` is at least the
`card`-th triangular number `card * (card − 1) / 2`. -/
theorem sum_finset_nat_ge_tri (s : Finset ℕ) :
    s.card * (s.card - 1) / 2 ≤ ∑ x ∈ s, x := by
  suffices h : s.card * (s.card - 1) ≤ 2 * ∑ x ∈ s, x by omega
  have key : ∀ n, ∀ t : Finset ℕ, t.card = n →
      n * (n - 1) ≤ 2 * ∑ x ∈ t, x := by
    intro n
    induction n using Nat.strongRecOn with
    | ind n ih =>
      intro t ht
      cases n with
      | zero => simp [Finset.card_eq_zero.mp ht]
      | succ n =>
        have hne : t.Nonempty := Finset.card_pos.mp (by omega)
        set m := t.max' hne
        have hm_mem : m ∈ t := Finset.max'_mem t hne
        set t' := t.erase m
        have ht'_card : t'.card = n := by
          rw [Finset.card_erase_of_mem hm_mem]; omega
        have hm_ge : n ≤ m := by
          have hle : t.card ≤ t.max' hne + 1 := by
            calc t.card
                ≤ (Finset.range (t.max' hne + 1)).card := by
                  apply Finset.card_le_card
                  intro x hx
                  exact Finset.mem_range.mpr (Nat.lt_succ_of_le (t.le_max' x hx))
              _ = t.max' hne + 1 := Finset.card_range _
          omega
        have ih_t' : n * (n - 1) ≤ 2 * ∑ x ∈ t', x :=
          ih n (by omega) t' ht'_card
        have hsum : ∑ x ∈ t, x = ∑ x ∈ t', x + m :=
          (Finset.sum_erase_add t (fun x => x) hm_mem).symm
        rw [hsum]
        simp only [Nat.succ_sub_one]
        have hrec : n * (n - 1) + 2 * n = (n + 1) * n := by
          cases n with
          | zero => simp
          | succ n => simp only [Nat.succ_sub_one]; ring
        linarith
  exact key s.card s rfl

/-- **Complement-max bound.** For a nonempty `Finset ℕ` with `0 ∉ s`,
the sum plus `card*(card − 1)/2` is at most `card * max`. -/
theorem sum_add_tri_le_card_mul_max (s : Finset ℕ) (hne : s.Nonempty) (_hpos : 0 ∉ s) :
    (∑ x ∈ s, x) + s.card * (s.card - 1) / 2 ≤ s.card * s.max' hne := by
  set m := s.max' hne
  have hinj : Set.InjOn (fun x => m - x) ↑s := by
    intro a ha b hb hab
    have ha' : a ≤ m := s.le_max' a (show a ∈ s from ha)
    have hb' : b ≤ m := s.le_max' b (show b ∈ s from hb)
    have h1 : m - a + a = m := Nat.sub_add_cancel ha'
    have h2 : m - b + b = m := Nat.sub_add_cancel hb'
    have hab' : m - a = m - b := hab
    omega
  set t := s.image (fun x => m - x)
  have ht_card : t.card = s.card := Finset.card_image_of_injOn hinj
  have hsum_rel : ∑ y ∈ t, y + ∑ x ∈ s, x = s.card * m := by
    show ∑ y ∈ s.image (fun x => m - x), y + ∑ x ∈ s, x = s.card * m
    rw [Finset.sum_image hinj, ← Finset.sum_add_distrib]
    rw [show s.card * m = ∑ _ ∈ s, m from
      (Finset.sum_const_nat (fun _ _ => rfl)).symm]
    apply Finset.sum_congr rfl
    intro x hx
    exact Nat.sub_add_cancel (s.le_max' x hx)
  have htri := sum_finset_nat_ge_tri t
  rw [ht_card] at htri
  linarith

/-- **Gap lower bound.** If `M` is partitioned into `k ≥ 1` distinct
positive parts, the largest part is at least `(M + k(k−1)/2) / k`. -/
theorem max_gap_lower_bound {M k : ℕ} (g : DistinctPosPartsOf M k)
    (hk : 0 < k) :
    (M + k * (k - 1) / 2) / k ≤ g.parts.max' (g.nonempty hk) := by
  have hbound := sum_add_tri_le_card_mul_max g.parts (g.nonempty hk) g.pos
  rw [g.card_eq, g.sum_eq] at hbound
  exact Nat.div_le_of_le_mul hbound

-- ============================================================
-- Part 8: Sidon gaps to DistinctPosPartsOf
-- ============================================================

private theorem mono_le_of_step (n : ℕ) (f : ℕ → ℕ)
    (hf : ∀ i, i < n → f i ≤ f (i + 1)) : f 0 ≤ f n := by
  induction n with
  | zero => omega
  | succ n ih =>
    exact le_trans (ih (fun i hi => hf i (by omega))) (hf n (by omega))

private theorem nat_telescope (n : ℕ) (f : ℕ → ℕ)
    (hf : ∀ i, i < n → f i ≤ f (i + 1)) :
    (∑ i ∈ Finset.range n, (f (i + 1) - f i)) + f 0 = f n := by
  induction n with
  | zero => simp
  | succ n ih =>
    rw [Finset.sum_range_succ]
    have h0n := mono_le_of_step n f (fun i hi => hf i (by omega))
    have hnn1 : f n ≤ f (n + 1) := hf n (by omega)
    have := ih (fun i hi => hf i (by omega))
    omega

private theorem mono_le_of_adjacent {k : ℕ} {r : ℕ → ℕ}
    (hr_mono : ∀ i, i + 1 < k → r i < r (i + 1))
    {a b : ℕ} (hab : a ≤ b) (hb : b < k) :
    r a ≤ r b := by
  let g : ℕ → ℕ := fun t => r (a + t)
  have hstep : ∀ i, i < b - a → g i ≤ g (i + 1) := by
    intro i hi
    dsimp [g]
    exact le_of_lt (hr_mono (a + i) (by omega))
  have hmono := mono_le_of_step (b - a) g hstep
  dsimp [g] at hmono
  simpa [Nat.add_sub_of_le hab] using hmono

theorem cyclicGapAt_pos {M k : ℕ} {r : ℕ → ℕ} (_hk : 2 ≤ k)
    (hMono : ∀ i, i + 1 < k → r i < r (i + 1))
    (hBound : r (k - 1) < M)
    {i : ℕ} (_hi : i < k) :
    0 < cyclicGapAt M k r i := by
  simp only [cyclicGapAt]
  split
  · have := hMono i (by omega); omega
  · omega

theorem cyclicGapAt_sum {M k : ℕ} {r : ℕ → ℕ} (hk : 2 ≤ k)
    (hMono : ∀ i, i + 1 < k → r i < r (i + 1))
    (hBound : r (k - 1) < M) :
    ∑ i ∈ Finset.range k, cyclicGapAt M k r i = M := by
  have hsplit : ∑ i ∈ Finset.range k, cyclicGapAt M k r i =
      (∑ i ∈ Finset.range (k - 1), cyclicGapAt M k r i) +
        cyclicGapAt M k r (k - 1) := by
    have h := Finset.sum_range_succ (cyclicGapAt M k r) (k - 1)
    rwa [show k - 1 + 1 = k from by omega] at h
  rw [hsplit]
  have hlast : cyclicGapAt M k r (k - 1) = M - r (k - 1) + r 0 := by
    simp only [cyclicGapAt]; split <;> omega
  rw [hlast]
  have hfirst : ∀ i ∈ Finset.range (k - 1),
      cyclicGapAt M k r i = r (i + 1) - r i := by
    intro i hi
    simp only [cyclicGapAt]
    have him : i < k - 1 := Finset.mem_range.mp hi
    split
    · rfl
    · omega
  rw [Finset.sum_congr rfl hfirst]
  have hle : ∀ i, i < k - 1 → r i ≤ r (i + 1) :=
    fun i hi => le_of_lt (hMono i (by omega))
  have := nat_telescope (k - 1) r hle
  have := mono_le_of_step (k - 1) r hle
  omega

/-- **Gap injectivity.** If a sorted sequence of residues comes from a
modular Sidon set, then the cyclic gap function is injective. -/
theorem cyclicGapAt_injective_of_sidon
    {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    {k : ℕ} (hk : 2 ≤ k)
    (f : ℕ → ℤ)
    (hf_mem : ∀ i, i < k → f i ∈ S)
    (hf_inj : ∀ i j, i < k → j < k → f i = f j → i = j)
    (r : ℕ → ℕ)
    (hr_def : ∀ i, i < k → (r i : ℤ) = f i % M_int)
    (hr_mono : ∀ i, i + 1 < k → r i < r (i + 1))
    (hr_bound : r (k - 1) < M_int.toNat) :
    ∀ i j, i < k → j < k →
      cyclicGapAt M_int.toNat k r i = cyclicGapAt M_int.toNat k r j → i = j := by
  set M := M_int.toNat
  intro i j hi hj heq
  have hM_cast : M_int = (M : ℤ) := (Int.toNat_of_nonneg (le_of_lt hM)).symm
  have hf_ne : ∀ a b, a < k → b < k → a ≠ b → f a ≠ f b :=
    fun a b ha hb hab hf => absurd (hf_inj a b ha hb hf) hab
  by_cases hi_wrap : i + 1 < k <;> by_cases hj_wrap : j + 1 < k
  · simp only [cyclicGapAt, hi_wrap, hj_wrap, ↓reduceIte] at heq
    have hri_le : r i ≤ r (i + 1) := le_of_lt (hr_mono i hi_wrap)
    have hrj_le : r j ≤ r (j + 1) := le_of_lt (hr_mono j hj_wrap)
    have h_resid : f (i + 1) % M_int - f i % M_int =
        f (j + 1) % M_int - f j % M_int := by
      rw [← hr_def i (by omega), ← hr_def (i + 1) (by omega),
          ← hr_def j (by omega), ← hr_def (j + 1) (by omega)]
      have h1 : (↑(r (i + 1) - r i) : ℤ) = ↑(r (i + 1)) - ↑(r i) :=
        Nat.cast_sub hri_le
      have h2 : (↑(r (j + 1) - r j) : ℤ) = ↑(r (j + 1)) - ↑(r j) :=
        Nat.cast_sub hrj_le
      have h3 : (r (i + 1) - r i : ℕ) = r (j + 1) - r j := heq
      have h3_cast : (↑(r (i + 1) - r i) : ℤ) = ↑(r (j + 1) - r j) := by exact_mod_cast h3
      linarith [h3_cast]
    have hdvd : M_int ∣ ((f (i + 1) - f i) - (f (j + 1) - f j)) :=
      dvd_diff_of_residue_diff (by rw [h_resid, sub_self]; exact dvd_zero _)
    have hab : f i ≠ f (i + 1) := hf_ne i (i + 1) (by omega) (by omega) (by omega)
    have ⟨_, h2⟩ := hS.diff_eq (hf_mem (i + 1) (by omega)) (hf_mem i (by omega))
      (hf_mem (j + 1) (by omega)) (hf_mem j (by omega))
      (Ne.symm hab) hdvd
    exact hf_inj i j (by omega) (by omega) h2
  · have hj_eq : j = k - 1 := by omega
    subst hj_eq
    simp only [cyclicGapAt, hi_wrap, show ¬(k - 1 + 1 < k) from by omega,
      ↓reduceIte] at heq
    have hri_le : r i ≤ r (i + 1) := le_of_lt (hr_mono i hi_wrap)
    have hrk_le : r (k - 1) ≤ M := le_of_lt hr_bound
    have h_resid : (f (i + 1) % M_int - f i % M_int) -
        (f 0 % M_int - f (k - 1) % M_int) = M_int := by
      rw [← hr_def i (by omega), ← hr_def (i + 1) (by omega),
          ← hr_def 0 (by omega), ← hr_def (k - 1) (by omega)]
      have h1 : (↑(r (i + 1) - r i) : ℤ) = ↑(r (i + 1)) - ↑(r i) :=
        Nat.cast_sub hri_le
      have h2 : (↑(M - r (k - 1)) : ℤ) = ↑M - ↑(r (k - 1)) :=
        Nat.cast_sub hrk_le
      have heq_cast : (↑(r (i + 1) - r i) : ℤ) = ↑(M - r (k - 1) + r 0) := by exact_mod_cast heq
      linarith [heq_cast, hM_cast]
    have hdvd : M_int ∣ ((f (i + 1) - f i) - (f 0 - f (k - 1))) :=
      dvd_diff_of_residue_diff ⟨1, by linarith [h_resid]⟩
    have hab : f i ≠ f (i + 1) := hf_ne i (i + 1) (by omega) (by omega) (by omega)
    have ⟨h1, _⟩ := hS.diff_eq (hf_mem (i + 1) (by omega)) (hf_mem i (by omega))
      (hf_mem 0 (by omega)) (hf_mem (k - 1) (by omega))
      (Ne.symm hab) hdvd
    exact absurd (hf_inj (i + 1) 0 (by omega) (by omega) h1) (by omega)
  · have hi_eq : i = k - 1 := by omega
    subst hi_eq
    simp only [cyclicGapAt, show ¬(k - 1 + 1 < k) from by omega, hj_wrap,
      ↓reduceIte] at heq
    have hrj_le : r j ≤ r (j + 1) := le_of_lt (hr_mono j hj_wrap)
    have hrk_le : r (k - 1) ≤ M := le_of_lt hr_bound
    have h_resid : (f 0 % M_int - f (k - 1) % M_int) -
        (f (j + 1) % M_int - f j % M_int) = -M_int := by
      rw [← hr_def 0 (by omega), ← hr_def (k - 1) (by omega),
          ← hr_def j (by omega), ← hr_def (j + 1) (by omega)]
      have h1 : (↑(r (j + 1) - r j) : ℤ) = ↑(r (j + 1)) - ↑(r j) :=
        Nat.cast_sub hrj_le
      have h2 : (↑(M - r (k - 1)) : ℤ) = ↑M - ↑(r (k - 1)) :=
        Nat.cast_sub hrk_le
      have heq_cast : (↑(M - r (k - 1) + r 0) : ℤ) = ↑(r (j + 1) - r j) := by exact_mod_cast heq
      linarith [heq_cast, hM_cast]
    have hdvd : M_int ∣ ((f 0 - f (k - 1)) - (f (j + 1) - f j)) :=
      dvd_diff_of_residue_diff ⟨-1, by linarith [h_resid]⟩
    have hab : f (k - 1) ≠ f 0 :=
      hf_ne (k - 1) 0 (by omega) (by omega) (by omega)
    have ⟨h1, _⟩ := hS.diff_eq (hf_mem 0 (by omega)) (hf_mem (k - 1) (by omega))
      (hf_mem (j + 1) (by omega)) (hf_mem j (by omega))
      (Ne.symm hab) hdvd
    exact absurd (hf_inj 0 (j + 1) (by omega) (by omega) h1) (by omega)
  · omega

/-- **Sidon gap-distinctness bridge.** Given a modular Sidon set `S` of
size `k ≥ 2` in `ℤ/Mℤ`, together with an explicit sorted enumeration of
its residues, the cyclic gaps form a `DistinctPosPartsOf M k`. -/
noncomputable def sidon_gaps_to_distinctPosPartsOf
    {M_int : ℤ} {S : Finset ℤ}
    (hS : IsSidonMod M_int S) (hM : 0 < M_int)
    {k : ℕ} (hk : 2 ≤ k) (_hcard : S.card = k)
    (f : ℕ → ℤ)
    (hf_mem : ∀ i, i < k → f i ∈ S)
    (hf_inj : ∀ i j, i < k → j < k → f i = f j → i = j)
    (r : ℕ → ℕ)
    (hr_def : ∀ i, i < k → (r i : ℤ) = f i % M_int)
    (hr_mono : ∀ i, i + 1 < k → r i < r (i + 1))
    (hr_bound : r (k - 1) < M_int.toNat) :
    DistinctPosPartsOf M_int.toNat k where
  parts := (Finset.range k).image (cyclicGapAt M_int.toNat k r)
  card_eq := by
    apply Eq.trans (Finset.card_image_of_injOn _) (Finset.card_range k)
    intro i hi j hj heq
    have hi' : i < k := Finset.mem_range.mp (Finset.mem_coe.mp hi)
    have hj' : j < k := Finset.mem_range.mp (Finset.mem_coe.mp hj)
    exact cyclicGapAt_injective_of_sidon hS hM hk f hf_mem hf_inj
      r hr_def hr_mono hr_bound i j hi' hj' heq
  pos := by
    intro h0
    rw [Finset.mem_image] at h0
    obtain ⟨i, hi, hgap⟩ := h0
    rw [Finset.mem_range] at hi
    have := cyclicGapAt_pos hk hr_mono hr_bound hi
    omega
  sum_eq := by
    rw [Finset.sum_image]
    · exact cyclicGapAt_sum hk hr_mono hr_bound
    · intro i hi j hj heq
      exact cyclicGapAt_injective_of_sidon hS hM hk f hf_mem hf_inj
        r hr_def hr_mono hr_bound i j (Finset.mem_range.mp hi) (Finset.mem_range.mp hj) heq

-- ============================================================
-- Part 9: Main theorem — exists_full_intervalSidon_of_quantitative_gap_bound
-- ============================================================


/-- Residue map s % M is injective on an IsSidonMod M set. -/
theorem IsSidonMod.residue_inj_on {M : ℤ} {S : Finset ℤ} (hS : IsSidonMod M S) (hM : M ≠ 0) (x y : ℤ) (hx : x ∈ S) (hy : y ∈ S) (hres : x % M = y % M) : x = y := by
  by_contra hne
  have hsub : x - y = M * (x / M - y / M) := by
    nlinarith [Int.ediv_add_emod x M, Int.ediv_add_emod y M, hres]
  have hmod : M ∣ (x + x) - (x + y) := by
    have : (x + x) - (x + y) = x - y := by ring
    rw [this, hsub]
    exact ⟨x / M - y / M, by ring⟩
  have hS' := hS hx hx hx hy hmod
  rcases hS' with ⟨hac, hbd⟩ | ⟨had, hbc⟩
  · exact hne hbd
  · exact hne had

/-- Singer set via residues: A = {s % M + 1 : s ∈ S} is IsIntervalSidon N, |A| = p+1. -/
theorem singerIntervalSidon (p N : ℕ) (hp : Nat.Prime p) (hN : p * p + p + 1 ≤ N) :
    ∃ A : Finset ℤ, IsIntervalSidon (N : ℤ) A ∧ A.card = p + 1 := by
  rcases singerFamilyHypothesis_holds p hp with ⟨S, hS, hcard⟩
  set M := (p * p + p + 1 : ℤ) with hM_def
  have hp_pos : 0 < p := hp.pos
  have hM_pos : M ≠ 0 := by positivity
  have hM_N : M ≤ (N : ℤ) := by
    simpa [hM_def] using mod_cast hN
  have hinj_on : ∀ (x y : ℤ), x ∈ S → y ∈ S → x % M = y % M → x = y :=
    hS.residue_inj_on hM_pos
  let A := Finset.image (fun (s : ℤ) => s % M + 1) S
  have hA_card : A.card = p + 1 := by
    have hinj_on_S : Set.InjOn (fun (s : ℤ) => s % M + 1) (S : Set ℤ) := by
      intro x hxS y hyS h
      apply hinj_on x y hxS hyS
      -- h : (x % M + 1) = (y % M + 1), so x % M = y % M
      linarith
    calc
      A.card = S.card := Finset.card_image_of_injOn hinj_on_S
      _ = p + 1 := hcard
  have hA_sub : ∀ a ∈ A, 1 ≤ a ∧ a ≤ (N : ℤ) := by
    intro a ha
    rcases Finset.mem_image.mp ha with ⟨s, hs, rfl⟩
    have h_nonneg : 0 ≤ s % M := Int.emod_nonneg s (by intro h; exact hM_pos (h.symm ▸ rfl))
    have h_lt : s % M < M := Int.emod_lt s hM_pos
    have h_bound : s % M + 1 ≤ (N : ℤ) := by
      have : s % M + 1 ≤ M := by omega
      omega
    exact ⟨by omega, h_bound⟩
  have hA_sidon : IsSidon A := by
    intro a b c d ha hb hc hd hsum
    rcases Finset.mem_image.mp ha with ⟨s₁, hs₁, rfl⟩
    rcases Finset.mem_image.mp hb with ⟨s₂, hs₂, rfl⟩
    rcases Finset.mem_image.mp hc with ⟨s₃, hs₃, rfl⟩
    rcases Finset.mem_image.mp hd with ⟨s₄, hs₄, rfl⟩
    have hsum_mod : s₁ % M + s₂ % M = s₃ % M + s₄ % M := by
      linarith
    have hM_dvd_each : ∀ (s : ℤ), M ∣ s - s % M := by
      intro s
      have : s - s % M = M * (s / M) := by
        nlinarith [Int.ediv_add_emod s M]
      rw [this]
      exact ⟨s / M, by ring⟩
    have hM_dvd : M ∣ (s₁ + s₂) - (s₃ + s₄) := by
      have h_sub : M ∣ (s₁ - s₁ % M) + (s₂ - s₂ % M) - (s₃ - s₃ % M) - (s₄ - s₄ % M) := by
        have h_sum : M ∣ (s₁ - s₁ % M) + (s₂ - s₂ % M) := dvd_add (hM_dvd_each s₁) (hM_dvd_each s₂)
        have h_sum' : M ∣ (s₃ - s₃ % M) + (s₄ - s₄ % M) := dvd_add (hM_dvd_each s₃) (hM_dvd_each s₄)
        have h_sub' : M ∣ (s₁ - s₁ % M) + (s₂ - s₂ % M) - ((s₃ - s₃ % M) + (s₄ - s₄ % M)) :=
          dvd_sub h_sum h_sum'
        simpa [sub_sub] using h_sub'
      have h_mod_zero : s₁ % M + s₂ % M - s₃ % M - s₄ % M = 0 := by
        rw [hsum_mod]; ring
      have h_eq : (s₁ + s₂) - (s₃ + s₄) = ((s₁ - s₁ % M) + (s₂ - s₂ % M) - (s₃ - s₃ % M) - (s₄ - s₄ % M)) := by
        calc
          (s₁ + s₂) - (s₃ + s₄) = ((s₁ - s₁ % M) + (s₂ - s₂ % M) - (s₃ - s₃ % M) - (s₄ - s₄ % M))
            + (s₁ % M + s₂ % M - s₃ % M - s₄ % M) := by ring
          _ = ((s₁ - s₁ % M) + (s₂ - s₂ % M) - (s₃ - s₃ % M) - (s₄ - s₄ % M)) + 0 := by rw [h_mod_zero]
          _ = ((s₁ - s₁ % M) + (s₂ - s₂ % M) - (s₃ - s₃ % M) - (s₄ - s₄ % M)) := by ring
      rw [h_eq]
      exact h_sub
    have hS' := hS hs₁ hs₂ hs₃ hs₄ hM_dvd
    rcases hS' with (⟨h₁, h₂⟩ | ⟨h₁, h₂⟩)
    · left; constructor
      · apply congrArg (fun x : ℤ => x % M + 1) h₁
      · apply congrArg (fun x : ℤ => x % M + 1) h₂
    · right; constructor
      · apply congrArg (fun x : ℤ => x % M + 1) h₁
      · apply congrArg (fun x : ℤ => x % M + 1) h₂
  exact ⟨A, ⟨hA_sub, hA_sidon⟩, hA_card⟩

/-! ## Conditional Erdős Problem 30 -/

/-- Helper: `Nat.sqrt N` is a lower bound for the real square root. -/
private lemma natSqrt_le_real_sqrt (N : ℕ) :
    (Nat.sqrt N : ℝ) ≤ Real.sqrt (N : ℝ) := by
  have hs : Nat.sqrt N * Nat.sqrt N ≤ N := Nat.le_sqrt.1 (le_refl (Nat.sqrt N))
  calc
    (Nat.sqrt N : ℝ) = Real.sqrt ((Nat.sqrt N : ℝ) * (Nat.sqrt N : ℝ)) :=
      (Real.sqrt_mul_self (Nat.cast_nonneg _)).symm
    _ ≤ Real.sqrt (N : ℝ) := Real.sqrt_le_sqrt (by exact_mod_cast hs)

/-- Helper: the real square root lies strictly below `Nat.sqrt N + 1`. -/
private lemma real_sqrt_lt_natSqrt_add_one (N : ℕ) :
    Real.sqrt (N : ℝ) < (Nat.sqrt N : ℝ) + 1 := by
  have hlt : N < (Nat.sqrt N + 1) * (Nat.sqrt N + 1) := Nat.lt_succ_sqrt N
  have hpos : (0 : ℝ) ≤ (Nat.sqrt N : ℝ) + 1 := by positivity
  calc
    Real.sqrt (N : ℝ)
        < Real.sqrt (((Nat.sqrt N : ℝ) + 1) * ((Nat.sqrt N : ℝ) + 1)) :=
          Real.sqrt_lt_sqrt (Nat.cast_nonneg N) (by exact_mod_cast hlt)
    _ = (Nat.sqrt N : ℝ) + 1 := Real.sqrt_mul_self hpos

/-- Helper: a prime `p` with `p + 2 ≤ Nat.sqrt N` has its Singer modulus
    `p² + p + 1` below `N`, so the Singer set fits inside `{1, …, N}`. -/
private lemma singer_modulus_fits {p N : ℕ}
    (hpN : p + 2 ≤ Nat.sqrt N) : p * p + p + 1 ≤ N := by
  have hs : Nat.sqrt N * Nat.sqrt N ≤ N := Nat.le_sqrt.1 (le_refl (Nat.sqrt N))
  have h2 : (p + 2) * (p + 2) ≤ Nat.sqrt N * Nat.sqrt N := Nat.mul_le_mul hpN hpN
  nlinarith

/-- **Conditional Erdős Problem 30.** A subpolynomial prime-gap hypothesis
    around `√N`, together with an upper-bound hypothesis
    `h(N) ≤ √N + C·N^ε`, implies the full Erdős Problem 30 statement.

    For ε ≥ 1/2 the bilateral bound is already unconditional
    (`erdos30_partial_half`). For ε < 1/2 the lower side is supplied by the
    Singer construction: the prime-gap hypothesis is applied at the shifted
    perfect square `(Nat.sqrt N - t)²` with `t = ⌈N^ε⌉ + 2`, which forces the
    resulting prime `p` to satisfy `p + 2 ≤ Nat.sqrt N` — so Singer's Sidon
    set mod `p² + p + 1` fits inside `{1, …, N}` — while still keeping
    `p ≥ √N - O(N^ε)`. -/
theorem conditional_erdos30
    (h_prime_gap : ∀ ε : ℝ, 0 < ε →
      ∃ N₀ : ℕ, ∀ N ≥ N₀, ∃ p : ℕ, Nat.Prime p ∧
        |(p : ℝ) - Real.sqrt (N : ℝ)| ≤ Real.rpow (N : ℝ) ε)
    (h_upper : ∀ ε : ℝ, 0 < ε →
      ∃ C : ℝ, ∃ N0 : ℕ, 0 < C ∧
        ∀ {N h : ℕ}, N0 ≤ N → IsSidonMaximum N h →
          (h : ℝ) ≤ Real.sqrt (N : ℝ) + C * Real.rpow (N : ℝ) ε) :
    Erdos30Statement := by
  intro ε hε_pos
  by_cases hε_half : (1 : ℝ) / 2 ≤ ε
  · -- For ε ≥ 1/2 the bilateral bound is unconditional.
    exact erdos30_partial_half ε hε_half hε_pos
  push_neg at hε_half
  -- ε < 1/2: combine the upper hypothesis with the Singer/prime-gap lower bound.
  obtain ⟨C_up, N_up, hC_up_pos, hC_up⟩ := h_upper ε hε_pos
  obtain ⟨N_gap, hgap⟩ := h_prime_gap ε hε_pos
  have hγ_pos : (0 : ℝ) < 1 / 2 - ε := by linarith
  -- Threshold beyond which N^(1/2 - ε) ≥ 2, i.e. N^ε ≤ √N / 2.
  obtain ⟨n1, hn1⟩ := exists_nat_ge ((2 : ℝ) ^ ((1 : ℝ) / (1 / 2 - ε)))
  obtain ⟨B, hB_def⟩ : ∃ B : ℕ, B = N_gap + 1 := ⟨_, rfl⟩
  refine ⟨C_up + 5, n1 + (2 * B + 8) * (2 * B + 8) + N_up + 1, by linarith, ?_⟩
  intro N h hN hmax
  -- Write `Real.rpow` multiplicatively as `^` throughout.
  have hrpow_def : ∀ x y : ℝ, Real.rpow x y = x ^ y := fun _ _ => rfl
  simp only [hrpow_def] at hgap hC_up ⊢
  -- Basic bounds packaged into the chosen N₀.
  have hN_up : N_up ≤ N := by omega
  have hN1 : 1 ≤ N := by omega
  have hn1N : n1 ≤ N := by omega
  have hB_sq_le : (2 * B + 8) * (2 * B + 8) ≤ N := by omega
  have hN_real_one : (1 : ℝ) ≤ (N : ℝ) := by exact_mod_cast hN1
  have hN_real_pos : (0 : ℝ) < (N : ℝ) := by linarith
  have hrpow_pos : (0 : ℝ) < (N : ℝ) ^ ε := Real.rpow_pos_of_pos hN_real_pos ε
  have hrpow_ge_one : (1 : ℝ) ≤ (N : ℝ) ^ ε := by
    have h1 : (1 : ℝ) ^ ε ≤ (N : ℝ) ^ ε :=
      Real.rpow_le_rpow (by norm_num) hN_real_one hε_pos.le
    simpa using h1
  -- Step 1: 2 · N^ε ≤ √N once N ≥ n1.
  have h_two_le_rpow_gap : (2 : ℝ) ≤ (N : ℝ) ^ (1 / 2 - ε) := by
    have h2pos : (0 : ℝ) ≤ (2 : ℝ) ^ ((1 : ℝ) / (1 / 2 - ε)) := by positivity
    have hbase : (2 : ℝ) ^ ((1 : ℝ) / (1 / 2 - ε)) ≤ (N : ℝ) :=
      le_trans hn1 (by exact_mod_cast hn1N)
    have hmono : ((2 : ℝ) ^ ((1 : ℝ) / (1 / 2 - ε))) ^ (1 / 2 - ε) ≤
        (N : ℝ) ^ (1 / 2 - ε) :=
      Real.rpow_le_rpow h2pos hbase hγ_pos.le
    have heq : ((2 : ℝ) ^ ((1 : ℝ) / (1 / 2 - ε))) ^ (1 / 2 - ε) = 2 := by
      rw [← Real.rpow_mul (by norm_num : (0 : ℝ) ≤ 2),
        one_div_mul_cancel (ne_of_gt hγ_pos), Real.rpow_one]
    linarith [hmono, heq.le, heq.symm.le]
  have h_rpow_le_half_sqrt : 2 * (N : ℝ) ^ ε ≤ Real.sqrt (N : ℝ) := by
    have hsplit : Real.sqrt (N : ℝ) = (N : ℝ) ^ ε * (N : ℝ) ^ (1 / 2 - ε) := by
      rw [Real.sqrt_eq_rpow, ← Real.rpow_add hN_real_pos]
      congr 1
      ring
    calc
      2 * (N : ℝ) ^ ε = (N : ℝ) ^ ε * 2 := by ring
      _ ≤ (N : ℝ) ^ ε * (N : ℝ) ^ (1 / 2 - ε) :=
        mul_le_mul_of_nonneg_left h_two_le_rpow_gap hrpow_pos.le
      _ = Real.sqrt (N : ℝ) := hsplit.symm
  -- Step 2: 2B + 8 ≤ √N once N ≥ (2B + 8)².
  have h_B_le_sqrt : 2 * (B : ℝ) + 8 ≤ Real.sqrt (N : ℝ) := by
    have hcast : ((2 * B + 8 : ℕ) : ℝ) ≤ Real.sqrt (N : ℝ) := by
      calc
        ((2 * B + 8 : ℕ) : ℝ)
            = Real.sqrt (((2 * B + 8 : ℕ) : ℝ) * ((2 * B + 8 : ℕ) : ℝ)) :=
              (Real.sqrt_mul_self (Nat.cast_nonneg _)).symm
        _ ≤ Real.sqrt (N : ℝ) := Real.sqrt_le_sqrt (by exact_mod_cast hB_sq_le)
    have hpush : ((2 * B + 8 : ℕ) : ℝ) = 2 * (B : ℝ) + 8 := by push_cast; ring
    linarith [hcast, hpush.le, hpush.symm.le]
  -- The shift amount t and the shifted square target M = (Nat.sqrt N − t)².
  obtain ⟨t, ht_def⟩ : ∃ t : ℕ, t = ⌈(N : ℝ) ^ ε⌉₊ + 2 := ⟨_, rfl⟩
  have ht_le : (t : ℝ) ≤ (N : ℝ) ^ ε + 3 := by
    have hceil : (⌈(N : ℝ) ^ ε⌉₊ : ℝ) < (N : ℝ) ^ ε + 1 :=
      Nat.ceil_lt_add_one hrpow_pos.le
    rw [ht_def]
    push_cast
    linarith
  have ht_ge : (N : ℝ) ^ ε + 2 ≤ (t : ℝ) := by
    have hceil : (N : ℝ) ^ ε ≤ (⌈(N : ℝ) ^ ε⌉₊ : ℝ) := Nat.le_ceil _
    rw [ht_def]
    push_cast
    linarith
  have hs_le_sqrt : (Nat.sqrt N : ℝ) ≤ Real.sqrt (N : ℝ) := natSqrt_le_real_sqrt N
  have hsqrt_lt : Real.sqrt (N : ℝ) < (Nat.sqrt N : ℝ) + 1 :=
    real_sqrt_lt_natSqrt_add_one N
  have h_tB_le_s : t + B ≤ Nat.sqrt N := by
    have hreal : (t : ℝ) + (B : ℝ) ≤ (Nat.sqrt N : ℝ) := by linarith
    exact_mod_cast hreal
  obtain ⟨m, hm_def⟩ : ∃ m : ℕ, m = Nat.sqrt N - t := ⟨_, rfl⟩
  have hm_add : m + t = Nat.sqrt N := by omega
  have hm_B : B ≤ m := by omega
  obtain ⟨M, hM_def⟩ : ∃ M : ℕ, M = m * m := ⟨_, rfl⟩
  -- The target M is large enough for the prime-gap hypothesis…
  have hM_gap : N_gap ≤ M := by
    have hB1 : 1 ≤ B := by omega
    have h1 : B ≤ B * B := by nlinarith
    have h2 : B * B ≤ m * m := Nat.mul_le_mul hm_B hm_B
    omega
  -- …and sits below N.
  have hM_le_N : M ≤ N := by
    have h1 : m ≤ Nat.sqrt N := by omega
    have h2 : m * m ≤ Nat.sqrt N * Nat.sqrt N := Nat.mul_le_mul h1 h1
    have h3 : Nat.sqrt N * Nat.sqrt N ≤ N := Nat.le_sqrt.1 (le_refl (Nat.sqrt N))
    omega
  obtain ⟨p, hp_prime, hp_gap⟩ := hgap M hM_gap
  have hsqrtM : Real.sqrt (M : ℝ) = (m : ℝ) := by
    rw [hM_def]
    push_cast
    exact Real.sqrt_mul_self (Nat.cast_nonneg m)
  have hMrpow_le : (M : ℝ) ^ ε ≤ (N : ℝ) ^ ε :=
    Real.rpow_le_rpow (Nat.cast_nonneg M) (by exact_mod_cast hM_le_N) hε_pos.le
  have hp_abs : |(p : ℝ) - (m : ℝ)| ≤ (N : ℝ) ^ ε := by
    rw [← hsqrtM]
    exact le_trans hp_gap hMrpow_le
  have hp_bounds := abs_le.mp hp_abs
  have hm_real : (m : ℝ) + (t : ℝ) = (Nat.sqrt N : ℝ) := by exact_mod_cast hm_add
  -- The prime sits safely below Nat.sqrt N…
  have hp_add_two_le_s : p + 2 ≤ Nat.sqrt N := by
    have hreal : (p : ℝ) + 2 ≤ (Nat.sqrt N : ℝ) := by
      have h1 : (p : ℝ) - (m : ℝ) ≤ (N : ℝ) ^ ε := hp_bounds.2
      linarith [ht_ge]
    exact_mod_cast hreal
  -- …so the Singer construction at p fits inside {1, …, N}.
  have h_fits : p * p + p + 1 ≤ N := singer_modulus_fits hp_add_two_le_s
  obtain ⟨A, hA_sidon, hA_card⟩ := singerIntervalSidon p N hp_prime h_fits
  have hp_succ_le_h : p + 1 ≤ h := by
    have hcard_le := hmax.2 hA_sidon
    omega
  -- Lower side: √N − h ≤ 5 · N^ε.
  have h_lower_abs : Real.sqrt (N : ℝ) - (h : ℝ) ≤ 5 * (N : ℝ) ^ ε := by
    have hp_ge : (m : ℝ) - (N : ℝ) ^ ε ≤ (p : ℝ) := by
      have h1 : -((N : ℝ) ^ ε) ≤ (p : ℝ) - (m : ℝ) := hp_bounds.1
      linarith
    have hh_real : (p : ℝ) + 1 ≤ (h : ℝ) := by exact_mod_cast hp_succ_le_h
    linarith [hsqrt_lt, ht_le, hrpow_ge_one]
  -- Upper side from the hypothesis.
  have h_upper_abs : (h : ℝ) - Real.sqrt (N : ℝ) ≤ C_up * (N : ℝ) ^ ε := by
    have hup := hC_up hN_up hmax
    linarith
  have hCupR_pos : (0 : ℝ) < C_up * (N : ℝ) ^ ε := mul_pos hC_up_pos hrpow_pos
  rw [abs_le]
  constructor
  · linarith
  · linarith

/-- **Lower bound on sidonMaximum: (√N + 1) / 2 < sidonMaximum N for N ≥ 5.** -/
theorem sidonMaximum_gt_sqrt_div_two (N : ℕ) (hN : 5 ≤ N) :
    (Nat.sqrt N + 1) / 2 < sidonMaximum N := by
  by_cases h9 : 9 ≤ N
  · -- N ≥ 9: use Singer + Bertrand
    set m := Nat.sqrt N with hm_def
    set n := (m - 1) / 2 with hn_def
    have hm3 : 3 ≤ m := by
      rw [hm_def]
      exact Nat.le_sqrt.2 h9
    have hn_ne : n ≠ 0 := by
      intro hnz
      have hm_lt3 : m < 3 := by
        rw [hn_def] at hnz
        omega
      have : 3 ≤ m := hm3
      omega
    rcases Nat.exists_prime_lt_and_le_two_mul n hn_ne with ⟨p, hp, hnp, hp2n⟩
    have hpn : p ≤ m := by
      have h2n_plus1 : 2 * n + 1 ≤ m := by
        rw [hn_def]
        omega
      omega
    have hp_bound : p * p + p + 1 ≤ N := by
      have hN_sq : m * m ≤ N := by
        have := Nat.sqrt_le' N
        simpa [hm_def, pow_two] using this
      have hpm1 : p ≤ m - 1 := by
        have h2n_plus1 : 2 * n + 1 ≤ m := by
          rw [hn_def]; omega
        omega
      have hp_sq : p * p ≤ (m - 1) * (m - 1) := Nat.mul_le_mul hpm1 hpm1
      have : p * p + p + 1 ≤ m * m := by
        have hm_pos : 0 < m := by
          have : 3 ≤ m := hm3
          omega
        have h_le : p * p + p + 1 ≤ (m - 1) * (m - 1) + (m - 1) + 1 := by
          nlinarith
        have h_bound : (m - 1) * (m - 1) + (m - 1) + 1 ≤ m * m := by
          have hm_pos : 0 < m := by
            have : 3 ≤ m := hm3
            omega
          have h_eq : (m - 1) * (m - 1) + (m - 1) = (m - 1) * m := by
            calc
              (m - 1) * (m - 1) + (m - 1) = (m - 1) * ((m - 1) + 1) := by ring
              _ = (m - 1) * m := by
                have hm1 : 1 ≤ m := by omega
                calc
                  (m - 1) * ((m - 1) + 1) = (m - 1) * m := by
                    rw [Nat.sub_add_cancel hm1]
                  _ = (m - 1) * m := rfl
          calc
            (m - 1) * (m - 1) + (m - 1) + 1 = (m - 1) * m + 1 := by rw [h_eq]
            _ ≤ m * m := by
              have hlt : (m - 1) * m < m * m := Nat.mul_lt_mul_of_pos_right (by omega) hm_pos
              omega
        omega
      omega
    rcases singerIntervalSidon p N hp hp_bound with ⟨A, hA, hA_card⟩
    have hmax := sidonMaximum_isSidonMaximum N
    have hle : A.card ≤ sidonMaximum N := hmax.2 hA
    have hp_gt : (m + 1) / 2 < A.card := by
      have : (m + 1) / 2 ≤ n + 1 := by
        rw [hn_def]
        omega
      have hn_lt_p : n < p := hnp
      calc
        (m + 1) / 2 ≤ n + 1 := by
          rw [hn_def]
          omega
        _ ≤ p := by omega
        _ < p + 1 := by omega
        _ = A.card := by symm; exact hA_card
    have hm_card : (Nat.sqrt N + 1) / 2 < A.card := by
      simpa [hm_def] using hp_gt
    omega
  · -- 5 ≤ N < 9: direct verification
    have hN_range : N = 5 ∨ N = 6 ∨ N = 7 ∨ N = 8 := by omega
    rcases hN_range with rfl | rfl | rfl | rfl
    · -- N = 5: (√5 + 1) / 2 = 1 < sidonMaximum 5
      have h5 : (Nat.sqrt 5 + 1) / 2 = 1 := by native_decide
      have h_gt_1 : sidonMaximum 5 > 1 := by
        have h_exists : ∃ A : Finset ℤ, IsIntervalSidon (5 : ℤ) A ∧ A.card = 2 := by
          refine ⟨{1, 2}, ?_, by simp⟩
          refine ⟨?_, ?_⟩
          · intro a ha; simp at ha; rcases ha with rfl | rfl <;> norm_num
          · intro a b c d ha hb hc hd hsum
            simp at ha hb hc hd
            rcases ha with rfl | rfl <;> rcases hb with rfl | rfl <;>
              rcases hc with rfl | rfl <;> rcases hd with rfl | rfl <;> omega
        rcases h_exists with ⟨A, hA, hA_card⟩
        have hmax := sidonMaximum_isSidonMaximum 5
        have hle' : 2 ≤ sidonMaximum 5 := by
          have : A.card = 2 := hA_card
          have hle'' : A.card ≤ sidonMaximum 5 := hmax.2 hA
          omega
        omega
      rw [h5]
      exact h_gt_1
    · -- N = 6: (√6 + 1) / 2 = 1 < sidonMaximum 6
      have h6 : (Nat.sqrt 6 + 1) / 2 = 1 := by native_decide
      have h_gt_1 : sidonMaximum 6 > 1 := by
        have h_exists : ∃ A : Finset ℤ, IsIntervalSidon (6 : ℤ) A ∧ A.card = 2 := by
          refine ⟨{1, 2}, ?_, by simp⟩
          refine ⟨?_, ?_⟩
          · intro a ha; simp at ha; rcases ha with rfl | rfl <;> norm_num
          · intro a b c d ha hb hc hd hsum
            simp at ha hb hc hd
            rcases ha with rfl | rfl <;> rcases hb with rfl | rfl <;>
              rcases hc with rfl | rfl <;> rcases hd with rfl | rfl <;> omega
        rcases h_exists with ⟨A, hA, hA_card⟩
        have hmax := sidonMaximum_isSidonMaximum 6
        have hle' : 2 ≤ sidonMaximum 6 := by
          have : A.card = 2 := hA_card
          have hle'' : A.card ≤ sidonMaximum 6 := hmax.2 hA
          omega
        omega
      rw [h6]
      exact h_gt_1
    · -- N = 7: use Singer p=2
      have hp2 : Nat.Prime 2 := by decide
      have h_bound : 2 * 2 + 2 + 1 ≤ 7 := by norm_num
      rcases singerIntervalSidon 2 7 hp2 h_bound with ⟨A, hA, hA_card⟩
      have hA_card_eq : A.card = 2 + 1 := hA_card
      have hmax := sidonMaximum_isSidonMaximum 7
      have hle : A.card ≤ sidonMaximum 7 := hmax.2 hA
      have h_small_case : (Nat.sqrt 7 + 1) / 2 < A.card := by
        -- (Nat.sqrt 7 + 1) / 2 = (2 + 1) / 2 = 1, A.card = 3
        have : A.card = 3 := by omega
        rw [this]
        native_decide
      omega
    · -- N = 8: use Singer p=2
      have hp2 : Nat.Prime 2 := by decide
      have h_bound : 2 * 2 + 2 + 1 ≤ 8 := by norm_num
      rcases singerIntervalSidon 2 8 hp2 h_bound with ⟨A, hA, hA_card⟩
      have hA_card_eq : A.card = 2 + 1 := hA_card
      have hmax := sidonMaximum_isSidonMaximum 8
      have hle : A.card ≤ sidonMaximum 8 := hmax.2 hA
      have h_small_case : (Nat.sqrt 8 + 1) / 2 < A.card := by
        have : A.card = 3 := by omega
        rw [this]
        native_decide
      omega

/-- Combined unconditional two-sided bound on sidonMaximum:
    (√N + 1) / 2 < h(N) ≤ √(2N) + 1 for all N ≥ 5. -/
theorem sidonMaximum_bounds (N : ℕ) (hN : 5 ≤ N) :
    (Nat.sqrt N + 1) / 2 < sidonMaximum N ∧
      sidonMaximum N ≤ Nat.sqrt (2 * N) + 1 := by
  constructor
  · exact sidonMaximum_gt_sqrt_div_two N hN
  · exact sidonMaximum_le_sqrt_two N (by omega)

end Semantics.SidonSets