/-
BlockCoprimeDensity.lean — C(n) block-coprime density (finite Euler product)

Defines the finite Euler product that appears in the block-coprime density
problem.  For (r, M) ∈ ℕ² satisfying gcd(r, M+j) = 1 for j = 0..n, the
natural density is (analytically) known to be:

  C(n) = ζ(2) · ∏_{p prime} (1 − min(n+1, p) / p²)

WHAT IS FORMALIZED:
  • Finite truncation D_G(n) = ∏_{p ≤ G} (1 − min(n+1, p) / p²)  (ℚ)
  • Saturation partition: when p ≤ n+1 the factor simplifies to 1−1/p
  • Decomposition into product over saturated × active primes
  • C_G(n) = (∏_{p ≤ G} (1−1/p²)⁻¹) · D_G(n), with C_G(0) = 1 exact
  • Positivity of D_G(n) and the local factors
  • Complementarity of saturated/active primes
  • Boundary value at n=1 (Feller-Tornier product)
  • Eval witnesses for small G

WHAT IS NOT FORMALIZED (analytic number theory, beyond scope):
  • The infinite product limit G → ∞ (convergence)
  • The identification as a natural density of (r, M) pairs
  • The Mertens asymptotic D(n) ∼ e^{-γ} / log(n+1)
  • The connection to ζ(2) = π²/6

Structure:
  §1  Local factor and saturation partition
  §2  Finite Euler product D_G(n)
  §3  Product decomposition into saturated/active primes
  §4  Boundary values at n=0, n=1 and the ζ(2) cancellation
  §5  Positivity theorems
  §6  Notes on analytic extensions (unformalized)
  §7  Eval witnesses

References:
  - Wessen Getachew, "C(n) — Block-Coprime Density"
    https://wessengetachew.github.io/smith/
  - OEIS A013661 (ζ(2)), A065474 (∏(1-2/p²)), A065469 (C(1))
-/

import Mathlib.Data.Nat.Prime.Defs
import Mathlib.Data.Finset.Basic
import Mathlib.Data.Rat.Defs
import Mathlib.Data.Rat.Lemmas
import Mathlib.Tactic

open Finset
open Nat

namespace SilverSight.BlockCoprimeDensity

-- ═══════════════════════════════════════════════════════════════════════════
-- §1  Local factor and saturation partition
-- ═══════════════════════════════════════════════════════════════════════════

/-- Local factor for prime p at block length n:
    1 − min(n+1, p) / p²  (in ℚ).

    This is the contribution of prime p to the finite Euler product D_G(n).
    Defined for any ℕ p, but meaningful only for primes p ≥ 2. -/
def localFactor (n p : ℕ) : ℚ :=
  1 - (min (n+1) p : ℚ) / ((p : ℚ) ^ 2)

/-- A prime p is **saturated** at block length n when p ≤ n+1.
    Saturated primes contribute factor (1 − 1/p) instead of (1 − (n+1)/p²). -/
def isSaturated (n p : ℕ) : Prop :=
  p ≤ n + 1

instance (n p : ℕ) : Decidable (isSaturated n p) :=
  inferInstanceAs (Decidable (p ≤ n + 1))

/-- A prime p is **active** at block length n when p > n+1.
    Active primes still depend on the block length. -/
def isActive (n p : ℕ) : Prop :=
  n + 1 < p

instance (n p : ℕ) : Decidable (isActive n p) :=
  inferInstanceAs (Decidable (n + 1 < p))

/-- Every natural number is either saturated or active at block length n. -/
theorem isSaturated_or_isActive (n p : ℕ) : isSaturated n p ∨ isActive n p := by
  by_cases h : p ≤ n + 1
  · left; exact h
  · right; exact Nat.lt_of_not_ge h

/-- Saturated and active are complementary. -/
theorem isSaturated_iff_not_isActive (n p : ℕ) : isSaturated n p ↔ ¬isActive n p := by
  unfold isSaturated isActive
  exact ⟨Nat.not_lt.mpr, Nat.le_of_not_gt⟩

/-- For saturated p (p ≤ n+1), the local factor simplifies to 1 − 1/p.
    Since min(n+1, p) = p, we have 1 − p/p² = 1 − 1/p. -/
theorem localFactor_saturated (n p : ℕ) (h : isSaturated n p) : localFactor n p = 1 - (1 : ℚ) / (p : ℚ) := by
  unfold isSaturated at h
  unfold localFactor
  have hmin : (min (n+1) p : ℚ) = (p : ℚ) := by exact_mod_cast Nat.min_eq_right h
  rw [hmin]
  by_cases hzero : (p : ℚ) = 0
  · simp [hzero]
  · field_simp [hzero]

/-- For active p (p > n+1), the local factor is 1 − (n+1)/p².
    Since min(n+1, p) = n+1, this is immediate from the definition. -/
theorem localFactor_active (n p : ℕ) (h : isActive n p) : localFactor n p = 1 - ((n+1 : ℕ) : ℚ) / ((p : ℚ) ^ 2) := by
  unfold isActive at h
  unfold localFactor
  have hmin : (min (n+1) p : ℚ) = ((n+1 : ℕ) : ℚ) := by
    exact_mod_cast Nat.min_eq_left (by omega : n+1 ≤ p)
  rw [hmin]

-- ═══════════════════════════════════════════════════════════════════════════
-- §2  Finite Euler product
-- ═══════════════════════════════════════════════════════════════════════════

/-- The set of primes ≤ G as a Finset ℕ.
    When G = 0 or G = 1 the result is empty (no primes ≤ 1). -/
def primesUpTo (G : ℕ) : Finset ℕ :=
  (Finset.range (G+1)).filter Nat.Prime

/-- D_G(n) = ∏_{p ≤ G} (1 − min(n+1, p) / p²).
    The finite Euler product truncation.  This is rational for any finite G.
    The infinite limit D(n) = lim_{G→∞} D_G(n) is the raw block-coprime
    density (real, unformalized). -/
def D_finite (n G : ℕ) : ℚ :=
  Finset.prod (primesUpTo G) (fun p => localFactor n p)

/-- C_G(n) = (∏_{p ≤ G} (1−1/p²)⁻¹) · D_G(n).
    This is the finite G-truncation of C(n).  The product (∏ (1−1/p²)⁻¹)
    is the G-truncated Euler factor of ζ(2); the full identity ζ(2) =
    ∏_p (1−1/p²)⁻¹ is analytic and not proven here.

    For finite G this is rational; the infinite limit C(n) = ζ(2) · D(n)
    is real and unformalized. -/
def C_finite (n G : ℕ) : ℚ :=
  (Finset.prod (primesUpTo G) (fun p => (1 - (1 : ℚ) / ((p : ℚ) ^ 2))⁻¹)) * D_finite n G

-- ═══════════════════════════════════════════════════════════════════════════
-- §3  Product decomposition into saturated/active primes
-- ═══════════════════════════════════════════════════════════════════════════

/-- The saturated subset of primesUpTo G at block length n. -/
def saturatedPart (n G : ℕ) : Finset ℕ :=
  (primesUpTo G).filter (fun p => decide (isSaturated n p))

/-- The active subset of primesUpTo G at block length n. -/
def activePart (n G : ℕ) : Finset ℕ :=
  (primesUpTo G).filter (fun p => decide (isActive n p))

/-- The saturated and active parts partition primesUpTo G. -/
lemma saturatedPart_union_activePart (n G : ℕ) : saturatedPart n G ∪ activePart n G = primesUpTo G := by
  apply Finset.Subset.antisymm
  · intro p hp
    rcases Finset.mem_union.1 hp with (hp' | hp')
    · exact (Finset.mem_filter.1 hp').1
    · exact (Finset.mem_filter.1 hp').1
  · intro p hp
    rcases isSaturated_or_isActive n p with (h_sat | h_act)
    · apply Finset.mem_union_left
      refine Finset.mem_filter.mpr ⟨hp, ?_⟩
      exact decide_eq_true h_sat
    · apply Finset.mem_union_right
      refine Finset.mem_filter.mpr ⟨hp, ?_⟩
      exact decide_eq_true h_act

/-- The saturated and active parts are disjoint. -/
lemma saturatedPart_disjoint_activePart (n G : ℕ) : Disjoint (saturatedPart n G) (activePart n G) := by
  rw [Finset.disjoint_iff_inter_eq_empty]
  by_contra hne
  have h_nonempty : (saturatedPart n G ∩ activePart n G).Nonempty := by
    rw [Finset.nonempty_iff_ne_empty]
    exact hne
  rcases h_nonempty with ⟨p, hp⟩
  rcases Finset.mem_inter.1 hp with ⟨hp_sat, hp_act⟩
  have hp_sat_pair := Finset.mem_filter.mp (by simpa [saturatedPart] using hp_sat)
  have hp_act_pair := Finset.mem_filter.mp (by simpa [activePart] using hp_act)
  have hp_sat_pure : isSaturated n p := hp_sat_pair.2
  have hp_act_pure : isActive n p := hp_act_pair.2
  unfold isSaturated at hp_sat_pure
  unfold isActive at hp_act_pure
  omega

/-- D_finite(n, G) decomposes into a product over saturated primes times a
    product over active primes.  Applying `localFactor_saturated` and
    `localFactor_active` to each factor gives the explicit form:

    D_finite n G = (∏_{saturated} (1 − 1/p)) · (∏_{active} (1 − (n+1)/p²)) -/
theorem D_finite_eq_prod_saturated_mul_active (n G : ℕ) : D_finite n G =
    (Finset.prod (saturatedPart n G) (fun p => localFactor n p)) *
    (Finset.prod (activePart n G) (fun p => localFactor n p)) := by
  unfold D_finite
  calc
    Finset.prod (primesUpTo G) (fun p => localFactor n p)
        = Finset.prod (saturatedPart n G ∪ activePart n G) (fun p => localFactor n p) := by
      rw [saturatedPart_union_activePart]
    _ = (Finset.prod (saturatedPart n G) (fun p => localFactor n p)) *
        (Finset.prod (activePart n G) (fun p => localFactor n p)) := by
      rw [Finset.prod_union (saturatedPart_disjoint_activePart n G)]

-- ═══════════════════════════════════════════════════════════════════════════
-- §4  Boundary values and the ζ(2) cancellation
-- ═══════════════════════════════════════════════════════════════════════════

/-- For n=0, every prime is active (since the smallest prime is 2 > 1).
    The local factor at every prime is 1 − 1/p². -/
lemma localFactor_zero (p : ℕ) (hp : p ≠ 0) : localFactor 0 p = 1 - (1 : ℚ) / ((p : ℚ) ^ 2) := by
  unfold localFactor
  norm_num
  have hmin : (min 1 p : ℚ) = (1 : ℚ) :=
    by exact_mod_cast Nat.min_eq_left (Nat.one_le_of_lt (Nat.pos_of_ne_zero hp))
  rw [hmin]; simp

/-- D_G(0) = ∏_{p ≤ G} (1 − 1/p²) — the G-truncated Euler factor of ζ(2)^{-1}. -/
lemma D_finite_zero_eq (G : ℕ) : D_finite 0 G = Finset.prod (primesUpTo G) (fun p => 1 - (1 : ℚ) / ((p : ℚ) ^ 2)) := by
  unfold D_finite
  refine Finset.prod_congr rfl fun p hp => ?_
  have hp_prime : Nat.Prime p := (Finset.mem_filter.mp hp).2
  simp [localFactor_zero p (Nat.Prime.ne_zero hp_prime)]

/-- C_G(0) = 1 exactly for any G: at n=0, each factor (1 − 1/p²)
    cancels its own inverse from the ζ(2) expansion, regardless of G.

    This is an algebraic identity that holds for every finite truncation.
    The infinite limit C(0) = 1 also holds, but is a corollary of this
    finite identity, not an independent analytic statement. -/
@[simp] theorem C_finite_zero_eq_one (G : ℕ) : C_finite 0 G = 1 := by
  unfold C_finite
  have h : ∀ p ∈ primesUpTo G, (1 - (1 : ℚ) / ((p : ℚ) ^ 2))⁻¹ * (1 - (1 : ℚ) / ((p : ℚ) ^ 2)) = 1 := by
    intro p hp
    have hp_prime : Nat.Prime p := (Finset.mem_filter.mp hp).2
    have hx : (1 - (1 : ℚ) / ((p : ℚ) ^ 2)) ≠ 0 := by
      have hp_gt_one : (p : ℚ) > 1 := by exact_mod_cast (Nat.Prime.one_lt hp_prime)
      have hp_sq_gt_one : (p : ℚ) ^ 2 > 1 := by nlinarith
      have hp_sq_pos : 0 < (p : ℚ) ^ 2 := by nlinarith
      have hdiv : (1 : ℚ) / ((p : ℚ) ^ 2) < 1 := (div_lt_one hp_sq_pos).mpr hp_sq_gt_one
      nlinarith
    field_simp [hx]
  calc
    (Finset.prod (primesUpTo G) (fun p => (1 - (1 : ℚ) / ((p : ℚ) ^ 2))⁻¹)) * D_finite 0 G
        = (Finset.prod (primesUpTo G) (fun p => (1 - (1 : ℚ) / ((p : ℚ) ^ 2))⁻¹)) *
          (Finset.prod (primesUpTo G) (fun p => (1 - (1 : ℚ) / ((p : ℚ) ^ 2)))) := by rw [D_finite_zero_eq]
    _ = Finset.prod (primesUpTo G) (fun p => ((1 - (1 : ℚ) / ((p : ℚ) ^ 2))⁻¹ * (1 - (1 : ℚ) / ((p : ℚ) ^ 2)))) := by
      rw [← Finset.prod_mul_distrib]
    _ = Finset.prod (primesUpTo G) (fun _ => (1 : ℚ)) := by
      refine Finset.prod_congr rfl fun p hp => ?_
      exact h p hp
    _ = 1 := by simp

/-- For n=1, every prime p ≥ 2 has min(2, p) = 2.
    The local factor simplifies to 1 − 2/p². -/
lemma localFactor_one_eq (p : ℕ) (hp : 2 ≤ p) : localFactor 1 p = 1 - (2 : ℚ) / ((p : ℚ) ^ 2) := by
  unfold localFactor
  norm_num
  have hmin : (min 2 p : ℚ) = (2 : ℚ) := by exact_mod_cast Nat.min_eq_left hp
  rw [hmin]

/-- Dedicated version of `localFactor_one_eq` for primes.
    The hypothesis `Nat.Prime p` supplies `2 ≤ p` via `Nat.Prime.two_le`. -/
theorem localFactor_one_prime (p : ℕ) (hp : Nat.Prime p) : localFactor 1 p = 1 - (2 : ℚ) / ((p : ℚ) ^ 2) :=
  localFactor_one_eq p (Nat.Prime.two_le hp)

/-- D_G(1) = ∏_{p ≤ G} (1 − 2/p²) — the G-truncated Feller-Tornier Euler product.
    The infinite limit ∏_p (1−2/p²) ≈ 0.322634 is OEIS A065474. -/
lemma D_finite_one_eq (G : ℕ) : D_finite 1 G = Finset.prod (primesUpTo G) (fun p => 1 - (2 : ℚ) / ((p : ℚ) ^ 2)) := by
  unfold D_finite
  refine Finset.prod_congr rfl fun p hp => ?_
  have hp_prime : Nat.Prime p := (Finset.mem_filter.mp hp).2
  simp [localFactor_one_prime p hp_prime]

-- ═══════════════════════════════════════════════════════════════════════════
-- §5  Positivity
-- ═══════════════════════════════════════════════════════════════════════════

/-- For a prime p ≥ 2, the local factor is strictly positive. -/
lemma localFactor_pos (n p : ℕ) (hp : Nat.Prime p) : localFactor n p > 0 := by
  unfold localFactor
  have hp_gt_one : (p : ℚ) > 1 := by exact_mod_cast (Nat.Prime.one_lt hp)
  have hp_sq_gt_one : (p : ℚ) ^ 2 > 1 := by nlinarith
  have hp_sq_pos : 0 < (p : ℚ) ^ 2 := by nlinarith
  have h_num_le_denom : (min (n+1) p : ℚ) < (p : ℚ) ^ 2 := by
    calc
      (min (n+1) p : ℚ) ≤ (p : ℚ) := by exact_mod_cast Nat.min_le_right (n+1) p
      _ < (p : ℚ) ^ 2 := by nlinarith
  have h_div_lt_one : (min (n+1) p : ℚ) / ((p : ℚ) ^ 2) < 1 :=
    (div_lt_one hp_sq_pos).mpr h_num_le_denom
  nlinarith

/-- D_G(n) is strictly positive for any G. -/
theorem D_finite_pos (n G : ℕ) : D_finite n G > 0 := by
  unfold D_finite
  refine Finset.prod_pos fun p hp => ?_
  have hp_prime : Nat.Prime p := (Finset.mem_filter.mp hp).2
  exact localFactor_pos n p hp_prime

-- ═══════════════════════════════════════════════════════════════════════════
-- §6  Notes on analytic extensions (unformalized)
-- ═══════════════════════════════════════════════════════════════════════════

/-
Mertens' third theorem:  ∏_{p ≤ x} (1 − 1/p) ∼ e^{-γ} / log x.

If the limit D(n) = lim_{G→∞} D_G(n) exists, then Mertens implies:

  D(n) ∼ e^{-γ} / log(n+1)   as n → ∞

because the product over active primes (p > n+1) is asymptotically
∏_{p > n+1} (1 − (n+1)/p²) → 1, and the saturated product
∏_{p ≤ n+1} (1 − 1/p) ∼ e^{-γ} / log(n+1) by Mertens.

ANALYTIC BOUNDARY: Both the convergence lim_{G→∞} D_G(n) and the
Mertens asymptotic are analytic number theory results.  They are
documented here for reference only; this module does not prove them.
-/

/-
Conceptual note: SilverSight.SieveLemmas.depth_token_coprime_intersect
proves the existence/uniqueness of coprime-sieve CRT reconstruction
(the "quality" side).  The density C(n) — if the limit exists — would
be the asymptotic frequency of such coprime-block pairs (the "quantity"
side).  This module provides the Euler product expression; the formal
identification as a density is not proven here.
-/

-- ═══════════════════════════════════════════════════════════════════════════
-- §7  Eval witnesses
-- ═══════════════════════════════════════════════════════════════════════════

-- D_G(0) at G=31 (first 11 primes): ∏_{p ≤ 31} (1 − 1/p²) ≈ 0.61174
-- The infinite limit is 6/π² ≈ 0.607927.
#eval D_finite 0 31

-- D_G(1) at G=31: ∏_{p ≤ 31} (1 − 2/p²) ≈ 0.32666
-- The infinite limit (Feller-Tornier) ≈ 0.322634 (OEIS A065474).
#eval D_finite 1 31

-- C_G(0) = 1 exactly for any G.
#eval C_finite 0 31

-- Saturated primes at n=2: p ≤ 3 → {2, 3}
#eval List.filter (λ p : ℕ => decide (isSaturated 2 p)) [2, 3, 5, 7, 11]

-- Active primes at n=2: p > 3 → {5, 7, 11}
#eval List.filter (λ p : ℕ => decide (isActive 2 p)) [2, 3, 5, 7, 11]

end SilverSight.BlockCoprimeDensity
