/-- The first-order case is multiplication plus multiplication after one derivative. -/
theorem differentialOperator_one_apply (a : ℕ → R[X]) (p : R[X]) :
    differentialOperator a 1 p = a 0 * p + a 1 * derivative p := by
  simp [differentialOperator_apply, Finset.sum_range_succ]

/-- Coefficient assembly reconstructs the exact polynomial residual. -/
theorem repr_symm_differentialOperator (b c : Module.Basis ℕ R R[X])
    (a : ℕ → R[X]) (N : ℕ) (v : ℕ →₀ R) :
    c.repr.symm (b.coefficientMap c (differentialOperator a N) v) =
      ∑ k ∈ Finset.range (N + 1), a k * derivative^[k] (b.repr.symm v) := by
  rw [b.repr_symm_coefficientMap, differentialOperator_apply]

/-- Uniform bandwidth bounds on the assembled terms bound their sum. In the
ultraspherical method these term bounds come from derivative, conversion, and
multiplication bandwidths, not from an arbitrary square truncation. -/
theorem differentialOperator_isBanded (b c : Module.Basis ℕ R R[X])
    (a : ℕ → R[X]) (N l u : ℕ)
    (h : ∀ k ≤ N, (c.toMatrix (fun j ↦ a k * derivative^[k] (b j))).IsBanded l u) :
    (c.toMatrix (fun j ↦ differentialOperator a N (b j))).IsBanded l u := by
  intro i j hij
  simp only [Module.Basis.toMatrix_apply, differentialOperator_apply, map_sum, Finsupp.finsetSum_apply]
  apply Finset.sum_eq_zero
  intro k hk
  exact h k (Nat.lt_succ_iff.mp (Finset.mem_range.mp hk)) i j hij

/-- A term `p D^k` represented at level N has the derivative/conversion bandwidth
plus the polynomial multiplication bandwidth. -/
theorem ultrasphericalBasis_differential_term_isBanded {k N : ℕ} (hN : 1 ≤ N) (hk : k ≤ N)
    (p : ℝ[X]) :
    ((ultrasphericalBasis N).toMatrix
      (fun j ↦ p * derivative^[k] (chebyshevBasis j))).IsBanded
        p.natDegree (p.natDegree + 2 * N - k) := by
  have h := chebyshevBasis.toMatrix_comp_isBanded (ultrasphericalBasis N)
    (ultrasphericalBasis N) (derivative ^ k) (LinearMap.mulLeft ℝ p)
    (by simpa only [Module.End.pow_apply] using ultrasphericalBasis_derivative_to_level_isBanded hk)
    (ultrasphericalBasis_mul_isBanded hN p)
  have he : (2 * N - k) + p.natDegree = p.natDegree + 2 * N - k := by omega
  simpa only [LinearMap.mulLeft_apply, Module.End.pow_apply, Nat.zero_add, he] using h

/-- Full all-order assembly, with separate bounds for each coefficient degree.
The analytic existence or uniqueness of an ODE solution is not assumed or asserted. -/
theorem ultraspherical_differentialOperator_isBanded
    (a : ℕ → ℝ[X]) {N l u : ℕ} (hN : 1 ≤ N)
    (hl : ∀ k ≤ N, (a k).natDegree ≤ l)
    (hu : ∀ k ≤ N, (a k).natDegree + 2 * N - k ≤ u) :
    ((ultrasphericalBasis N).toMatrix
      (fun j ↦ differentialOperator a N (chebyshevBasis j))).IsBanded l u := by
  apply differentialOperator_isBanded
  intro k hk
  exact (ultrasphericalBasis_differential_term_isBanded hN hk (a k)).mono (hl k hk) (hu k hk)

/-- A common coefficient-degree bound d gives residual bandwidths (d,d+2N). -/
theorem ultraspherical_differentialOperator_isBanded_of_degree_le
    (a : ℕ → ℝ[X]) {N d : ℕ} (hN : 1 ≤ N)
    (hd : ∀ k ≤ N, (a k).natDegree ≤ d) :
    ((ultrasphericalBasis N).toMatrix
      (fun j ↦ differentialOperator a N (chebyshevBasis j))).IsBanded d (d + 2 * N) := by
  apply ultraspherical_differentialOperator_isBanded a hN hd
  intro k hk
  have := hd k hk
  omega

/-- With m input coefficients and coefficient degrees at most d, m+d residual
rows capture the whole polynomial result. This justifies rectangular padding. -/
theorem ultraspherical_differentialOperator_coeff_eq_zero
    (a : ℕ → ℝ[X]) {N d m : ℕ} (hN : 1 ≤ N)
    (hd : ∀ k ≤ N, (a k).natDegree ≤ d) (v : ℕ →₀ ℝ)
    (hv : ∀ j ∈ v.support, j < m) {i : ℕ} (hi : m + d ≤ i) :
    chebyshevBasis.coefficientMap (ultrasphericalBasis N) (differentialOperator a N) v i = 0 :=
  chebyshevBasis.coefficientMap_eq_zero_of_bandwidth (ultrasphericalBasis N)
    (differentialOperator a N) (ultraspherical_differentialOperator_isBanded_of_degree_le a hN hd)
    v hv hi

/-- Evaluation at the right endpoint is the row of ones. -/
theorem chebyshevBasis_eval_one (v : ℕ →₀ ℝ) :
    (chebyshevBasis.repr.symm v).eval 1 = ∑ j ∈ v.support, v j := by
  simpa [chebyshevBasis_apply, Chebyshev.T_eval_one] using
    chebyshevBasis_functional_apply (leval (1 : ℝ)) v

/-- Evaluation at the left endpoint is the alternating row. -/
theorem chebyshevBasis_eval_neg_one (v : ℕ →₀ ℝ) :
    (chebyshevBasis.repr.symm v).eval (-1) =
      ∑ j ∈ v.support, v j * (-1 : ℝ) ^ j := by
  simpa [chebyshevBasis_apply, Chebyshev.T_eval_neg_one] using
    chebyshevBasis_functional_apply (leval (-1 : ℝ)) v

/-- The right endpoint derivative row is the sequence of squared degrees. -/
theorem chebyshevBasis_derivative_eval_one (v : ℕ →₀ ℝ) :
    (derivative (chebyshevBasis.repr.symm v)).eval 1 =
      ∑ j ∈ v.support, v j * (j : ℝ) ^ 2 := by
  simpa [chebyshevBasis_apply, Chebyshev.derivative_T_eval_one] using
    chebyshevBasis_functional_apply ((leval (1 : ℝ)).comp derivative) v

/-- The left endpoint derivative row includes the zero constant entry. -/
theorem chebyshevBasis_derivative_eval_neg_one (v : ℕ →₀ ℝ) :
    (derivative (chebyshevBasis.repr.symm v)).eval (-1) =
      ∑ j ∈ v.support, v j * ((-1 : ℝ) ^ (j + 1) * (j : ℝ) ^ 2) := by
  have h (j : ℕ) : (derivative (chebyshevBasis j)).eval (-1) =
      (-1 : ℝ) ^ (j + 1) * (j : ℝ) ^ 2 := by
    cases j with
    | zero => simp [chebyshevBasis_apply, Chebyshev.T_zero]
    | succ j =>
      rw [chebyshevBasis_apply, Chebyshev.derivative_T_succ_eval_neg_one]
      simp [pow_succ, Nat.cast_add, Nat.cast_one]
  simpa only [LinearMap.comp_apply, leval_apply, h] using
    chebyshevBasis_functional_apply ((leval (-1 : ℝ)).comp derivative) v

/-- A bordered coefficient equation is equivalent to the original boundary
and residual equations; it makes no invertibility assumption. -/
theorem borderedCoefficientMap_eq_iff (b : Basis ι R M) (c : Basis κ R N)
    {q : ℕ} (B : M →ₗ[R] (Fin q → R)) (L : M →ₗ[R] N)
    (v : ι →₀ R) (β : Fin q → R) (f : N) :
    b.borderedCoefficientMap c B L v = (β, c.repr f) ↔
      B (b.repr.symm v) = β ∧ L (b.repr.symm v) = f := by
  rw [borderedCoefficientMap_apply, Prod.mk.injEq, c.repr.injective.eq_iff]

/-- Termwise application represents the actual differential action. -/
theorem hasSum_differentialOperator_series {P : ℕ → ℂ[X]} {a : ℕ → ℂ}
    {u : ℕ → ℝ} {U : TopologicalSpace.Opens ℂ} (hu : Summable u)
    (hb : ∀ n z, z ∈ U → ‖a n * (P n).eval z‖ ≤ u n)
    (A : ℕ → ℂ[X]) (m : ℕ) (f : IteratedDerivInput U m)
    (hf : ∀ z : U, HasSum (fun n ↦ a n * (P n).eval z.val) (f.val z)) (z : U) :
    HasSum (fun n ↦ a n * (differentialOperator A m (P n)).eval z.val)
      (differentialOperatorOnFunctions A m f.val (.ofIterated f z)) := by
  rw [differentialOperatorOnFunctions_eq_extension A m f
    (fun z ↦ ∑' n, a n * (P n).eval z) (analyticOnNhd_polynomial_tsum hu hb)
    (fun z ↦ (hf z).tsum_eq)]
  exact hasSum_differentialOperator_series_tsum U.isOpen hu hb A m z.property

/-- Absolute coefficient reconstruction agrees with the classical differential action. -/
theorem hasSum_seriesCoefficientMap_differentialOperator
    (b c : Module.Basis ℕ ℂ ℂ[X]) {a : ℕ → ℂ} {u : ℕ → ℝ}
    {U : TopologicalSpace.Opens ℂ} (hu : Summable u)
    (hb : ∀ n z, z ∈ U → ‖a n * (b n).eval z‖ ≤ u n)
    (A : ℕ → ℂ[X]) (m : ℕ) (f : IteratedDerivInput U m)
    (hf : ∀ z : U, HasSum (fun n ↦ a n * (b n).eval z.val) (f.val z)) (z : U)
    (habs : Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * c.repr (differentialOperator A m (b p.1)) p.2 * (c p.2).eval z.val‖)) :
    HasSum (fun i ↦ regroupCoefficients (RegroupCoefficientInput.ofSummableNorm habs Prod.snd i))
      (differentialOperatorOnFunctions A m f.val (.ofIterated f z)) := by
  rw [differentialOperatorOnFunctions_eq_extension A m f
    (fun z ↦ ∑' n, a n * (b n).eval z) (analyticOnNhd_polynomial_tsum hu hb)
    (fun z ↦ (hf z).tsum_eq)]
  exact hasSum_seriesCoefficientMap_differentialOperator_tsum b c U.isOpen hu hb A m z.property habs

/-- The reconstruction also respects coefficient-ring evaluation maps. -/
theorem hasSum_seriesCoefficientMap_differentialOperator_eval₂
    {R : Type*} [CommRing R] (b c : Module.Basis ℕ R R[X]) (σ : R →+* ℂ)
    {a : ℕ → ℂ} {u : ℕ → ℝ} {U : TopologicalSpace.Opens ℂ} (hu : Summable u)
    (hb : ∀ n z, z ∈ U → ‖a n * (b n).eval₂ σ z‖ ≤ u n)
    (A : ℕ → R[X]) (m : ℕ) (f : IteratedDerivInput U m)
    (hf : ∀ z : U, HasSum (fun n ↦ a n * (b n).eval₂ σ z.val) (f.val z)) (z : U)
    (habs : Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * σ (c.repr (differentialOperator A m (b p.1)) p.2) *
        (c p.2).eval₂ σ z.val‖)) :
    HasSum (fun i ↦ regroupCoefficients (RegroupCoefficientInput.ofSummableNorm habs Prod.snd i))
      (differentialOperatorOnFunctions (fun k ↦ (A k).map σ) m f.val (.ofIterated f z)) := by
  have hb' : ∀ n z, z ∈ U → ‖a n * ((b n).map σ).eval z‖ ≤ u n := by
    simpa only [eval_map] using hb
  rw [differentialOperatorOnFunctions_eq_extension _ m f
    (fun z ↦ ∑' n, a n * ((b n).map σ).eval z) (analyticOnNhd_polynomial_tsum hu hb')
    (fun z ↦ by simpa only [eval_map] using (hf z).tsum_eq)]
  simpa only [eval_map] using
    hasSum_seriesCoefficientMap_differentialOperator_eval₂_tsum b c σ U.isOpen hu hb A m z.property habs

/-- Convergent coefficient rows give the coordinate form of analytic residual
reconstruction. Weighted contributions alone need not have convergent rows. -/
theorem hasSum_seriesCoefficientMap_differentialOperator_eval₂_of_rows
    {R : Type*} [CommRing R] (b c : Module.Basis ℕ R R[X]) (σ : R →+* ℂ)
    {a : ℕ → ℂ} {u : ℕ → ℝ} {U : TopologicalSpace.Opens ℂ} (hu : Summable u)
    (hb : ∀ n z, z ∈ U → ‖a n * (b n).eval₂ σ z‖ ≤ u n)
    (A : ℕ → R[X]) (m : ℕ) (f : IteratedDerivInput U m)
    (hf : ∀ z : U, HasSum (fun n ↦ a n * (b n).eval₂ σ z.val) (f.val z)) (z : U)
    (habs : Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * σ (c.repr (differentialOperator A m (b p.1)) p.2) *
        (c p.2).eval₂ σ z.val‖))
    (rows : ∀ i, b.SeriesCoefficientInput c (differentialOperator A m) σ i)
    (hrows : ∀ i, (rows i).val = a) :
    HasSum (fun i ↦ b.seriesCoefficientMap c (differentialOperator A m) σ (rows i) *
      (c i).eval₂ σ z.val)
      (differentialOperatorOnFunctions (fun k ↦ (A k).map σ) m f.val (.ofIterated f z)) := by
  have hb' : ∀ n z, z ∈ U → ‖a n * ((b n).map σ).eval z‖ ≤ u n := by
    simpa only [eval_map] using hb
  have hf' : ∀ z : U, HasSum (fun n ↦ a n * ((b n).map σ).eval z.val) (f.val z) := by
    simpa only [eval_map] using hf
  have hs := hasSum_differentialOperator_series hu hb' (fun k ↦ (A k).map σ) m f hf' z
  simp only [← differentialOperator_map, eval_map] at hs
  have h := b.hasSum_seriesCoefficientMap_of_rows c (differentialOperator A m) σ a z.val habs rows hrows
  rw [hs.tsum_eq] at h
  exact h

/-- The normally convergent expansion's ultraspherical residual represents its
actual differential action, also for order-zero multiplication. -/
theorem hasSum_ultraspherical_differentialOperator_series
    {a : ℕ → ℂ} {u : ℕ → ℝ} {U : TopologicalSpace.Opens ℂ} (hu : Summable u)
    (hb : ∀ n z, z ∈ U → ‖a n * (Chebyshev.T ℂ n).eval z‖ ≤ u n)
    (A : ℕ → ℝ[X]) (m : ℕ) (f : IteratedDerivInput U m)
    (hf : ∀ z : U, HasSum (fun n ↦ a n * (Chebyshev.T ℂ n).eval z.val) (f.val z))
    (z : U) (w : ℕ → ℝ)
    (hw : ∀ i, ‖(ultrasphericalBasis m i).eval₂ ofRealHom z.val‖ ≤ w i)
    (ha : Summable (fun j ↦ ‖a j‖ *
      ∑ i ∈ ((ultrasphericalBasis m).repr (differentialOperator A m (chebyshevBasis j))).support,
        ‖(((ultrasphericalBasis m).repr (differentialOperator A m (chebyshevBasis j))) i : ℂ)‖ * w i)) :
    HasSum (fun i ↦ chebyshevBasis.seriesCoefficientMap (ultrasphericalBasis m)
      (differentialOperator A m) ofRealHom (ultrasphericalDifferentialSeriesInput A m a i) *
      (ultrasphericalBasis m i).eval₂ ofRealHom z.val)
      (differentialOperatorOnFunctions (fun k ↦ (A k).map ofRealHom) m f.val (.ofIterated f z)) := by
  have he (n : ℕ) (v : ℂ) : (chebyshevBasis n).eval₂ ofRealHom v =
      (Chebyshev.T ℂ n).eval v := by
    rw [chebyshevBasis_apply, eval₂_eq_eval_map, Chebyshev.map_T]
  exact hasSum_seriesCoefficientMap_differentialOperator_eval₂_of_rows chebyshevBasis
    (ultrasphericalBasis m) ofRealHom hu (by simpa only [he] using hb) A m f
    (by simpa only [he] using hf) z
    (chebyshevBasis.summable_norm_seriesCoefficientMap_terms (ultrasphericalBasis m)
      (differentialOperator A m) ofRealHom a z.val w hw ha)
    (ultrasphericalDifferentialSeriesInput A m a) (fun _ ↦ rfl)

/-- Boundary errors involve certified derivative values at exactly the used points. -/
theorem norm_boundary_complexChebyshevSeries_sub_sum_le
    {a : ℕ → ℂ} {C ρ r : ℝ} (hρ : 1 < ρ) (hr : 1 < r) (hrρ : r < ρ)
    (ha : ∀ n, ‖a n‖ ≤ C / ρ ^ n)
    {ι : Type*} (s : Finset ι) (c : ι → ℂ) (k : ι → ℕ) (z : ι → ℂ) (δ : ι → ℝ)
    (hz : ∀ i ∈ s, z i ∈ bernsteinEllipse ⟨r, hr⟩) (hδ : ∀ i ∈ s, 0 < δ i)
    (hball : ∀ i ∈ s, Metric.sphere (z i) (δ i) ⊆ bernsteinEllipse ⟨r, hr⟩) (N : ℕ) :
    ‖(∑ i : s, c i.val * iteratedDerivOnOpen (complexChebyshevSeriesDerivInput hρ ha (k i.val))
        ⟨z i.val, (closure_bernsteinEllipse_subset hr hrρ) (subset_closure (hz i.val i.property))⟩) -
      ∑ i ∈ s, c i * ∑ n ∈ Finset.range N, a n * (derivative^[k i] (T ℂ n)).eval (z i)‖ ≤
      (∑ i ∈ s, ‖c i‖ * (k i).factorial / (δ i) ^ k i) *
        (C * (r / ρ) ^ N / (1 - r / ρ)) := by
  rw [← s.sum_coe_sort (fun i ↦ c i * ∑ n ∈ Finset.range N,
    a n * (derivative^[k i] (T ℂ n)).eval (z i)),
    ← s.sum_coe_sort (fun i ↦ ‖c i‖ * (k i).factorial / (δ i) ^ k i),
    ← Finset.sum_sub_distrib, Finset.sum_mul]
  apply (norm_sum_le _ _).trans
  apply Finset.sum_le_sum
  intro i _
  rw [← mul_sub, norm_mul]
  calc
    _ ≤ ‖c i.val‖ * ((k i.val).factorial * (C * (r / ρ) ^ N / (1 - r / ρ)) /
        (δ i.val) ^ k i.val) :=
      mul_le_mul_of_nonneg_left
        (norm_iterated_derivative_complexChebyshevSeries_sub_sum_le hρ hr hrρ ha (k i.val) _
          (hz i.val i.property) (hδ i.val i.property) (hball i.val i.property) N) (norm_nonneg _)
    _ = _ := by ring