/-
================================================================================
Waring–Chebyshev Structural Development

This file develops a division-free form of Waring’s formula for power sums
A^n + B^n by showing that:

  • The sequence Uₙ = A^n + B^n satisfies a Lucas-type recurrence
        U_{n+2} = (A+B) U_{n+1} − (AB) U_n

  • The Chebyshev polynomials Tₙ satisfy the same recurrence under the
    transformation x = (A + B)/2 in the unit case AB = 1.

  • Therefore power sums are determined purely by the symmetric quantities
        S = A + B
        P = A B

  • We explicitly unroll this recurrence to obtain Waring’s formula
    in a division-free binomial form.

  • Finally, we compress the coefficients into a single closed-form
    binomial expression:

        A^n + B^n =
        ∑_{j ≤ n/2} (-1)^j *
        ( n * C(n-j, j) / (n-j) ) *
        S^(n-2j) * P^j

Thus the power sum is expressed purely as a polynomial in S and P,
with no reference to A and B individually.

This completes the structural identification of Lucas sequences,
Chebyshev polynomials, and Waring’s formula in a single algebraic framework.
Code provided by ChatGPT-5.2, Claude 4.6 Sonnet, Claude Opus 4.6, and Gemini 3.1_pro_preview
with human assistance.
================================================================================
-/

import Mathlib.Analysis.SpecialFunctions.Trigonometric.Chebyshev.Basic
import Mathlib.Tactic

open Polynomial
open Chebyshev

/-- Core identity: `T_n` evaluated at `(x + x⁻¹)/2` equals `(x^n + x⁻¹^n)/2`. -/
lemma chebyshev_inv_identity {R : Type*} [Field R] [CharZero R]
    (x : R) (hx : x ≠ 0) (n : ℕ) :
    (T R n).eval ((x + x⁻¹) / 2) = (x ^ n + x⁻¹ ^ n) / 2 := by
  suffices key : ∀ k : ℕ,
      (T R k).eval ((x + x⁻¹) / 2) = (x ^ k + x⁻¹ ^ k) / 2 ∧
      (T R (k + 1)).eval ((x + x⁻¹) / 2) = (x ^ (k + 1) + x⁻¹ ^ (k + 1)) / 2 from
    (key n).1
  intro k
  induction k with
  | zero =>
      exact ⟨by simp [T_zero], by simp [T_one]⟩
  | succ k ih =>
      obtain ⟨ih_k, ih_succ⟩ := ih
      refine ⟨ih_succ, ?_⟩
      -- Shift the goal natively to `k + 2` to perfectly align with `T_add_two`
      change (T R (k + 2)).eval ((x + x⁻¹) / 2) = (x ^ (k + 2) + x⁻¹ ^ (k + 2)) / 2
      have hrec := congr_arg (Polynomial.eval ((x + x⁻¹) / 2)) (T_add_two R k)
      simp only [eval_sub, eval_mul, eval_ofNat, eval_X] at hrec
      rw [hrec, ih_k, ih_succ]
      -- Rearrange so x*x⁻¹ and x⁻¹*x appear explicitly, then cancel
      have h1 : x * x⁻¹ ^ (k + 1) = x⁻¹ ^ k := by
        calc x * x⁻¹ ^ (k + 1) = x⁻¹ ^ k * (x * x⁻¹) := by ring
        _ = x⁻¹ ^ k := by rw [mul_inv_cancel₀ hx, mul_one]
      have h2 : x⁻¹ * x ^ (k + 1) = x ^ k := by
        calc x⁻¹ * x ^ (k + 1) = x ^ k * (x⁻¹ * x) := by ring
        _ = x ^ k := by rw [inv_mul_cancel₀ hx, mul_one]
      linear_combination (1 / 2 : R) * h1 + (1 / 2 : R) * h2
-- end of chebyshev_inv_identity lemma


/-- Step 1: `A^n + B^n` satisfies the Lucas recurrence. -/
lemma sum_pow_lucas_rec {R : Type*} [CommRing R] (A B : R) (n : ℕ) :
    A ^ (n + 2) + B ^ (n + 2) =
      (A + B) * (A ^ (n + 1) + B ^ (n + 1)) - A * B * (A ^ n + B ^ n) := by
  ring
-- end of sum_pow_lucas_rec lemma


/-- Step 2: Unit case Chebyshev specialisation. -/
theorem lucas_chebyshev_unit {R : Type*} [Field R] [CharZero R]
    (A B : R) (hA : A ≠ 0) (hAB : A * B = 1) (n : ℕ) :
    A ^ n + B ^ n = 2 * (T R n).eval ((A + B) / 2) := by
  have hB : B = A⁻¹ := by
    calc B = A⁻¹ * (A * B) := by rw [← mul_assoc, inv_mul_cancel₀ hA, one_mul]
         _ = A⁻¹ := by rw [hAB, mul_one]
  subst hB
  rw [chebyshev_inv_identity A hA n]
  ring
-- end of lucas_chebyshev_unit theorem


/-- Corollary: trace formula for monic quadratic with unit constant term. -/
theorem chebyshev_trace_monic_unit {R : Type*} [Field R] [CharZero R]
    (A B : R) (hA : A ≠ 0) (hAB : A * B = 1) (n : ℕ) :
    A ^ n + B ^ n = 2 * (T R n).eval ((A + B) / 2) :=
  lucas_chebyshev_unit A B hA hAB n
-- end of chebyshev_trace_monic_unit theorem


/-- Power sums depend only on S = A+B and P = AB. -/
theorem sum_pow_determined_by_SP {R : Type*} [CommRing R] (A B : R) (n : ℕ) :
    ∃ f : R → R → R, A ^ n + B ^ n = f (A + B) (A * B) := by
  induction n using Nat.strongRecOn with
  | ind n ih =>
    rcases n with _ | _ | n
    · exact ⟨fun _ _ => 2, by ring⟩
    · exact ⟨fun s _ => s, by ring⟩
    · obtain ⟨fn, hfn⟩ := ih n (by omega)
      obtain ⟨gn, hgn⟩ := ih (n + 1) (by omega)
      refine ⟨fun s p => s * gn s p - p * fn s p, ?_⟩
      rw [sum_pow_lucas_rec, hfn, hgn]
-- end of sum_pow_determined_by_SP theorem

variable {R : Type*} [CommRing R]

/-- Unrolled coefficients for Waring's formula. -/
def waringCoeff' (n j : ℕ) : ℕ :=
  if j = 0 then 1
  else if n < 2 * j then 0
  else Nat.choose (n - j) j + Nat.choose (n - 1 - j) (j - 1)
-- end of waringCoeff' definition


/-- Core Waring sum formulation. -/
def waringSum' (n : ℕ) (S P : R) : R :=
  if n = 0 then 2 else
  ∑ j ∈ Finset.range (n / 2 + 1),
    (-1 : R)^j * (waringCoeff' n j : R) * S^(n - 2 * j) * P^j
-- end of waringSum' definition


/-- Coefficient recurrence alignment lemma. -/
lemma waringCoeff'_succ_add (n j : ℕ) (hn : 1 ≤ n) (hj : 1 ≤ j) :
    waringCoeff' (n + 2) j = waringCoeff' (n + 1) j + waringCoeff' n (j - 1) := by
  by_cases hj1 : j = 1
  · subst hj1
    have hw0 : waringCoeff' n (1 - 1) = 1 := rfl
    rw [hw0]
    unfold waringCoeff'
    have h_ne : 1 ≠ 0 := by decide
    have h1 : ¬(n + 2 < 2 * 1) := by omega
    have h2 : ¬(n + 1 < 2 * 1) := by omega
    simp only [h_ne, h1, h2, ↓reduceIte]
    -- Sub 1 - 1 to 0 explicitly to trigger zero_right
    have h_sub : 1 - 1 = 0 := rfl
    rw [h_sub]
    -- Apply Pascal's terminal conditions safely and let omega solve the algebra
    simp only [Nat.choose_one_right, Nat.choose_zero_right]
    omega
  · obtain ⟨k, rfl⟩ : ∃ k, j = k + 2 := ⟨j - 2, by omega⟩
    have hk1 : k + 2 - 1 = k + 1 := by omega
    rw [hk1]
    unfold waringCoeff'
    have hj0 : k + 2 ≠ 0 := by omega
    have hj1_ne : k + 1 ≠ 0 := by omega
    simp only [hj0, hj1_ne, ↓reduceIte]
    have eq_cond : n + 2 < 2 * (k + 2) ↔ n < 2 * (k + 1) := by omega
    by_cases c1 : n + 2 < 2 * (k + 2)
    · have c2 : n + 1 < 2 * (k + 2) := by omega
      have c3 : n < 2 * (k + 1) := eq_cond.mp c1
      simp only [c1, c2, c3, ↓reduceIte]
    · have c1_not : ¬(n + 2 < 2 * (k + 2)) := c1
      have c3_not : ¬(n < 2 * (k + 1)) := eq_cond.not.mp c1_not
      by_cases c2 : n + 1 < 2 * (k + 2)
      · -- Resolve `if` blocks *before* applying the specific variable
        --substitution to ensure syntactic match
        simp only [c1_not, c2, c3_not, ↓reduceIte, zero_add]
        have eq : n = 2 * k + 2 := by omega
        rw [eq]
        -- Normalise the left arguments
        have e1 : 2 * k + 2 + 2 - (k + 2) = k + 2 := by omega
        have e2 : 2 * k + 2 + 2 - 1 - (k + 2) = k + 1 := by omega
        have e3 : 2 * k + 2 - (k + 1) = k + 1 := by omega
        have e4 : 2 * k + 2 - 1 - (k + 1) = k := by omega
        -- Normalise the right arguments
        have e5 : k + 2 - 1 = k + 1 := by omega
        have e6 : k + 1 - 1 = k := by omega
        rw [e1, e2, e3, e4, e5, e6]
        -- Safely collapse all identical bounds to 1 universally
        simp only [Nat.choose_self]
      · have c2_not : ¬(n + 1 < 2 * (k + 2)) := c2
        simp only [c1_not, c2_not, c3_not, ↓reduceIte]
        -- Isolated exact indexing normalisation for left hand side arguments
        rw [show n + 2 - (k + 2) = n - k by omega]
        rw [show n + 2 - 1 - (k + 2) = n - k - 1 by omega]
        rw [show n + 1 - (k + 2) = n - k - 1 by omega]
        rw [show n + 1 - 1 - (k + 2) = n - k - 2 by omega]
        rw [show n - (k + 1) = n - k - 1 by omega]
        rw [show n - 1 - (k + 1) = n - k - 2 by omega]
        -- Isolated exact indexing normalisation for right hand side arguments
        have e_right1 : k + 2 - 1 = k + 1 := by omega
        have e_right2 : k + 1 - 1 = k := by omega
        rw [e_right1, e_right2]
        -- Prove structural identity substitutions directly
        --mapping to `.succ` to ensure perfect syntactic matches
        have H1 : Nat.choose (n - k) (k + 2) = Nat.choose (n - k - 1)
         (k + 1) + Nat.choose (n - k - 1) (k + 2) := by
          have h := Nat.choose_succ_succ (n - k - 1) (k + 1)
          have e1 : (n - k - 1).succ = n - k := by omega
          have e2 : (k + 1).succ = k + 2 := by omega
          rw [e1, e2] at h
          exact h
        have H2 : Nat.choose (n - k - 1) (k + 1) = Nat.choose (n - k - 2)
         k + Nat.choose (n - k - 2) (k + 1) := by
          have h := Nat.choose_succ_succ (n - k - 2) k
          have e1 : (n - k - 2).succ = n - k - 1 := by omega
          have e2 : k.succ = k + 1 := by omega
          rw [e1, e2] at h
          exact h
        rw [H1, H2]
        omega
-- end of waringCoeff'_succ_add lemma


/-- Range extension alignment lemma. -/
lemma waringSum'_eq_sum_range (n : ℕ) (S P : R) :
    waringSum' n S P = if n = 0 then 2 else ∑ j ∈ Finset.range (n + 2), (-1 : R)^j *
    (waringCoeff' n j : R) * S^(n - 2 * j) * P^j := by
  unfold waringSum'
  split_ifs with h
  · rfl
  · apply Finset.sum_subset
    · intro j hj; rw [Finset.mem_range] at hj ⊢; omega
    · intro j hj hj'
      rw [Finset.mem_range] at hj hj'
      have hj0 : j ≠ 0 := by omega
      have hn2j : n < 2 * j := by omega
      simp [waringCoeff', hj0, hn2j]
-- end of waringSum'_eq_sum_range lemma


/-- Full recurrence identity over polynomial coefficients. -/
lemma waringSum'_rec_eq (n : ℕ) (hn : 1 ≤ n) (S P : R) :
    waringSum' (n + 2) S P = S * waringSum' (n + 1) S P - P * waringSum' n S P := by
  have hn0 : (n = 0) ↔ False := iff_false_intro (by omega)
  have hn1 : (n + 1 = 0) ↔ False := iff_false_intro (by omega)
  have hn2 : (n + 2 = 0) ↔ False := iff_false_intro (by omega)
  rw [waringSum'_eq_sum_range (n + 2), waringSum'_eq_sum_range (n + 1), waringSum'_eq_sum_range n]
  simp only [hn0, hn1, hn2, ite_false]
  rw [Finset.mul_sum, Finset.mul_sum]
    -- Align backwards shifting sequences natively via
    --explicit hypothesis instantiation to avoid `+ 0` bleeding
  have sum_P_shift : ∑ j ∈ Finset.range (n + 3), (if j = 0 then (0 : R) else -
   ((-1 : R) ^ j * (waringCoeff' n (j - 1) : R) * S ^ (n + 2 - 2 * j) * P ^ j)) =
      ∑ j ∈ Finset.range (n + 2), P * ((-1 : R) ^ j * (waringCoeff' n j : R) *
       S ^ (n - 2 * j) * P ^ j) := by
    have h_split := Finset.sum_range_succ' (fun j => if j = 0 then (0 : R)
     else - ((-1 : R) ^ j * (waringCoeff' n (j - 1) : R) * S ^ (n + 2 - 2 * j) * P ^ j)) (n + 2)
    rw [h_split]
    have eq0 : (0 = 0) ↔ True := iff_true_intro rfl
    simp only [ite_true, add_zero]
    apply Finset.sum_congr rfl
    intro j _
    have hj0 : j + 1 ≠ 0 := by omega
    have hj0_iff : (j + 1 = 0) ↔ False := iff_false_intro hj0
    simp only [hj0_iff, ite_false, Nat.add_sub_cancel_right]
    have hS : n + 2 - 2 * (j + 1) = n - 2 * j := by omega
    rw [hS]
    ring
  have lhs_eq : ∑ j ∈ Finset.range (n + 4), (-1 : R) ^ j *
   (waringCoeff' (n + 2) j : R) * S ^ (n + 2 - 2 * j) * P ^ j =
      ∑ j ∈ Finset.range (n + 3), (-1 : R) ^ j * (waringCoeff' (n + 2) j : R)
       * S ^ (n + 2 - 2 * j) * P ^ j := by
    have h_split := Finset.sum_range_succ (fun j => (-1 : R) ^ j * (waringCoeff' (n + 2) j : R)
     * S ^ (n + 2 - 2 * j) * P ^ j) (n + 3)
    rw [h_split]
    have hz : waringCoeff' (n + 2) (n + 3) = 0 := by
      unfold waringCoeff'
      have h1 : (n + 3 = 0) ↔ False := iff_false_intro (by omega)
      have h2 : (n + 2 < 2 * (n + 3)) ↔ True := iff_true_intro (by omega)
      simp only [h1, h2, ite_false, ite_true]
    rw [hz]
    simp only [Nat.cast_zero, mul_zero, zero_mul, add_zero]
  rw [← sum_P_shift, lhs_eq, ← Finset.sum_sub_distrib]
  apply Finset.sum_congr rfl
  intro j _
  rcases j with _ | j
  · have eq0 : (0 = 0) ↔ True := iff_true_intro rfl
    simp only [ ite_true, sub_zero, pow_zero, mul_one]
    have c1 : waringCoeff' (n + 2) 0 = 1 := rfl
    have c2 : waringCoeff' (n + 1) 0 = 1 := rfl
    simp only [c1, c2, Nat.cast_one, one_mul]
    have e1 : n + 2 - 2 * 0 = n + 2 := rfl
    have e2 : n + 1 - 2 * 0 = n + 1 := rfl
    rw [e1, e2]
    ring
  · have hsucc : (j + 1 = 0) ↔ False := iff_false_intro (by omega)
    simp only [hsucc, ite_false, sub_neg_eq_add]
    rw [waringCoeff'_succ_add n (j + 1) hn (by omega)]
    push_cast
    by_cases h_bounds : n + 1 < 2 * (j + 1)
    · have z1 : waringCoeff' (n + 1) (j + 1) = 0 := by
        unfold waringCoeff'
        have h_bounds_true : (n + 1 < 2 * (j + 1)) ↔ True := iff_true_intro h_bounds
        simp only [hsucc, h_bounds_true, ite_false, ite_true]
      rw [z1]
      ring
    · have h_pow : n + 2 - 2 * (j + 1) = n + 1 - 2 * (j + 1) + 1 := by omega
      have hS : S ^ (n + 2 - 2 * (j + 1)) = S ^ (n + 1 - 2 * (j + 1))
       * S := by rw [h_pow, pow_add, pow_one]
      rw [hS]
      ring
-- end of waringSum'_rec_eq lemma


/-- Main division-free Waring formula. -/
theorem waring_formula_div_free (A B : R) (n : ℕ) :
    A^n + B^n = waringSum' n (A + B) (A * B) := by
  induction n using Nat.strongRecOn with
  | ind n ih =>
    rcases n with _ | _ | _ | n
    · simp [waringSum']; norm_num        -- ← add norm_num here
    · simp [waringSum', waringCoeff']
    · rw [waringSum']
      simp [waringCoeff', Finset.sum_range_succ]
      ring
    · have hrec : A ^ (n + 3) + B ^ (n + 3) =
          (A + B) * (A ^ (n + 2) + B ^ (n + 2)) - A * B * (A ^ (n + 1) + B ^ (n + 1)) := by ring
      rw [hrec, ih (n + 2) (by omega), ih (n + 1) (by omega)]
      exact (waringSum'_rec_eq (n + 1) (by omega) (A + B) (A * B)).symm
-- end of waring_formula_div_free theorem


/-- Coefficient compression lemma. -/
lemma waringCoeff_single_eq (n j : ℕ) (hn : n ≠ 0) (hj : j ≤ n / 2) :
    (if j = 0 then 1 else Nat.choose (n - j) j + Nat.choose (n - 1 - j) (j - 1)) =
    n * Nat.choose (n - j) j / (n - j) := by
  split_ifs with h
  · simp only [h, Nat.sub_zero, Nat.choose_zero_right, mul_one]
    exact (Nat.div_self (by omega)).symm
  · have h_nj_pos : 0 < n - j := by omega
    symm
    apply Nat.div_eq_of_eq_mul_right h_nj_pos
    have h1 : n - 1 - j + 1 = n - j := by omega
    have h2 : j - 1 + 1 = j := by omega
    have H := Nat.add_one_mul_choose_eq (n - 1 - j) (j - 1)
    rw [h1, h2] at H
    have hn_eq : n = n - j + j := (Nat.sub_add_cancel (by omega)).symm
    have expand_lhs : n * Nat.choose (n - j) j =
        (n - j) * Nat.choose (n - j) j + j * Nat.choose (n - j) j := by
      conv_lhs => rw [hn_eq]
      simp only [Nat.add_sub_cancel]
      ring
    have expand_rhs : (Nat.choose (n - j) j + Nat.choose (n - 1 - j) (j - 1)) * (n - j) =
        (n - j) * Nat.choose (n - j) j + (n - j) * Nat.choose (n - 1 - j) (j - 1) := by
      ring
    linarith [H]
-- end of waringCoeff_single_eq lemma


/-- Final Theorem: Clean single-sum Waring formula. -/
theorem waring_formula_div_free_sum_single (A B : R) {n : ℕ} (hn : n ≠ 0) :
    A^n + B^n = ∑ j ∈ Finset.range (n / 2 + 1),
      (-1 : R)^j * (((n * Nat.choose (n - j) j / (n - j)) : ℕ) : R)
       * (A + B)^(n - 2 * j) * (A * B)^j := by
  rw [waring_formula_div_free A B n, waringSum', if_neg hn]
  apply Finset.sum_congr rfl
  intro j hj
  rw [waringCoeff']
  have hj_bound : j ≤ n / 2 := by
    rw [Finset.mem_range] at hj; omega
  have h_not_lt : ¬(n < 2 * j) := by omega
  -- Collapse the leftover conditional structure completely safely
  have h_simplify : (if j = 0 then 1 else if n < 2 * j then 0 else
   Nat.choose (n - j) j + Nat.choose (n - 1 - j) (j - 1)) =
                    (if j = 0 then 1 else Nat.choose (n - j) j
                     + Nat.choose (n - 1 - j) (j - 1)) := by
    split_ifs with hj0
    · rfl
    · simp
  rw [h_simplify, waringCoeff_single_eq n j hn hj_bound]
-- end of waring_formula_div_free_sum_single theorem