import Mathlib

set_option linter.style.header false
set_option linter.style.longLine false

/-!
WP28 — machine-checked algebraic core of the proofs (Lean 4 + Mathlib).

What is formalized (deterministic statements only; the probabilistic O_p arguments are not):
  S1  Vandermonde moment constraints for the patterns (1,2), (1,2,4), (1,2,4,8), (1,2,3) and the constants beta_r.
  S3  Exact bias identities of J_2 and J_3 in the exponential model (Theorem 2(ii)).
  S4  Newton map and error identity in the exponential model (Theorem 2(iv)).
  S5  q-fold Newton error identity (Theorem 2(iv)) and the boundary constant algebra (Theorem 2(v)).
  S7  Newton basin: |u_j| = |u_0|^(2^j); convergence iff |u_0| < 1 (Theorem 2(iv)).
  S8  Exponent algebra of the frontier: with alpha(r,q) = 1 - 1/(2^(q+1) r) the scaled bias term is exactly 1 (Theorem 1 proof).
  S2/S6 Gamma inverse moments, E R_n, E R_n^2, Var(theta_hat_n), and the exact variance of J_2 from the Beta/Gamma
        decomposition (Theorem 2(i),(iii)); the moment values are derived from Real.Gamma_add_one.
-/


namespace WP28

/-! ### S1 Vandermonde weights -/

theorem vandermonde_12 : ((2:ℚ) + (-1) = 1) ∧ ((2:ℚ) * 1 + (-1) * 2 = 0) := by norm_num

theorem vandermonde_124 :
    ((8/3:ℚ) + (-2) + 1/3 = 1) ∧ ((8/3:ℚ) * 1 + (-2) * 2 + (1/3) * 4 = 0) ∧
    ((8/3:ℚ) * 1 + (-2) * 4 + (1/3) * 16 = 0) := by norm_num

theorem vandermonde_1248 :
    ((64/21:ℚ) - 8/3 + 2/3 - 1/21 = 1) ∧ ((64/21:ℚ) * 1 - (8/3) * 2 + (2/3) * 4 - (1/21) * 8 = 0) ∧
    ((64/21:ℚ) * 1 - (8/3) * 4 + (2/3) * 16 - (1/21) * 64 = 0) ∧
    ((64/21:ℚ) * 1 - (8/3) * 8 + (2/3) * 64 - (1/21) * 512 = 0) := by norm_num

theorem vandermonde_123 :
    ((3:ℚ) - 3 + 1 = 1) ∧ ((3:ℚ) * 1 - 3 * 2 + 1 * 3 = 0) ∧ ((3:ℚ) * 1 - 3 * 4 + 1 * 9 = 0) := by norm_num

/-- Leading bias coefficients beta_r = sum_j w_j s_j^r. -/
theorem beta_2 : (2:ℚ) * 1 + (-1) * 4 = -2 := by norm_num
theorem beta_3 : (8/3:ℚ) * 1 + (-2) * 8 + (1/3) * 64 = 8 := by norm_num
theorem beta_4 : (64/21:ℚ) * 1 - (8/3) * 16 + (2/3) * 256 - (1/21) * 4096 = -64 := by norm_num
theorem beta_3_pattern123 : (3:ℚ) * 1 - 3 * 8 + 1 * 27 = 6 := by norm_num

/-- Remainder amplification constants A_r = sum_j |w_j| s_j. -/
theorem amp_2 : |(2:ℚ)| * 1 + |(-1:ℚ)| * 2 = 4 := by norm_num
theorem amp_3 : |(8/3:ℚ)| * 1 + |(-2:ℚ)| * 2 + |(1/3:ℚ)| * 4 = 8 := by norm_num

/-! ### S3 Exact biases in the exponential model (Theorem 2(ii)) -/

/-- E J_2 - θ0 = θ0 Σ_j w_j s_j/(m - s_j) = -2θ0/((m-1)(m-2)). -/
theorem exact_bias_J2 (θ m : ℝ) (h1 : m - 1 ≠ 0) (h2 : m - 2 ≠ 0) :
    θ * (2 * 1 / (m - 1) + (-1) * 2 / (m - 2)) = -2 * θ / ((m - 1) * (m - 2)) := by
  field_simp <;> ring

/-- E J_3 - θ0 = 8θ0/((m-1)(m-2)(m-4)) for the pattern (1,2,4). -/
theorem exact_bias_J3 (θ m : ℝ) (h1 : m - 1 ≠ 0) (h2 : m - 2 ≠ 0) (h4 : m - 4 ≠ 0) :
    θ * ((8/3) * 1 / (m - 1) + (-2) * 2 / (m - 2) + (1/3) * 4 / (m - 4))
      = 8 * θ / ((m - 1) * (m - 2) * (m - 4)) := by
  field_simp <;> ring

/-! ### S4 Newton identities in the exponential model (Theorem 2(iv)) -/

/-- Newton map: x - g(x)/H(x) with g = 1/x - xbar, H = -1/x^2. -/
theorem newton_map (x xbar : ℝ) (hx : x ≠ 0) :
    x - (1 / x - xbar) / (-(1 / x ^ 2)) = 2 * x - x ^ 2 * xbar := by
  field_simp <;> ring

/-- Error identity: N(x) - θ̂ = -(x - θ̂)^2/θ̂ with θ̂ = 1/xbar. -/
theorem newton_error (x xbar : ℝ) (h : xbar ≠ 0) :
    (2 * x - x ^ 2 * xbar) - 1 / xbar = -((x - 1 / xbar) ^ 2) / (1 / xbar) := by
  field_simp <;> ring

/-! ### S5 q-fold Newton error and boundary constant -/

/-- Error recursion e_{j+1} = -e_j^2 / t with t = θ̂_N. -/
noncomputable def newtonErr (t : ℝ) : ℕ → ℝ → ℝ
  | 0, e => e
  | n + 1, e => -(newtonErr t n e) ^ 2 / t

/-- t^(2^q) · e_q = -t · e_0^(2^q) for q ≥ 1, i.e. e_q = -e_0^(2^q)/t^(2^q - 1). -/
theorem newtonErr_closed (t e : ℝ) (ht : t ≠ 0) :
    ∀ q : ℕ, t ^ (2 ^ (q + 1)) * newtonErr t (q + 1) e = -t * e ^ (2 ^ (q + 1)) := by
  intro q
  induction q with
  | zero =>
    simp [newtonErr]
    field_simp <;> ring
  | succ n ih =>
    have h2 : (2:ℕ) ^ (n + 2) = 2 ^ (n + 1) * 2 := by ring
    rw [newtonErr, h2, pow_mul, pow_mul]
    have : (t ^ 2 ^ (n + 1)) ^ 2 * (-(newtonErr t (n + 1) e) ^ 2 / t)
        = -((t ^ 2 ^ (n + 1) * newtonErr t (n + 1) e) ^ 2) / t := by
      field_simp <;> ring
    rw [this, ih]
    field_simp <;> ring

/-- Boundary constant (Theorem 2(v)), deterministic part: if the pilot error is
    e0 = β θ c^r N^(-1/2^(q+1)) then √N · (−e0^(2^q)/θ^(2^q−1)) = −θ (β c^r)^(2^q). -/
theorem boundary_constant (N θ β c : ℝ) (r q : ℕ) (hN : 0 < N) (hθ : θ ≠ 0) :
    let e0 := β * θ * c ^ r * N ^ (-(1 : ℝ) / 2 ^ (q + 1))
    N ^ ((1:ℝ) / 2) * (-(e0 ^ (2 ^ q)) / θ ^ (2 ^ q - 1)) = -θ * (β * c ^ r) ^ (2 ^ q) := by
  intro e0
  have hq : (2:ℕ) ^ q ≥ 1 := Nat.one_le_two_pow
  have hpow : (N ^ (-(1 : ℝ) / 2 ^ (q + 1))) ^ (2 ^ q) = N ^ (-(1:ℝ) / 2) := by
    rw [← Real.rpow_natCast, ← Real.rpow_mul hN.le]
    congr 1
    push_cast
    field_simp <;> ring
  have hθp : θ ^ (2 ^ q - 1) ≠ 0 := pow_ne_zero _ hθ
  have hroot : N ^ ((1:ℝ) / 2) * N ^ (-(1:ℝ) / 2) = 1 := by
    rw [← Real.rpow_add hN]; norm_num
  have hθs : θ ^ (2 ^ q) = θ * θ ^ (2 ^ q - 1) := by
    rw [← pow_succ']; congr 1; omega
  calc N ^ ((1:ℝ) / 2) * (-(e0 ^ (2 ^ q)) / θ ^ (2 ^ q - 1))
      = N ^ ((1:ℝ) / 2) * (-((β * θ * c ^ r) ^ (2 ^ q) * (N ^ (-(1 : ℝ) / 2 ^ (q + 1))) ^ (2 ^ q)) / θ ^ (2 ^ q - 1)) := by
          simp only [e0, mul_pow]
    _ = -(β * θ * c ^ r) ^ (2 ^ q) / θ ^ (2 ^ q - 1) * (N ^ ((1:ℝ) / 2) * N ^ (-(1:ℝ) / 2)) := by
          rw [hpow]; ring
    _ = -(β * θ * c ^ r) ^ (2 ^ q) / θ ^ (2 ^ q - 1) := by rw [hroot, mul_one]
    _ = -θ * (β * c ^ r) ^ (2 ^ q) := by
          rw [show (β * θ * c ^ r) ^ (2 ^ q) = θ ^ (2 ^ q) * (β * c ^ r) ^ (2 ^ q) by ring, hθs]
          field_simp <;> ring

/-! ### S7 Newton basin -/

/-- Normalized error recursion u_{j+1} = -u_j^2. -/
def basinIter : ℕ → ℝ → ℝ
  | 0, u => u
  | n + 1, u => -(basinIter n u) ^ 2

theorem basinIter_abs (u : ℝ) : ∀ j : ℕ, |basinIter j u| = |u| ^ (2 ^ j) := by
  intro j
  induction j with
  | zero => simp [basinIter]
  | succ n ih =>
    rw [basinIter, abs_neg, abs_pow, ih, ← pow_mul, pow_succ]

/-- Convergence of the Newton map in the exponential model when the start lies in (0, 2θ̂):
    |u_j| = |u|^(2^j) ≤ |u|^j → 0 for |u| < 1. -/
theorem basin_converges (u : ℝ) (hu : |u| < 1) :
    Filter.Tendsto (fun j => basinIter j u) Filter.atTop (nhds 0) := by
  rw [tendsto_zero_iff_abs_tendsto_zero]
  refine squeeze_zero (fun j => abs_nonneg _) (fun j => ?_)
    (tendsto_pow_atTop_nhds_zero_of_lt_one (abs_nonneg u) hu)
  show |basinIter j u| ≤ |u| ^ j
  rw [basinIter_abs]
  exact pow_le_pow_of_le_one (abs_nonneg u) hu.le (Nat.lt_two_pow_self).le

/-- Divergence when the start lies outside [0, 2θ̂]: |u_j| = |u|^(2^j) ≥ |u|^j → ∞ for |u| > 1. -/
theorem basin_diverges (u : ℝ) (hu : 1 < |u|) :
    Filter.Tendsto (fun j => |basinIter j u|) Filter.atTop Filter.atTop := by
  refine Filter.tendsto_atTop_mono (fun j => ?_) (tendsto_pow_atTop_atTop_of_one_lt hu)
  rw [basinIter_abs]
  exact pow_le_pow_right₀ hu.le (Nat.lt_two_pow_self).le

/-! ### S8 Exponent algebra of the frontier (proof of Theorem 1) -/

/-- With K = N^α, α = 1 − 1/(2^(q+1) r) and m = N/K = N^(1−α), the scaled bias term √N·m^(−r2^q) equals 1:
    the frontier exponent is exactly where the deterministic bias term stops vanishing. -/
theorem frontier_exponent (N : ℝ) (r q : ℕ) (hN : 0 < N) (hr : 0 < r) :
    let α : ℝ := 1 - 1 / (2 ^ (q + 1) * r)
    N ^ ((1:ℝ) / 2) * (N ^ ((1:ℝ) - α)) ^ (-((r : ℝ) * 2 ^ q)) = 1 := by
  intro α
  have hr' : (r : ℝ) ≠ 0 := by exact_mod_cast hr.ne'
  rw [← Real.rpow_mul hN.le, ← Real.rpow_add hN]
  have : (1:ℝ) / 2 + ((1:ℝ) - α) * (-((r : ℝ) * 2 ^ q)) = 0 := by
    simp only [α]
    field_simp <;> ring
  rw [this, Real.rpow_zero]

/-- Below the frontier, K = N^α' with α' < α, the scaled bias term is N^(−(α−α') r 2^q), a negative power of N. -/
theorem frontier_below (N : ℝ) (r q : ℕ) (α' : ℝ) (hN : 0 < N) (hr : 0 < r) :
    let α : ℝ := 1 - 1 / (2 ^ (q + 1) * r)
    N ^ ((1:ℝ) / 2) * (N ^ ((1:ℝ) - α')) ^ (-((r : ℝ) * 2 ^ q)) = N ^ (-((α - α') * (r * 2 ^ q))) := by
  intro α
  have hr' : (r : ℝ) ≠ 0 := by exact_mod_cast hr.ne'
  rw [← Real.rpow_mul hN.le, ← Real.rpow_add hN]
  congr 1
  simp only [α]
  field_simp <;> ring

/-! ### S2 / S6 Gamma moments and the exact variance of J_2 -/

/-- E[Y^{-1}] = n Γ(n−1)/Γ(n) = n/(n−1) for n > 1 (nY ~ Gamma(n,1)). -/
theorem gamma_inv_moment1 (n : ℝ) (hn : 1 < n) :
    n * Real.Gamma (n - 1) / Real.Gamma n = n / (n - 1) := by
  have h : Real.Gamma n = (n - 1) * Real.Gamma (n - 1) := by
    have := Real.Gamma_add_one (s := n - 1) (by linarith)
    rw [show n - 1 + 1 = n by ring] at this; exact this
  have hpos : Real.Gamma (n - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have h1 : n - 1 ≠ 0 := by linarith
  rw [h]; field_simp

/-- E[Y^{-2}] = n² Γ(n−2)/Γ(n) = n²/((n−1)(n−2)) for n > 2. -/
theorem gamma_inv_moment2 (n : ℝ) (hn : 2 < n) :
    n ^ 2 * Real.Gamma (n - 2) / Real.Gamma n = n ^ 2 / ((n - 1) * (n - 2)) := by
  have h1 : Real.Gamma n = (n - 1) * Real.Gamma (n - 1) := by
    have := Real.Gamma_add_one (s := n - 1) (by linarith)
    rw [show n - 1 + 1 = n by ring] at this; exact this
  have h2 : Real.Gamma (n - 1) = (n - 2) * Real.Gamma (n - 2) := by
    have := Real.Gamma_add_one (s := n - 2) (by linarith)
    rw [show n - 2 + 1 = n - 1 by ring] at this; exact this
  have hpos : Real.Gamma (n - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have e1 : n - 1 ≠ 0 := by linarith
  have e2 : n - 2 ≠ 0 := by linarith
  rw [h1, h2]; field_simp

/-- Theorem 2(i): with E Y = 1, E Y² = 1 + 1/n, E Y⁻¹ = n/(n−1), E Y⁻² = n²/((n−1)(n−2)):
    E R_n = θ0/(n−1), E R_n² = θ0²(3n+2)/(n(n−1)(n−2)), Var θ̂_n = θ0² n²/((n−1)²(n−2)). -/
theorem exponential_remainder_moments (θ n : ℝ) (hn : 2 < n) :
    θ * (n / (n - 1) - 2 + 1) = θ / (n - 1) ∧
    θ ^ 2 * (n ^ 2 / ((n - 1) * (n - 2)) - 4 * (n / (n - 1)) + 6 - 4 * 1 + (1 + 1 / n))
      = θ ^ 2 * (3 * n + 2) / (n * (n - 1) * (n - 2)) ∧
    θ ^ 2 * (n ^ 2 / ((n - 1) * (n - 2)) - (n / (n - 1)) ^ 2) = θ ^ 2 * n ^ 2 / ((n - 1) ^ 2 * (n - 2)) := by
  have h1 : n - 1 ≠ 0 := by linarith
  have h2 : n - 2 ≠ 0 := by linarith
  have h0 : n ≠ 0 := by linarith
  refine ⟨?_, ?_, ?_⟩ <;> (field_simp; ring)

/-- Beta(h,h) inverse moments as Gamma ratios reduce to rational functions of h (h > 2):
    E W⁻¹ = (2h−1)/(h−1), E W⁻² = (2h−1)(2h−2)/((h−1)(h−2)), E W⁻¹(1−W)⁻¹ = (2h−1)(2h−2)/(h−1)². -/
theorem beta_inverse_moments (h : ℝ) (hh : 2 < h) :
    Real.Gamma (h - 1) * Real.Gamma h * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 1))
        = (2 * h - 1) / (h - 1) ∧
    Real.Gamma (h - 2) * Real.Gamma h * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 2))
        = (2 * h - 1) * (2 * h - 2) / ((h - 1) * (h - 2)) ∧
    Real.Gamma (h - 1) * Real.Gamma (h - 1) * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 2))
        = (2 * h - 1) * (2 * h - 2) / (h - 1) ^ 2 := by
  have gh : Real.Gamma h = (h - 1) * Real.Gamma (h - 1) := by
    have := Real.Gamma_add_one (s := h - 1) (by linarith)
    rw [show h - 1 + 1 = h by ring] at this; exact this
  have gh1 : Real.Gamma (h - 1) = (h - 2) * Real.Gamma (h - 2) := by
    have := Real.Gamma_add_one (s := h - 2) (by linarith)
    rw [show h - 2 + 1 = h - 1 by ring] at this; exact this
  have g2h : Real.Gamma (2 * h) = (2 * h - 1) * Real.Gamma (2 * h - 1) := by
    have := Real.Gamma_add_one (s := 2 * h - 1) (by linarith)
    rw [show 2 * h - 1 + 1 = 2 * h by ring] at this; exact this
  have g2h1 : Real.Gamma (2 * h - 1) = (2 * h - 2) * Real.Gamma (2 * h - 2) := by
    have := Real.Gamma_add_one (s := 2 * h - 2) (by linarith)
    rw [show 2 * h - 2 + 1 = 2 * h - 1 by ring] at this; exact this
  have p1 : Real.Gamma (h - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have p2 : Real.Gamma (2 * h - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have p3 : Real.Gamma (h - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have p4 : Real.Gamma (2 * h - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'
  have p1' : Real.Gamma (-1 + h) ≠ 0 := by rw [show (-1:ℝ) + h = h - 1 by ring]; exact p3
  have p2' : Real.Gamma (-1 + h * 2) ≠ 0 := by rw [show (-1:ℝ) + h * 2 = 2 * h - 1 by ring]; exact p4
  have p3' : Real.Gamma (-2 + h) ≠ 0 := by rw [show (-2:ℝ) + h = h - 2 by ring]; exact p1
  have p4' : Real.Gamma (-2 + h * 2) ≠ 0 := by rw [show (-2:ℝ) + h * 2 = 2 * h - 2 by ring]; exact p2
  have f1 : (-1:ℝ) + h ≠ 0 := by linarith
  have f2 : (-2:ℝ) + h ≠ 0 := by linarith
  have f3 : (-1:ℝ) + h * 2 ≠ 0 := by linarith
  have f4 : (-2:ℝ) + h * 2 ≠ 0 := by linarith
  have q1 : h * 2 - 1 ≠ 0 := by linarith
  have q2 : h * 2 - 2 ≠ 0 := by linarith
  have q3 : Real.Gamma (h * 2 - 1) ≠ 0 := by rw [show h * 2 - 1 = 2 * h - 1 by ring]; exact p4
  have q4 : Real.Gamma (h * 2 - 2) ≠ 0 := by rw [show h * 2 - 2 = 2 * h - 2 by ring]; exact p2
  have q5 : Real.Gamma (2 * (h - 1)) ≠ 0 := by rw [show 2 * (h - 1) = 2 * h - 2 by ring]; exact p2
  have e1 : h - 1 ≠ 0 := by linarith
  have e2 : h - 2 ≠ 0 := by linarith
  have e3 : 2 * h - 1 ≠ 0 := by linarith
  have e4 : 2 * h - 2 ≠ 0 := by linarith
  refine ⟨?_, ?_, ?_⟩
  · rw [gh, g2h]; field_simp
  · rw [gh, gh1, g2h, g2h1]; field_simp
  · rw [gh, gh1, g2h, g2h1]; field_simp

/-- Theorem 2(iii): exact variance of J_2 (θ0 = 1, h = m/2). With E S⁻¹ = 1/(m−1), E S⁻² = 1/((m−1)(m−2)) and the
    Beta moments above, Var J_2 = E S⁻² E G² − (E S⁻¹ E G)² equals m²(m²−6m+17)/((m−4)(m−2)²(m−1)²),
    where G(W) = 2m − (h/2)(W⁻¹ + (1−W)⁻¹). -/
theorem var_J2_exact (m : ℝ) (hm : 4 < m) :
    (1 / ((m - 1) * (m - 2))) *
      (4 * m ^ 2 - 2 * (2 * m) * (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1)))
        + ((m / 2) ^ 2 / 4) * (2 * ((2 * (m / 2) - 1) * (2 * (m / 2) - 2) / ((m / 2 - 1) * (m / 2 - 2)))
                                + 2 * ((2 * (m / 2) - 1) * (2 * (m / 2) - 2) / (m / 2 - 1) ^ 2)))
      - ((1 / (m - 1)) * (2 * m - (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1))))) ^ 2
      = m ^ 2 * (m ^ 2 - 6 * m + 17) / ((m - 4) * (m - 2) ^ 2 * (m - 1) ^ 2) := by
  have e1 : m - 1 ≠ 0 := by linarith
  have e2 : m - 2 ≠ 0 := by linarith
  have e4 : m - 4 ≠ 0 := by linarith
  have e5 : m / 2 - 1 ≠ 0 := by intro h; linarith
  have e6 : m / 2 - 2 ≠ 0 := by intro h; linarith
  have f1 : (-1:ℝ) + m ≠ 0 := by linarith
  have f2 : (-2:ℝ) + m ≠ 0 := by linarith
  have f4 : (-4:ℝ) + m ≠ 0 := by linarith
  have g4 : m - 2 ^ 2 ≠ 0 := by norm_num; linarith
  field_simp
  ring

/-- Mean of J_2 from the same decomposition: E S⁻¹ · E G = 1 − 2/((m−1)(m−2)). -/
theorem mean_J2_exact (m : ℝ) (hm : 4 < m) :
    (1 / (m - 1)) * (2 * m - (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1))))
      = 1 - 2 / ((m - 1) * (m - 2)) := by
  have e1 : m - 1 ≠ 0 := by linarith
  have e2 : m - 2 ≠ 0 := by linarith
  have e5 : m / 2 - 1 ≠ 0 := by intro h; linarith
  have f1 : (-1:ℝ) + m ≠ 0 := by linarith
  have f2 : (-2:ℝ) + m ≠ 0 := by linarith
  field_simp
  ring

end WP28