/-- Positive-order cardinal polynomials are the classical Lagrange products. -/
theorem extremaCardinal_positive (N : {n : ℕ // 0 < n}) (j : Fin (N.val + 1)) :
    extremaCardinal N.val j = Lagrange.classicalBasis (Finset.range (N.val + 1))
      (fun i ↦ (classicalNode N i : ℂ)) j
      (Lagrange.BasisInput.of_injOn (Complex.ofReal_injective.comp_injOn (strictAntiOn_classicalNode N).injOn) (Finset.mem_range.mpr j.isLt)) := by
  simp [extremaCardinal, Nat.ne_of_gt N.property]

/-- Cardinal values on the entire finite grid. -/
theorem eval_extremaCardinal (N : ℕ) (i j : Fin (N + 1)) :
    (extremaCardinal N i).eval (extremaGridNode N j : ℂ) = if i = j then 1 else 0 := by
  by_cases hN : N = 0
  · subst N
    have hij : i = j := Fin.ext (by omega)
    simp [extremaCardinal_zero, hij]
  · rw [extremaCardinal_positive ⟨N,Nat.pos_of_ne_zero hN⟩,
      extremaGridNode_positive ⟨N,Nat.pos_of_ne_zero hN⟩]
    by_cases hij : i = j
    · subst j
      rw [ite_eq_left rfl]
      exact Lagrange.eval_classicalBasis_self (Finset.range (N+1))
        (fun k ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩ k : ℂ)) i.val _
    · rw [ite_eq_right hij]
      exact Lagrange.eval_classicalBasis_of_ne
        (v := fun k ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩ k : ℂ)) _ (fun h ↦ hij (Fin.ext h))
        (Finset.mem_range.mpr j.isLt)

/-- The extrema interpolant has degree at most `N`. -/
theorem degree_extremaInterpolant_lt (N : ℕ) (f : ℝ → ℂ) :
    (extremaInterpolant N f).degree < (N + 1 : ℕ) := by
  by_cases h : N = 0
  · subst N
    rw [extremaInterpolant_zero]
    simpa using (degree_C_lt (a := f 1))
  · rw [extremaInterpolant_positive ⟨N,Nat.pos_of_ne_zero h⟩]
    simpa only [Lagrange.classicalInterpolate,Finset.card_range] using
      Lagrange.degree_interpolate_lt (v := fun j ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero h⟩ j : ℂ))
        (fun j ↦ f (classicalNode ⟨N,Nat.pos_of_ne_zero h⟩ j)) (Complex.ofReal_injective.comp_injOn (strictAntiOn_classicalNode ⟨N,Nat.pos_of_ne_zero h⟩).injOn)

/-- Endpoint modes have discrete squared norm `N`, and all interior modes `N/2`. -/
theorem extremaSum_sq_eval_T {N k : ℕ} (hN : N ≠ 0) (hk : k ≤ N) :
    extremaSum ⟨N, by omega⟩ (fun m ↦ (T ℝ k).eval (node N m) ^ 2) =
      if k = 0 ∨ k = N then (N : ℝ) else N / 2 := by
  have h := extremaSum_product (N := N) (j := k) (k := k) (Nat.pos_of_ne_zero hN) le_rfl
  simp only [Nat.sub_self, Nat.cast_zero, T_zero, eval_one, extremaSum_one] at h
  by_cases hzero : k = 0
  · subst k
    simp
  · by_cases hlast : k = N
    · subst k
      rw [show N + N = 2 * N by omega, extremaSum_eval_T_two_mul hN] at h
      rw [ite_eq_left (Or.inr rfl)]
      simpa only [pow_two] using (show extremaSum ⟨N, by omega⟩
        (fun m ↦ (T ℝ N).eval (node N m) * (T ℝ N).eval (node N m)) = N by linarith)
    · rw [extremaSum_eval_T (by omega : k + k < 2 * N),
        ite_eq_right (by omega : k + k ≠ 0)] at h
      rw [ite_eq_right (by simp [hzero, hlast])]
      simpa only [pow_two] using (show extremaSum ⟨N, by omega⟩
        (fun m ↦ (T ℝ k).eval (node N m) * (T ℝ k).eval (node N m)) = (N : ℝ) / 2 by linarith)

/-- The DCT-I recovers each resolved Chebyshev mode exactly, including both endpoints. -/
theorem extremaCoefficient_eval_T {N j k : ℕ} (hN : N ≠ 0) (hj : j ≤ N) (hk : k ≤ N) :
    extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩ (fun x ↦ (T ℝ j).eval x)
      ⟨k, Nat.lt_succ_of_le hk⟩ = if j = k then 1 else 0 := by
  rw [extremaCoefficient_eq_sum hN]
  by_cases h : j = k
  · subst j
    rw [ite_eq_left rfl]
    simp only [← pow_two]
    rw [extremaSum_sq_eval_T hN hk]
    have hn : (N : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hN
    split_ifs <;> field_simp
  · rw [extremaSum_mul_eval_T_eq_zero hj hk h, mul_zero, ite_eq_right h]

/-- The DCT-I reconstructs classical distinct-node interpolation. -/
theorem dctInterpolant_eq_lagrange (N : {n : ℕ // 0 < n}) (f : ℝ → ℝ) :
    dctInterpolant N.val f =
      Lagrange.classicalInterpolate (Finset.range (N.val + 1)) (classicalNode N)
        ⟨(strictAntiOn_classicalNode N).injOn⟩ (fun j ↦ f (classicalNode N j)) := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  change dctInterpolant N f = Lagrange.interpolate (Finset.range (N+1)) (node N)
    (fun j ↦ f (node N j))
  let p := Lagrange.interpolate (Finset.range (N + 1)) (node N) (fun j ↦ f (node N j))
  have hv := (strictAntiOn_node N).injOn
  have hp : p.degree ≤ N := by
    have h := Lagrange.degree_interpolate_lt (fun j ↦ f (node N j)) hv
    simp only [Finset.card_range] at h
    by_cases hz : p = 0
    · simp [hz]
    · exact degree_le_of_natDegree_le
        (Nat.le_of_lt_succ ((natDegree_lt_iff_degree_lt hz).2 h))
  have he (j : ℕ) (hj : j ≤ N) : p.eval (node N j) = f (node N j) :=
    Lagrange.eval_interpolate_at_node _ hv (Finset.mem_range.mpr (by omega))
  change dctInterpolant N f = p
  rw [← dctInterpolant_eval hN hp]
  simp only [dctInterpolant, dite_eq_right hN]
  apply Finset.sum_congr rfl
  intro k hk
  congr 2
  simp only [extremaCoefficient_eq_sum hN]
  unfold extremaSum
  dsimp only
  rw [he N le_rfl, he 0 (Nat.zero_le N)]
  congr 2
  apply Finset.sum_congr rfl
  intro j hj
  rw [he j (by simpa using (Finset.mem_range.mp hj).le)]

/-- The inverse DCT-I takes exactly the supplied values at every extrema node. -/
theorem eval_dctInterpolant_node {N : ℕ} (hN : N ≠ 0) (f : ℝ → ℝ)
    {j : ℕ} (hj : j ≤ N) : (dctInterpolant N f).eval (node N j) = f (node N j) := by
  rw [dctInterpolant_eq_lagrange ⟨N,Nat.pos_of_ne_zero hN⟩]
  exact Lagrange.eval_classicalInterpolate
    (s := Finset.range (N+1)) (v := classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩)
    _ _ ⟨j,Finset.mem_range.mpr (Nat.lt_succ_of_le hj)⟩

/-- The DCT normalization is exactly the project's weighted-integral coefficient
normalization on the interpolating polynomial, with zero coefficients above `N`. -/
theorem coefficient_dctInterpolant (N : ℕ) (f : ℝ → ℝ) (k : ℕ) :
    (coefficient (n := k) ⟨(fun x ↦ (dctInterpolant N f).eval x), by first | assumption | apply integrable_classicalMeasureT; fun_prop⟩) =
      if hN : N = 0 then (if k = 0 then f 1 else 0)
      else if hk : k ≤ N then
        extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩ f ⟨k, Nat.lt_succ_of_le hk⟩ else 0 := by
  by_cases hN : N = 0
  · subst N
    simpa only [dctInterpolant_zero, eval_C, dite_true, coefficient,
      CoefficientInput.ofContinuousOn_val] using coefficient_const (f 1) k
  simp only [dite_eq_right hN, dctInterpolant, eval_finsetSum, eval_mul, eval_C]
  simp only [coefficient, Finset.sum_mul]
  rw [MeasureTheory.integral_finsetSum Finset.univ
      (fun i _ ↦ integrable_classicalMeasureT (by fun_prop)), Finset.mul_sum]
  simp only [mul_assoc, MeasureTheory.integral_const_mul]
  have hc (i : Fin (N + 1)) :
      (if k = 0 then 1 else 2) / Real.pi *
        (extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩ f i *
          ∫ x, (T ℝ i.val).eval x * (T ℝ k).eval x ∂classicalMeasureT) =
        extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩ f i *
          (if i.val = k then 1 else 0) := by
    rw [← coefficient_eval_T i.val k]
    simp only [coefficient, CoefficientInput.ofContinuousOn, CoefficientInput.ofL1]
    ring
  simp_rw [hc]
  by_cases hk : k ≤ N
  · rw [dite_eq_left hk, Finset.sum_eq_single (⟨k, Nat.lt_succ_of_le hk⟩ : Fin (N + 1))]
    · simp
    · intro b _ hb
      simp [show b.val ≠ k from fun h ↦ hb (Fin.ext h)]
    · simp
  · rw [dite_eq_right hk]
    apply Finset.sum_eq_zero
    intro b _
    simp [show b.val ≠ k by have := b.isLt; omega]

/-- Classical frequency folding uses a strictly positive modulus. -/
theorem aliasIndex_positive (N : {n : ℕ // 0 < n}) (k : ℕ) :
    aliasIndex N.val k = min (k % (2 * N.val)) (2 * N.val - k % (2 * N.val)) := by
  simp [aliasIndex, Nat.ne_of_gt N.property]

/-- A folded frequency has degree at most the interpolation degree. -/
theorem aliasIndex_le {N : ℕ} (hN : N ≠ 0) (k : ℕ) : aliasIndex N k ≤ N := by
  have hmod := Nat.mod_lt k (show 0 < 2 * N by omega)
  simp only [aliasIndex, ite_eq_right hN]
  omega

/-- Frequencies at or below the interpolation degree are unchanged. -/
theorem aliasIndex_eq_self {N k : ℕ} (hN : N ≠ 0) (hk : k ≤ N) : aliasIndex N k = k := by
  simp only [aliasIndex, ite_eq_right hN]
  rw [Nat.mod_eq_of_lt (by omega)]
  omega

/-- Node-level aliasing with the positive order in the node type (ATAP, Chapter 4). -/
theorem eval_T_node_aliasIndex (N : {n : ℕ // 0 < n}) (k j : ℕ) :
    (T ℝ k).eval (classicalNode N j) =
      (T ℝ (aliasIndex N.val k)).eval (classicalNode N j) := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  change (T ℝ k).eval (node N j) = (T ℝ (aliasIndex N k)).eval (node N j)
  rw [eval_T_node_mod hN]
  have hmod := Nat.mod_lt k (show 0 < 2 * N by omega)
  rw [aliasIndex_positive ⟨N,hNp⟩]
  dsimp only
  by_cases hk : k % (2 * N) ≤ N
  · rw [min_eq_left (by omega)]
  · rw [min_eq_right (by omega)]
    simp only [node, T_real_cos, Int.cast_natCast]
    have hn : (N : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hN
    have hang : ((2 * N - k % (2 * N) : ℕ) : ℝ) * (j * π / N) =
        (j : ℝ) * (2 * π) - (k % (2 * N) : ℕ) * (j * π / N) := by
      rw [Nat.cast_sub hmod.le]
      push_cast
      field_simp
    rw [hang, cos_nat_mul_two_pi_sub]

/-- Coefficient-level aliasing of a finite Chebyshev series (ATAP, Theorem 4.2). -/
theorem extremaCoefficient_sum_eval_T {N k : ℕ} (hN : N ≠ 0) (hk : k ≤ N)
    (s : Finset ℕ) (a : ℕ → ℝ) :
    extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩
      (fun x ↦ ∑ m ∈ s, a m * (T ℝ m).eval x) ⟨k, Nat.lt_succ_of_le hk⟩ =
      ∑ m ∈ s, if aliasIndex N m = k then a m else 0 := by
  rw [extremaCoefficient_sum]
  simp only [extremaCoefficient_mul, extremaCoefficient_eval_T_alias hN hk,
    mul_ite, mul_one, mul_zero]

/-- Coefficient aliasing for summable Chebyshev data, grouped by the exact
folding index. The nodal hypotheses use `HasSum` (unconditional convergence);
absolute summability of the source coefficients suffices (ATAP, Theorem 4.2). -/
theorem hasSum_alias_coefficients {N k : ℕ} (hN : N ≠ 0) (hk : k ≤ N)
    {a : ℕ → ℝ} {f : ℝ → ℝ}
    (hf : ∀ j : Fin (N + 1), HasSum (fun m ↦ a m * (T ℝ m).eval (extremaGridNode N j))
      (f (extremaGridNode N j))) :
    HasSum (fun m ↦ if aliasIndex N m = k then a m else 0) (extremaCoefficient ⟨N, Nat.pos_of_ne_zero hN⟩ f ⟨k, Nat.lt_succ_of_le hk⟩) := by
  have h := hasSum_extremaCoefficient (u := fun m x ↦ a m * (T ℝ m).eval x) (N := ⟨N, Nat.pos_of_ne_zero hN⟩) hf ⟨k, Nat.lt_succ_of_le hk⟩
  simpa only [extremaCoefficient_mul, extremaCoefficient_eval_T_alias hN hk,
    mul_ite, mul_one, mul_zero] using h

/-- Interpolating a single Chebyshev polynomial has norm at most one on the interval. -/
theorem norm_eval_extremaInterpolant_T_le_one {N : ℕ} (hN : N ≠ 0) (k : ℕ)
    {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    ‖(extremaInterpolant N (fun t : ℝ ↦ (((T ℝ k).eval t : ℝ) : ℂ))).eval (x : ℂ)‖ ≤ 1 := by
  rw [extremaInterpolant_T hN, ← complex_ofReal_eval_T, Complex.norm_real, Real.norm_eq_abs]
  exact abs_eval_T_real_le_one _ (abs_le.mpr hx)

/-- Evaluation as a finite sum of classical cardinal polynomials. -/
theorem eval_extremaInterpolant (N : ℕ) (f : ℝ → ℂ) (z : ℂ) :
    (extremaInterpolant N f).eval z = ∑ j : Fin (N + 1),
      f (extremaGridNode N j) * (extremaCardinal N j).eval z := by
  by_cases hN : N = 0
  · subst N
    simp [extremaInterpolant_zero, extremaCardinal_zero, extremaGridNode_zero]
  · rw [extremaInterpolant_positive ⟨N,Nat.pos_of_ne_zero hN⟩]
    simp only [extremaGridNode_positive ⟨N,Nat.pos_of_ne_zero hN⟩,
      extremaCardinal_positive ⟨N,Nat.pos_of_ne_zero hN⟩,
      Lagrange.classicalInterpolate, Lagrange.classicalBasis,
      Lagrange.interpolate_apply, Polynomial.eval_finsetSum,
      Polynomial.eval_mul, Polynomial.eval_C]
    exact (Fin.sum_univ_eq_sum_range (fun j : ℕ ↦
      f (classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩ j) *
        (Lagrange.basis (Finset.range (N+1))
          (fun k ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩ k : ℂ)) j).eval z) (N+1)).symm

/-- Interpolation commutes with an unconditionally convergent series at each node. -/
theorem hasSum_eval_extremaInterpolant {N : ℕ} {a : ℕ → ℂ} {f : ℝ → ℂ}
    (hs : ∀ j : Fin (N + 1), HasSum (fun k : ℕ ↦ a k * (((T ℝ k).eval (extremaGridNode N j) : ℝ) : ℂ))
      (f (extremaGridNode N j))) (z : ℂ) :
    HasSum (fun k : ℕ ↦ a k *
      (extremaInterpolant N (fun t : ℝ ↦ (((T ℝ k).eval t : ℝ) : ℂ))).eval z)
      ((extremaInterpolant N f).eval z) := by
  simp_rw [eval_extremaInterpolant, Finset.mul_sum, ← mul_assoc]
  exact hasSum_sum fun j hj ↦ (hs j).mul_right _

/-- Interpolation commutes with ordinary partial-sum limits at its finitely many
nodes; absolute convergence is not required (ATAP, Chapters 4 and 8). -/
theorem tendsto_sum_eval_extremaInterpolant {N : ℕ} {a : ℕ → ℂ} {f : ℝ → ℂ}
    (hs : ∀ j : Fin (N + 1), Filter.Tendsto
      (fun m : ℕ ↦ ∑ k ∈ Finset.range m,
        a k * (((T ℝ k).eval (extremaGridNode N j) : ℝ) : ℂ))
      Filter.atTop (nhds (f (extremaGridNode N j)))) (z : ℂ) :
    Filter.Tendsto
      (fun m : ℕ ↦ ∑ k ∈ Finset.range m, a k *
        (extremaInterpolant N (fun t : ℝ ↦ (((T ℝ k).eval t : ℝ) : ℂ))).eval z)
      Filter.atTop (nhds ((extremaInterpolant N f).eval z)) := by
  simp_rw [eval_extremaInterpolant, Finset.mul_sum, ← mul_assoc]
  simp_rw [Finset.sum_comm (s := Finset.range _), ← Finset.sum_mul]
  exact tendsto_finsetSum _ fun j _ ↦ (hs j).mul_const _

/-- Aliasing bounds interpolation error by twice the absolute coefficient tail
(ATAP, Chapter 4). -/
theorem norm_sub_extremaInterpolant_le_tsum_tail {N : ℕ} (hN : N ≠ 0)
    {a : ℕ → ℂ} {f : ℝ → ℂ} (ha : Summable (fun k ↦ ‖a k‖))
    (hf : ∀ x ∈ Icc (-1) 1,
      HasSum (fun k : ℕ ↦ a k * (((T ℝ k).eval x : ℝ) : ℂ)) (f x))
    {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    ‖f x - (extremaInterpolant N f).eval (x : ℂ)‖ ≤
      2 * ∑' k, ‖a (k + (N + 1))‖ := by
  let d (k : ℕ) : ℂ := a k * (((T ℝ k).eval x : ℝ) : ℂ) - a k *
    (extremaInterpolant N (fun t : ℝ ↦ (((T ℝ k).eval t : ℝ) : ℂ))).eval (x : ℂ)
  have hd : HasSum d (f x - (extremaInterpolant N f).eval (x : ℂ)) :=
    (hf x hx).sub (hasSum_eval_extremaInterpolant (fun j ↦ hf _ (extremaGridNode_mem_Icc N j)) (x : ℂ))
  have hz (k : ℕ) (hk : k ≤ N) : d k = 0 := by
    dsimp [d]
    rw [extremaInterpolant_T hN, aliasIndex_eq_self hN hk, ← complex_ofReal_eval_T,
      sub_self]
  have hb (k : ℕ) : ‖d k‖ ≤ 2 * ‖a k‖ := by
    dsimp [d]
    apply (norm_sub_le _ _).trans
    rw [norm_mul, norm_mul, Complex.norm_real, Real.norm_eq_abs]
    have ht := abs_eval_T_real_le_one (k : ℤ) (abs_le.mpr hx)
    have hi := norm_eval_extremaInterpolant_T_le_one hN k hx
    nlinarith [norm_nonneg (a k)]
  have ht := ((summable_nat_add_iff (N + 1)).2 ha).mul_left 2
  have hn : Summable (fun k ↦ ‖d (k + (N + 1))‖) :=
    ht.of_norm_bounded (fun k ↦ by simpa only [norm_norm] using hb (k + (N + 1)))
  have he := hd.summable.sum_add_tsum_nat_add (N + 1)
  rw [hd.tsum_eq, Finset.sum_eq_zero (fun k hk ↦ hz k (by simpa using hk)), zero_add] at he
  rw [← he, ← tsum_mul_left]
  exact (norm_tsum_le_tsum_norm hn).trans (hn.tsum_le_tsum (fun k ↦ hb _) ht)

/-- Classical product weights at positive-order extrema (ATAP, Chapter 5). -/
theorem nodalWeight_extrema (N : {n : ℕ // 0 < n}) (j : Fin (N.val + 1)) :
    Lagrange.classicalNodalWeight (Finset.range (N.val + 1)) (classicalNode N) j
      (Lagrange.BasisInput.of_injOn (strictAntiOn_classicalNode N).injOn
        (Finset.mem_range.mpr j.isLt)) =
      (extremaNodal N).leadingCoeff / (N.val : ℝ) ^ 2 *
        extremaBarycentricWeight N j := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  rcases j with ⟨j,hjp⟩
  change j < N+1 at hjp
  have hj := Nat.le_of_lt_succ hjp
  change Lagrange.nodalWeight (Finset.range (N + 1)) (node N) j =
      (extremaNodal ⟨N, by omega⟩).leadingCoeff / (N : ℝ) ^ 2 * extremaBarycentricWeight ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩
  have hmem : j ∈ Finset.range (N + 1) := Finset.mem_range.mpr (by omega)
  let d := (derivative (Lagrange.nodal (Finset.range (N + 1)) (node N))).eval (node N j)
  have hd : d ≠ 0 := by
    dsimp only [d]
    rw [Lagrange.eval_nodal_derivative_eval_node_eq hmem]
    apply Lagrange.eval_nodal_not_at_node
    intro i hi hji
    exact (Finset.mem_erase.mp hi).1
      ((strictAntiOn_classicalNode ⟨N,Nat.pos_of_ne_zero hN⟩).injOn
        (Finset.mem_erase.mp hi).2 hmem hji.symm)
  have he : (extremaNodal ⟨N, by omega⟩).leadingCoeff * d =
      extremaEndpointFactor ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩ * (N : ℝ) ^ 2 * (-1) ^ j := by
    have h := congrArg (fun p : ℝ[X] ↦ (derivative p).eval (node N j))
      (extremaNodal_eq_leadingCoeff_mul_nodal hN)
    simpa only [eval_derivative_extremaNodal_node hN hj, derivative_C_mul,
      eval_mul, eval_C] using h.symm
  rw [Lagrange.nodalWeight_eq_eval_derivative_nodal hmem]
  change d⁻¹ = _
  apply (mul_right_cancel₀ hd)
  rw [inv_mul_cancel₀ hd]
  have hNr : (N : ℝ) ≠ 0 := by exact_mod_cast hN
  have hE := (extremaEndpointFactor_pos ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩).ne'
  unfold extremaBarycentricWeight
  field_simp [hNr, hE]
  have hs : (-1 : ℝ) ^ j * (-1) ^ j = 1 := by rw [← mul_pow]; simp
  linear_combination -(-1 : ℝ) ^ j * he -
    extremaEndpointFactor ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩ * (N : ℝ) ^ 2 * hs

/-- The common scaling factor in the barycentric weights is nonzero (ATAP, Chapter 5). -/
theorem extremaBarycentricScale_ne_zero {N : ℕ} (hN : N ≠ 0) :
    (extremaNodal ⟨N, by omega⟩).leadingCoeff / (N : ℝ) ^ 2 ≠ 0 :=
  div_ne_zero (leadingCoeff_ne_zero.mpr (extremaNodal_ne_zero hN))
    (pow_ne_zero _ (by exact_mod_cast hN))

/-- Complex classical product weights (ATAP, Chapter 5). -/
theorem nodalWeight_complex_extrema (N : {n : ℕ // 0 < n}) (j : Fin (N.val + 1)) :
    Lagrange.classicalNodalWeight (Finset.range (N.val + 1))
      (fun i ↦ (classicalNode N i : ℂ)) j
      (complex_basisInput (Lagrange.BasisInput.of_injOn
        (strictAntiOn_classicalNode N).injOn (Finset.mem_range.mpr j.isLt))) =
      (((extremaNodal N).leadingCoeff / (N.val : ℝ) ^ 2 : ℝ) : ℂ) *
        (extremaBarycentricWeight N j : ℂ) := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  rcases j with ⟨j,hjp⟩
  change j < N+1 at hjp
  have hj := Nat.le_of_lt_succ hjp
  change Lagrange.nodalWeight (Finset.range (N + 1)) (fun i ↦ (node N i : ℂ)) j =
      (((extremaNodal ⟨N, by omega⟩).leadingCoeff / (N : ℝ) ^ 2 : ℝ) : ℂ) *
        (extremaBarycentricWeight ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩ : ℂ)
  have h := congrArg Complex.ofReal (nodalWeight_extrema ⟨N,hNp⟩ ⟨j,hjp⟩)
  simpa [Lagrange.classicalNodalWeight, classicalNode, Lagrange.nodalWeight] using h

/-- Nonvanishing on the classical off-grid domain (ATAP, Chapter 5). -/
theorem extremaBarycentric_denominator_ne_zero (N : {n : ℕ // 0 < n})
    (z : {z : ℂ // ∀ j ≤ N.val, z ≠ (classicalNode N j : ℂ)}) :
    (∑ j : Fin (N.val + 1), (extremaBarycentricWeight N j : ℂ) /
      (z.val - (classicalNode N j : ℂ))) ≠ 0 := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  rcases z with ⟨z,hz⟩
  change ∀ j ≤ N, z ≠ (node N j : ℂ) at hz
  have hv : Set.InjOn (fun j ↦ (node N j : ℂ)) (Finset.range (N+1)) :=
    Complex.ofReal_injective.comp_injOn (strictAntiOn_classicalNode ⟨N,hNp⟩).injOn
  change (∑ j : Fin (N + 1), (extremaBarycentricWeight ⟨N, by omega⟩ j : ℂ) /
      (z - (node N j : ℂ))) ≠ 0
  have h := Lagrange.sum_nodalWeight_mul_inv_sub_ne_zero hv
    (fun j hj ↦ hz j (by simpa using hj)) (Finset.nonempty_range_iff.mpr (by omega))
  have he : (∑ j ∈ Finset.range (N + 1),
      Lagrange.nodalWeight (Finset.range (N + 1)) (fun i ↦ (node N i : ℂ)) j *
        (z - (node N j : ℂ))⁻¹) =
      (((extremaNodal ⟨N, by omega⟩).leadingCoeff / (N : ℝ) ^ 2 : ℝ) : ℂ) *
        ∑ j : Fin (N + 1), (extremaBarycentricWeight ⟨N, by omega⟩ j : ℂ) /
          (z - (node N j : ℂ)) := by
    rw [Finset.mul_sum, ← Fin.sum_univ_eq_sum_range]
    apply Finset.sum_congr rfl
    intro j _
    have hw := nodalWeight_complex_extrema ⟨N,hNp⟩ j
    change Lagrange.nodalWeight (Finset.range (N+1)) (fun i ↦ (node N i : ℂ)) j = _ at hw
    rw [hw]
    simp only [div_eq_mul_inv, mul_assoc]
  rw [he] at h
  exact (mul_ne_zero_iff.mp h).2

/-- Barycentric evaluation on its classical domain (ATAP, Chapter 5). -/
theorem eval_extremaInterpolant_barycentric (N : {n : ℕ // 0 < n}) (f : ℝ → ℂ)
    (z : {z : ℂ // ∀ j ≤ N.val, z ≠ (classicalNode N j : ℂ)}) :
    (extremaInterpolant N.val f).eval z.val =
      (∑ j : Fin (N.val + 1), (extremaBarycentricWeight N j : ℂ) /
        (z.val - (classicalNode N j : ℂ)) * f (classicalNode N j)) /
      (∑ j : Fin (N.val + 1), (extremaBarycentricWeight N j : ℂ) /
        (z.val - (classicalNode N j : ℂ))) := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  rcases z with ⟨z,hz⟩
  change ∀ j ≤ N, z ≠ (node N j : ℂ) at hz
  have hv : Set.InjOn (fun j ↦ (node N j : ℂ)) (Finset.range (N+1)) :=
    Complex.ofReal_injective.comp_injOn (strictAntiOn_classicalNode ⟨N,hNp⟩).injOn
  change (extremaInterpolant N f).eval z =
      (∑ j : Fin (N + 1), (extremaBarycentricWeight ⟨N, by omega⟩ j : ℂ) /
        (z - (node N j : ℂ)) * f (node N j)) /
      (∑ j : Fin (N + 1), (extremaBarycentricWeight ⟨N, by omega⟩ j : ℂ) /
        (z - (node N j : ℂ)))
  rw [extremaInterpolant_positive ⟨N,Nat.pos_of_ne_zero hN⟩]
  change (Lagrange.interpolate (Finset.range (N+1)) (fun j ↦ (node N j : ℂ))
    (fun j ↦ f (node N j))).eval z = _
  rw [Lagrange.eval_interpolate_not_at_node'
    _ hv (Finset.nonempty_range_iff.mpr (by omega))
      (fun j hj ↦ hz j (by simpa using hj))]
  have hsum (r : ℕ → ℂ) :
      (∑ j ∈ Finset.range (N + 1),
        Lagrange.nodalWeight (Finset.range (N + 1)) (fun i ↦ (node N i : ℂ)) j *
          (z - (node N j : ℂ))⁻¹ * r j) =
      (((extremaNodal ⟨N, by omega⟩).leadingCoeff / (N : ℝ) ^ 2 : ℝ) : ℂ) *
        ∑ j : Fin (N + 1), (extremaBarycentricWeight ⟨N, by omega⟩ j : ℂ) /
          (z - (node N j : ℂ)) * r j := by
    rw [Finset.mul_sum, ← Fin.sum_univ_eq_sum_range]
    apply Finset.sum_congr rfl
    intro j _
    have hw := nodalWeight_complex_extrema ⟨N,hNp⟩ j
    change Lagrange.nodalWeight (Finset.range (N+1)) (fun i ↦ (node N i : ℂ)) j = _ at hw
    rw [hw]
    simp only [div_eq_mul_inv, mul_assoc]
  rw [hsum]
  have hone := hsum (fun _ ↦ 1)
  simp only [mul_one] at hone
  rw [hone, mul_div_mul_left]
  exact_mod_cast extremaBarycentricScale_ne_zero hN

/-- Classical real cardinal evaluation off its own node (ATAP, Chapter 5). -/
theorem eval_extremaBasis (N : {n : ℕ // 0 < n}) (j : Fin (N.val + 1))
    (x : {x : ℝ // x ≠ classicalNode N j}) :
    (Lagrange.classicalBasis (Finset.range (N.val + 1)) (classicalNode N) j
      (Lagrange.BasisInput.of_injOn (strictAntiOn_classicalNode N).injOn
        (Finset.mem_range.mpr j.isLt))).eval x.val =
      (extremaNodal N).eval x.val * (-1) ^ (j : ℕ) /
        ((N.val : ℝ) ^ 2 * extremaEndpointFactor N j * (x.val - classicalNode N j)) := by
  rcases N with ⟨N,hNp⟩
  have hN := Nat.ne_of_gt hNp
  rcases j with ⟨j,hjp⟩
  change j < N+1 at hjp
  have hj := Nat.le_of_lt_succ hjp
  rcases x with ⟨x,hx⟩
  change x ≠ node N j at hx
  change (Lagrange.basis (Finset.range (N + 1)) (node N) j).eval x =
      (extremaNodal ⟨N, by omega⟩).eval x * (-1) ^ j /
        ((N : ℝ) ^ 2 * extremaEndpointFactor ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩ * (x - node N j))
  have hlc := leadingCoeff_ne_zero.mpr (extremaNodal_ne_zero hN)
  have hn : (N : ℝ) ≠ 0 := by exact_mod_cast hN
  have he := congrArg (fun p : ℝ[X] ↦ p.eval x)
    (extremaNodal_eq_leadingCoeff_mul_nodal hN)
  simp only [eval_mul, eval_C] at he
  have hnval : (Lagrange.nodal (Finset.range (N + 1)) (node N)).eval x =
      (extremaNodal ⟨N, by omega⟩).eval x / (extremaNodal ⟨N, by omega⟩).leadingCoeff := by
    apply (eq_div_iff hlc).mpr
    simpa [mul_comm] using he.symm
  have hw := nodalWeight_extrema ⟨N,hNp⟩ ⟨j,hjp⟩
  change Lagrange.nodalWeight (Finset.range (N+1)) (node N) j = _ at hw
  rw [Lagrange.eval_basis_not_at_node (Finset.mem_range.mpr (by omega)) hx,
    hw, hnval, extremaBarycentricWeight]
  field_simp [hlc, hn, (extremaEndpointFactor_pos ⟨N, by omega⟩ ⟨j, by change j < N + 1; omega⟩).ne', sub_ne_zero.mpr hx]

/-- Every pointwise value is below the genuine bounded supremum (ATAP, Chapter 15). -/
theorem extremaLebesgueFunction_le_constant (N : ℕ) {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    extremaLebesgueFunction N x ≤ extremaLebesgueConstant N :=
  le_csSup (bddAbove_extremaLebesgueFunction N) ⟨x, hx, rfl⟩

/-- The order-zero Lebesgue constant is one (ATAP, Chapter 15). -/
@[simp] theorem extremaLebesgueConstant_zero : extremaLebesgueConstant 0 = 1 := by
  apply le_antisymm
  · apply csSup_le (Set.Nonempty.image _ (show (Icc (-1 : ℝ) 1).Nonempty from ⟨0, by norm_num⟩))
    rintro _ ⟨x, hx, rfl⟩
    simp
  · simpa using extremaLebesgueFunction_le_constant 0 (x := 0) (by norm_num)

/-- An explicit logarithmic Lebesgue-constant bound for positive order
(ATAP, Chapter 15, nonsharp constants). -/
theorem extremaLebesgueConstant_le_log {N : ℕ} (hN : N ≠ 0) :
    extremaLebesgueConstant N ≤ 2 * π + 2 + 2 * log N := by
  apply csSup_le (Set.Nonempty.image _ (show (Icc (-1 : ℝ) 1).Nonempty from ⟨0, by norm_num⟩))
  rintro _ ⟨x, hx, rfl⟩
  exact extremaLebesgueFunction_le_log hN hx

/-- Interpolation reproduces every polynomial of degree at most `N`
(ATAP, Chapter 4). -/
theorem extremaInterpolant_eq_self {N : ℕ} {p : ℂ[X]} (hp : p.degree < (N + 1 : ℕ)) :
    extremaInterpolant N (fun x : ℝ ↦ p.eval (x : ℂ)) = p := by
  by_cases h : N = 0
  · subst N
    have hp0 : p.degree ≤ 0 := by
      have hp' : p.degree < (1 : WithBot ℕ) := by simpa using hp
      revert hp'
      cases p.degree
      · exact fun _ ↦ bot_le
      · exact fun h ↦ WithBot.coe_le_coe.mpr (Nat.le_of_lt_succ (WithBot.coe_lt_coe.mp h))
    have he := eq_C_of_degree_le_zero (p := p) hp0
    rw [extremaInterpolant_zero, he]
    simp
  · rw [extremaInterpolant_positive ⟨N,Nat.pos_of_ne_zero h⟩]
    symm
    apply Lagrange.eq_interpolate (v := fun j ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero h⟩ j : ℂ))
      (Complex.ofReal_injective.comp_injOn (strictAntiOn_classicalNode ⟨N,Nat.pos_of_ne_zero h⟩).injOn)
    simpa only [Finset.card_range] using hp

/-- Interpolation commutes with subtraction (ATAP, Chapter 4). -/
theorem extremaInterpolant_sub (N : ℕ) (f g : ℝ → ℂ) :
    extremaInterpolant N (fun x ↦ f x - g x) =
      extremaInterpolant N f - extremaInterpolant N g := by
  by_cases h : N = 0
  · subst N
    simp only [extremaInterpolant_zero, map_sub]
  · simp only [extremaInterpolant_positive ⟨N,Nat.pos_of_ne_zero h⟩,
      Lagrange.classicalInterpolate]
    exact map_sub (Lagrange.interpolate (Finset.range (N+1))
      (fun j ↦ (classicalNode ⟨N,Nat.pos_of_ne_zero h⟩ j : ℂ))) _ _

/-- The pointwise near-best estimate for any uniformly accurate polynomial competitor
(ATAP, Chapter 16). -/
theorem norm_sub_extremaInterpolant_le_competitor {N : ℕ} {f : ℝ → ℂ} {p : ℂ[X]}
    (hp : p.degree < (N + 1 : ℕ)) {ε : ℝ}
    (hf : ∀ x ∈ Icc (-1) 1, ‖f x - p.eval (x : ℂ)‖ ≤ ε)
    {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    ‖f x - (extremaInterpolant N f).eval (x : ℂ)‖ ≤
      (1 + extremaLebesgueFunction N x) * ε := by
  have he := extremaInterpolant_sub N f (fun x : ℝ ↦ p.eval (x : ℂ))
  rw [extremaInterpolant_eq_self hp] at he
  have hi := norm_eval_extremaInterpolant_le N (fun t ↦ f t - p.eval (t : ℂ))
    (fun j ↦ hf _ (extremaGridNode_mem_Icc N j)) x
  rw [he, eval_sub] at hi
  have ht := norm_sub_le (f x - p.eval (x : ℂ))
    ((extremaInterpolant N f).eval (x : ℂ) - p.eval (x : ℂ))
  rw [sub_sub_sub_cancel_right] at ht
  nlinarith [hf x hx]

/-- Uniform near-best estimate using the bounded Lebesgue constant
(ATAP, Chapter 16). -/
theorem norm_sub_extremaInterpolant_le_lebesgueConstant {N : ℕ} {f : ℝ → ℂ} {p : ℂ[X]}
    (hp : p.degree < (N + 1 : ℕ)) {ε : ℝ}
    (hf : ∀ x ∈ Icc (-1) 1, ‖f x - p.eval (x : ℂ)‖ ≤ ε)
    {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    ‖f x - (extremaInterpolant N f).eval (x : ℂ)‖ ≤
      (1 + extremaLebesgueConstant N) * ε := by
  apply (norm_sub_extremaInterpolant_le_competitor hp hf hx).trans
  exact mul_le_mul_of_nonneg_right
    (add_le_add (le_refl 1) (extremaLebesgueFunction_le_constant N hx))
    ((norm_nonneg _).trans (hf x hx))

/-- An explicit logarithmic near-best factor for positive interpolation order
(ATAP, Chapter 16; nonsharp Lebesgue bound). -/
theorem norm_sub_extremaInterpolant_le_log_competitor {N : ℕ} (hN : N ≠ 0)
    {f : ℝ → ℂ} {p : ℂ[X]} (hp : p.degree < (N + 1 : ℕ)) {ε : ℝ}
    (hf : ∀ x ∈ Icc (-1) 1, ‖f x - p.eval (x : ℂ)‖ ≤ ε)
    {x : ℝ} (hx : x ∈ Icc (-1) 1) :
    ‖f x - (extremaInterpolant N f).eval (x : ℂ)‖ ≤
      (2 * Real.pi + 3 + 2 * Real.log N) * ε := by
  apply (norm_sub_extremaInterpolant_le_lebesgueConstant hp hf hx).trans
  apply mul_le_mul_of_nonneg_right _ ((norm_nonneg _).trans (hf x hx))
  linarith [extremaLebesgueConstant_le_log hN]