/-- Columns are the target coordinates of images of source basis vectors
(Olver--Townsend 2013, §2). -/
theorem coefficientMap_single (b : Basis ι R M) (c : Basis κ R N)
    (L : M →ₗ[R] N) (j : ι) :
    b.coefficientMap c L (Finsupp.single j 1) = c.repr (L (b j)) := by
  simp [coefficientMap]

/-- Matrix entries act by a finite sum over the input support
(Olver--Townsend 2013, §2). -/
theorem coefficientMap_apply (b : Basis ι R M) (c : Basis κ R N)
    (L : M →ₗ[R] N) (v : ι →₀ R) (i : κ) :
    b.coefficientMap c L v i =
      v.sum (fun j a ↦ c.toMatrix (fun j ↦ L (b j)) i j * a) := by
  classical
  conv_lhs => rw [← v.sum_single]
  simp only [map_finsuppSum, Finsupp.sum_apply]
  simp only [coefficientMap, LinearMap.comp_apply, LinearEquiv.coe_coe,
    Finsupp.sum, map_smul,
    b.repr_symm_single, Finsupp.coe_smul, Pi.smul_apply,
    smul_eq_mul, toMatrix_apply, mul_comm]

/-- Composition is exact on finitely supported coefficients, with no infinite
matrix multiplication (Olver--Townsend 2013, §2). -/
theorem coefficientMap_comp (b : Basis ι R M) (c : Basis κ R N) (d : Basis τ R P)
    (L : M →ₗ[R] N) (K : N →ₗ[R] P) :
    (c.coefficientMap d K).comp (b.coefficientMap c L) =
      b.coefficientMap d (K.comp L) := by
  ext v i
  simp [coefficientMap]

/-- With lower bandwidth `l`, the first `N+l` rows retain every output coefficient
of an input supported in the first `N` columns. No square-section assumption. -/
theorem coefficientMap_eq_zero_of_bandwidth (b : Basis ℕ R M) (c : Basis ℕ R P)
    (L : M →ₗ[R] P) {l u : ℕ}
    (hL : (c.toMatrix (fun j ↦ L (b j))).IsBanded l u)
    (v : ℕ →₀ R) {N : ℕ} (hv : ∀ j ∈ v.support, j < N)
    {i : ℕ} (hi : N + l ≤ i) : b.coefficientMap c L v i = 0 := by
  rw [coefficientMap_apply_range b c L v hv]
  apply Finset.sum_eq_zero
  intro j hj
  rw [hL i j (Or.inl (by have := Finset.mem_range.mp hj; omega)), zero_mul]

/-- The existing same-basis coefficient map is conjugation by the basis equivalence. -/
theorem coefficientMap_eq_conj {ι : Type*} (b : Module.Basis ι R R[X])
    (L : Module.End R R[X]) : b.coefficientMap b L = b.repr.conjAlgEquiv R L := by
  rfl

/-- In any polynomial basis, multiplication coefficients satisfy `M[p]=p(J)`.
The operator acts on finitely supported coefficients, even for an infinite basis. -/
theorem coefficientMap_mulLeft_eq_aeval {ι : Type*} (b : Module.Basis ι R R[X]) (p : R[X]) :
    b.coefficientMap b (LinearMap.mulLeft R p) =
      aeval (b.coefficientMap b (LinearMap.mulLeft R (X : R[X]))) p := by
  rw [coefficientMap_eq_conj, coefficientMap_eq_conj,
    aeval_algHom_apply, aeval_mulLeft_X]

/-- Polynomial functional calculus reconstructs exact multiplication of every
finite coefficient vector (Olver--Townsend, Sections 2--3). -/
theorem repr_symm_aeval_mulLeft {ι : Type*} (b : Module.Basis ι R R[X])
    (p : R[X]) (v : ι →₀ R) :
    b.repr.symm ((aeval (b.coefficientMap b (LinearMap.mulLeft R (X : R[X]))) p) v) =
      p * b.repr.symm v := by
  rw [← coefficientMap_mulLeft_eq_aeval, b.repr_symm_coefficientMap]
  rfl

/-- Clenshaw evaluation commutes with a ring homomorphism, in the full-constant
convention (Clenshaw 1955). -/
theorem clenshaw_map {S : Type*} [Ring S] (f : R →+* S)
    (a b c : ℕ → R) (p₀ p₁ : R) (n : ℕ) :
    clenshaw (f ∘ a) (f ∘ b) (f ∘ c) (f p₀) (f p₁) n =
      f (clenshaw a b c p₀ p₁ n) := by
  simp [clenshaw, clenshawState_map]

/-- Polynomial-valued Clenshaw reconstructs the full Chebyshev expansion. -/
theorem clenshaw_polynomial_T (c : ℕ → R) (n : ℕ) :
    clenshaw (fun _ ↦ 2 * (X : R[X])) (fun _ ↦ -1) (fun k ↦ Polynomial.C (c k)) 1 X n =
      ∑ k ∈ Finset.range n, Polynomial.C (c k) * T R k := by
  have hp (k : ℕ) : T R (k + 2 : ℕ) =
      (2 * X) * T R (k + 1 : ℕ) + (-1) * T R k := by
    simp only [Nat.cast_add, Nat.cast_one, Nat.cast_ofNat, T_add_two, neg_one_mul,
      sub_eq_add_neg]
  simpa only [Nat.cast_zero, Nat.cast_one, T_zero, T_one] using
    clenshaw_eq_sum (fun _ ↦ 2 * (X : R[X])) (fun _ ↦ -1)
      (fun k ↦ Polynomial.C (c k)) (fun k ↦ T R k) hp n

/-- The same backward recurrence evaluates a Chebyshev expansion at one element
of an arbitrary algebra. This includes linear operators without assuming that
the whole endomorphism algebra is commutative. -/
theorem clenshaw_aeval_T {A : Type*} [Ring A] [Algebra R A] (J : A)
    (c : ℕ → R) (n : ℕ) :
    clenshaw (fun _ ↦ 2 * J) (fun _ ↦ -1) (fun k ↦ algebraMap R A (c k)) 1 J n =
      aeval J (∑ k ∈ Finset.range n, Polynomial.C (c k) * T R k) := by
  have h := clenshaw_map (aeval J).toRingHom (fun _ ↦ 2 * (X : R[X])) (fun _ ↦ -1)
    (fun k ↦ Polynomial.C (c k)) 1 X n
  rw [clenshaw_polynomial_T] at h
  simpa only [Function.comp_def, AlgHom.toRingHom_eq_coe, RingHom.coe_coe,
    map_mul, map_ofNat, map_neg, map_one, aeval_X, Polynomial.aeval_C] using h

/-- Operator Clenshaw computes multiplication by the Chebyshev expansion exactly
on the polynomial module (Olver--Townsend, Sections 2--3). -/
theorem clenshaw_mulLeft_T (c : ℕ → R) (n : ℕ) :
    clenshaw (fun _ ↦ 2 * LinearMap.mulLeft R (X : R[X])) (fun _ ↦ -1)
      (fun k ↦ algebraMap R (Module.End R R[X]) (c k)) 1
      (LinearMap.mulLeft R (X : R[X])) n =
      LinearMap.mulLeft R (∑ k ∈ Finset.range n, Polynomial.C (c k) * T R k) := by
  rw [clenshaw_aeval_T, aeval_mulLeft_X]

/-- In any polynomial basis, operator Clenshaw produces the multiplication
coefficient map on finitely supported vectors, with no finite-section claim. -/
theorem clenshaw_coefficientMap_T {ι : Type*} (b : Module.Basis ι R R[X])
    (c : ℕ → R) (n : ℕ) :
    clenshaw (fun _ ↦ 2 * b.coefficientMap b (LinearMap.mulLeft R (X : R[X])))
      (fun _ ↦ -1) (fun k ↦ algebraMap R (Module.End R (ι →₀ R)) (c k)) 1
      (b.coefficientMap b (LinearMap.mulLeft R (X : R[X]))) n =
      b.coefficientMap b
        (LinearMap.mulLeft R (∑ k ∈ Finset.range n, Polynomial.C (c k) * T R k)) := by
  rw [clenshaw_aeval_T, ← coefficientMap_mulLeft_eq_aeval]

/-- Absolute convergence permits regrouping weighted contributions, including
zero output modes, without asserting convergence of unweighted rows. -/
theorem hasSum_seriesCoefficientMap (b c : Basis ℕ R R[X])
    (L : R[X] →ₗ[R] R[X]) (σ : R →+* ℂ) (a : ℕ → ℂ) (z : ℂ)
    (habs : Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * σ (c.repr (L (b p.1)) p.2) * (c p.2).eval₂ σ z‖)) :
    HasSum (fun i ↦ regroupCoefficients
      (RegroupCoefficientInput.ofSummableNorm habs Prod.snd i))
      (∑' j, a j * (L (b j)).eval₂ σ z) := by
  have hs := habs.of_norm
  have hcol (j : ℕ) : (∑' i, a j * σ (c.repr (L (b j)) i) * (c i).eval₂ σ z) =
      a j * (L (b j)).eval₂ σ z := by
    simpa only [mul_assoc] using ((c.hasSum_repr_eval₂ σ (L (b j)) z).mul_left (a j)).tsum_eq
  have he : (∑' p : ℕ × ℕ, a p.1 * σ (c.repr (L (b p.1)) p.2) * (c p.2).eval₂ σ z) =
      ∑' j, a j * (L (b j)).eval₂ σ z := by
    rw [hs.tsum_prod]
    exact tsum_congr hcol
  simpa only [he, regroupCoefficients, RegroupCoefficientInput.ofSummableNorm,
    RegroupCoefficientInput.ofSummable, Set.preimage_singleton, Set.coe_ofPred] using hs.hasSum.tsum_fiberwise Prod.snd

/-- Certified coefficient rows can be factored from their output modes. -/
theorem hasSum_seriesCoefficientMap_of_rows (b c : Basis ℕ R R[X])
    (L : R[X] →ₗ[R] R[X]) (σ : R →+* ℂ) (a : ℕ → ℂ) (z : ℂ)
    (habs : Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * σ (c.repr (L (b p.1)) p.2) * (c p.2).eval₂ σ z‖))
    (rows : ∀ i, SeriesCoefficientInput b c L σ i)
    (hrows : ∀ i, (rows i).val = a) :
    HasSum (fun i ↦ b.seriesCoefficientMap c L σ (rows i) * (c i).eval₂ σ z)
      (∑' j, a j * (L (b j)).eval₂ σ z) := by
  have hs := habs.of_norm
  have hcol (j : ℕ) : (∑' i, a j * σ (c.repr (L (b j)) i) * (c i).eval₂ σ z) =
      a j * (L (b j)).eval₂ σ z := by
    simpa only [mul_assoc] using ((c.hasSum_repr_eval₂ σ (L (b j)) z).mul_left (a j)).tsum_eq
  have he : (∑' p : ℕ × ℕ, a p.1 * σ (c.repr (L (b p.1)) p.2) * (c p.2).eval₂ σ z) =
      ∑' j, a j * (L (b j)).eval₂ σ z := by
    rw [hs.tsum_prod]
    exact tsum_congr hcol
  have h := hs.prod_symm.prod.hasSum
  simp only [Prod.swap_prod_mk] at h
  have he' : (∑' i, ∑' j, a j * σ (c.repr (L (b j)) i) * (c i).eval₂ σ z) =
      ∑' j, a j * (L (b j)).eval₂ σ z := by
    rw [hs.tsum_comm]
    exact tsum_congr hcol
  rw [he'] at h
  simpa only [seriesCoefficientMap, hrows, tsum_mul_right] using h

/-- A single weighted coefficient sum suffices for absolute convergence of
matrix reconstruction; every polynomial matrix column has finite support. -/
theorem summable_norm_seriesCoefficientMap_terms (b c : Basis ℕ R R[X])
    (L : R[X] →ₗ[R] R[X]) (σ : R →+* ℂ) (a : ℕ → ℂ) (z : ℂ) (w : ℕ → ℝ)
    (hw : ∀ i, ‖(c i).eval₂ σ z‖ ≤ w i)
    (ha : Summable (fun j ↦ ‖a j‖ * ∑ i ∈ (c.repr (L (b j))).support,
      ‖σ (c.repr (L (b j)) i)‖ * w i)) :
    Summable (fun p : ℕ × ℕ ↦
      ‖a p.1 * σ (c.repr (L (b p.1)) p.2) * (c p.2).eval₂ σ z‖) := by
  classical
  have hw0 (i : ℕ) : 0 ≤ w i := (norm_nonneg _).trans (hw i)
  have hrow (j : ℕ) : HasSum (fun i ↦ ‖a j‖ * ‖σ (c.repr (L (b j)) i)‖ * w i)
      (‖a j‖ * ∑ i ∈ (c.repr (L (b j))).support, ‖σ (c.repr (L (b j)) i)‖ * w i) := by
    rw [Finset.mul_sum]
    simp only [← mul_assoc]
    apply hasSum_sum_of_ne_finset_zero
    intro i hi
    rw [Finsupp.notMem_support_iff.mp hi, map_zero, norm_zero, mul_zero, zero_mul]
  have hs : Summable (fun p : ℕ × ℕ ↦ ‖a p.1‖ *
      ‖σ (c.repr (L (b p.1)) p.2)‖ * w p.2) := by
    apply (summable_prod_of_nonneg (fun p ↦
      mul_nonneg (mul_nonneg (norm_nonneg _) (norm_nonneg _)) (hw0 p.2))).2
    refine ⟨fun j ↦ (hrow j).summable, ?_⟩
    simpa only [(hrow _).tsum_eq] using ha
  apply Summable.of_nonneg_of_le (fun _ ↦ norm_nonneg _) _ hs
  intro p
  simp only [norm_mul]
  exact mul_le_mul_of_nonneg_left (hw p.2) (by positivity)