import Mathlib

open Nat

/-
Theorem: gaps between squarefree numbers are unbounded.
For every k there is N such that none of N+1, ..., N+k is squarefree.

We prove by induction on k the strengthened statement that we may also produce
a modulus M ≥ 1 such that each witnessing prime square divides M. This lets the
inductive step merge a fresh prime square via the (two-modulus) CRT.
-/

theorem unbounded_squarefree_gaps (k : ℕ) :
    ∃ N : ℕ, ∀ i, 1 ≤ i → i ≤ k → ¬ Squarefree (N + i) := by
  suffices h : ∃ N M : ℕ, 1 ≤ M ∧
      ∀ i, 1 ≤ i → i ≤ k → ∃ p : ℕ, p.Prime ∧ p * p ∣ (N + i) ∧ p * p ∣ M by
    obtain ⟨N, M, _, hsq⟩ := h
    refine ⟨N, ?_⟩
    intro i hi1 hik
    obtain ⟨p, hp, hpd, _⟩ := hsq i hi1 hik
    rw [Nat.squarefree_iff_prime_squarefree]
    push_neg
    exact ⟨p, hp, hpd⟩
  induction k with
  | zero =>
      refine ⟨0, 1, le_refl 1, ?_⟩
      intro i hi1 hik; omega
  | succ n ih =>
      obtain ⟨N, M, hM1, hsq⟩ := ih
      obtain ⟨p, hpge, hp⟩ := Nat.exists_infinite_primes (M + n + 2)
      have hpM : ¬ p ∣ M := fun hdvdM => by
        have hple : p ≤ M := Nat.le_of_dvd (by omega) hdvdM
        omega
      have hcop : Nat.Coprime M (p * p) := by
        have h2 : Nat.Coprime M (p ^ 2) := hp.coprime_pow_of_not_dvd hpM
        simpa [pow_two] using h2
      have hpge1 : (n + 1) ≤ p := by omega
      have hpp_big : (n + 1) ≤ p * p := by nlinarith [hpge1, hp.two_le]
      set r : ℕ := p * p - (n + 1) with hr
      obtain ⟨N', hN'M, hN'p⟩ := Nat.chineseRemainder hcop (N % M) r
      refine ⟨N', M * (p * p), ?_, ?_⟩
      · have hp0 : 0 < p := hp.pos
        nlinarith [hM1, hp0]
      intro i hi1 hik
      rcases Nat.lt_or_ge i (n + 1) with hlt | hge
      · -- i ≤ n : reuse old witness q with q*q ∣ M
        have hile : i ≤ n := by omega
        obtain ⟨q, hq, hqd, hqM⟩ := hsq i hi1 hile
        refine ⟨q, hq, ?_, hqM.mul_right (p * p)⟩
        have e1 : N' ≡ N [MOD M] := hN'M.trans (Nat.mod_modEq N M)
        have e2 : N' + i ≡ N + i [MOD M] := e1.add_right i
        have e3 : N' + i ≡ N + i [MOD (q * q)] := e2.of_dvd hqM
        have hz : (N' + i) ≡ 0 [MOD (q * q)] :=
          e3.trans ((Nat.modEq_zero_iff_dvd).2 hqd)
        exact (Nat.modEq_zero_iff_dvd).1 hz
      · -- i = n+1 : new prime square
        have hie : i = n + 1 := by omega
        subst hie
        refine ⟨p, hp, ?_, (dvd_refl (p * p)).mul_left M⟩
        have hb : N' + (n + 1) ≡ r + (n + 1) [MOD (p * p)] := hN'p.add_right (n + 1)
        have hsum : r + (n + 1) = p * p := by omega
        rw [hsum] at hb
        have hz : (N' + (n + 1)) ≡ 0 [MOD (p * p)] :=
          hb.trans ((Nat.modEq_zero_iff_dvd).2 (dvd_refl _))
        exact (Nat.modEq_zero_iff_dvd).1 hz