/-- Jacobi squared norms are positive throughout the orthogonality range. -/
theorem integral_sq_jacobi_pos {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (n : ℕ) :
    0 < ∫ x, (jacobi α β n).eval x ^ 2 ∂jacobiMeasure α β := by
  simpa only [jacobiSequence_apply] using
    hasPositiveDefiniteIntegral_jacobiMeasure hα hβ _ ((jacobiSequence hα hβ).ne_zero n)

/-- The global reciprocal-norm constant is positive for every integrable Jacobi weight. -/
theorem jacobiReciprocalNormBound_pos {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (k : ℕ) :
    0 < jacobiReciprocalNormBound ⟨α, hα⟩ ⟨β, hβ⟩ k := by
  have h0 := integral_sq_jacobi_pos hα hβ 0
  have h1 := integral_sq_jacobi_pos hα hβ 1
  have h2 := integral_sq_jacobi_pos hα hβ 2
  unfold jacobiReciprocalNormBound
  positivity

/-- The kernel-growth constant is nonnegative. -/
theorem jacobiKernelBound_nonneg {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (k : ℕ) :
    0 ≤ jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k :=
  mul_nonneg (sq_nonneg _) (jacobiReciprocalNormBound_pos hα hβ k).le

/-- Scalar embedding of the classical finite kernel vector (DLMF 18.2(v)). -/
theorem kernelVectorScalar_apply
    (N : {N : ℕ // ∀ i ≤ N, 0 < P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _)}) (x : ℝ) :
    P.kernelVectorScalar 𝕂 hμ N x =
      (@RCLike.ofRealCLM 𝕂 _).compLp (P.kernelVector hμ N x) := rfl

/-- The finite scalar formula on its classical input domain. -/
theorem partialSumValueScalar_apply (μ : Measure ℝ) {N : ℕ}
    (f : IntegralCoefficientFamilyInput 𝕂 μ P N) (x : ℝ) :
    P.partialSumValueScalar 𝕂 μ f x =
      ∑ i : Fin (N + 1), P.integralCoefficientScalar 𝕂 μ (f.at i) * (((P i).eval x : ℝ) : 𝕂) := rfl

/-- The classical complex weighted formula (DLMF 18.18). -/
theorem gegenbauerComplexCoefficient_apply (lam : ℝ) {n : ℕ}
    (f : IntegralCoefficientInput ℂ (gegenbauerMeasure lam) (gegenbauer lam n)) :
    gegenbauerComplexCoefficient lam f =
      ((∫ x, (gegenbauer lam n).eval x ^ 2 ∂gegenbauerMeasure lam : ℝ) : ℂ)⁻¹ *
        ∫ x, f.val x * (((gegenbauer lam n).eval x : ℝ) : ℂ) ∂gegenbauerMeasure lam := rfl

/-- Complex coefficients agree with real coefficients on admissible real data. -/
theorem jacobiComplexCoefficient_ofReal (α β : ℝ) {n : ℕ}
    (f : IntegralCoefficientInput ℝ (jacobiMeasure α β) (jacobi α β n)) :
    jacobiComplexCoefficient α β f.ofReal = (jacobiCoefficient α β f : ℂ) := by
  simp only [jacobiComplexCoefficient, jacobiCoefficient, IntegralCoefficientInput.ofReal_val,
    ofReal_mul, ofReal_inv]
  congr 1
  simp_rw [← ofReal_mul]
  exact integral_complex_ofReal

/-- Complex coefficients agree with real coefficients on admissible real data. -/
theorem gegenbauerComplexCoefficient_ofReal (lam : ℝ) {n : ℕ}
    (f : IntegralCoefficientInput ℝ (gegenbauerMeasure lam) (gegenbauer lam n)) :
    gegenbauerComplexCoefficient lam f.ofReal = (gegenbauerCoefficient lam f : ℂ) := by
  simp only [gegenbauerComplexCoefficient, gegenbauerCoefficient, IntegralCoefficientInput.ofReal_val,
    ofReal_mul, ofReal_inv]
  congr 1
  simp_rw [← ofReal_mul]
  exact integral_complex_ofReal

/-- Orthogonal coefficients minimise the error among all combinations of a
fixed finite orthonormal family (Hilbert-space best approximation). -/
theorem Orthonormal.norm_sub_sum_inner_le {v : ι → E} (hv : Orthonormal 𝕂 v)
    (s : Finset ι) (x : E) (c : ι → 𝕂) :
    ‖x - ∑ i ∈ s, inner 𝕂 (v i) x • v i‖ ≤ ‖x - ∑ i ∈ s, c i • v i‖ := by
  let a := ∑ i ∈ s, inner 𝕂 (v i) x • v i
  let b := ∑ i ∈ s, c i • v i
  have ho : inner 𝕂 (x - a) (a - b) = 0 := by
    rw [inner_sub_right, hv.inner_sub_sum_inner_sum_eq_zero s x,
      hv.inner_sub_sum_inner_sum_eq_zero s x, sub_self]
  have hp := norm_add_sq_eq_norm_sq_add_norm_sq_of_inner_eq_zero (x - a) (a - b) ho
  rw [sub_add_sub_cancel] at hp
  change ‖x - a‖ ≤ ‖x - b‖
  nlinarith [norm_nonneg (x - a), norm_nonneg (x - b), sq_nonneg ‖a - b‖]

/-- A finite orthogonal projection cannot increase the residual norm. -/
theorem Orthonormal.norm_sub_sum_inner_le_norm {v : ι → E} (hv : Orthonormal 𝕂 v)
    (s : Finset ι) (x : E) :
    ‖x - ∑ i ∈ s, inner 𝕂 (v i) x • v i‖ ≤ ‖x‖ := by
  simpa using hv.norm_sub_sum_inner_le s x (fun _ ↦ (0 : 𝕂))

/-- Best approximation holds for every vector in the selected span, without
choosing coordinates for the competitor. -/
theorem Orthonormal.norm_sub_sum_inner_le_of_mem_span {v : ι → E}
    (hv : Orthonormal 𝕂 v) (s : Finset ι) (x y : E)
    (hy : y ∈ Submodule.span 𝕂 (v '' (s : Set ι))) :
    ‖x - ∑ i ∈ s, inner 𝕂 (v i) x • v i‖ ≤ ‖x - y‖ := by
  let a := ∑ i ∈ s, inner 𝕂 (v i) x • v i
  have ho : inner 𝕂 (x - a) y = 0 := by
    induction hy using Submodule.span_induction with
    | mem z hz =>
      obtain ⟨i, hi, rfl⟩ := hz
      apply inner_eq_zero_symm.mp
      rw [inner_sub_right, hv.inner_right_sum _ hi, sub_self]
    | zero => simp
    | add u w _ _ hu hw => rw [inner_add_right, hu, hw, add_zero]
    | smul c u _ hu => rw [inner_smul_right, hu, mul_zero]
  have ha : inner 𝕂 (x - a) (a - y) = 0 := by
    rw [inner_sub_right, hv.inner_sub_sum_inner_sum_eq_zero s x, ho, sub_self]
  have hp := norm_add_sq_eq_norm_sq_add_norm_sq_of_inner_eq_zero (x - a) (a - y) ha
  rw [sub_add_sub_cancel] at hp
  change ‖x - a‖ ≤ ‖x - y‖
  nlinarith [norm_nonneg (x - a), norm_nonneg (x - y), sq_nonneg ‖a - y‖]

/-- The partial-sum polynomial has degree at most its truncation index. -/
theorem degree_partialSum_le {N : ℕ} (f : IntegralCoefficientFamilyInput ℝ μ P N) :
    (P.partialSum μ f).degree ≤ N := by
  apply (degree_sum_le _ _).trans
  apply Finset.sup_le
  intro i hi
  exact (degree_smul_le _ _).trans (by simpa only [P.degree_eq, Nat.cast_le] using
    (Nat.lt_succ_iff.mp i.isLt))

/-- Polynomial partial sums represent finite orthonormal projections in L². -/
theorem polynomialToL2_partialSum (hpos : HasPositiveDefiniteIntegral μ)
    (f : Lp ℝ 2 μ) (N : ℕ) :
    polynomialToL2 hμ (P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P f N)) =
      ∑ i ∈ Finset.range (N + 1), inner ℝ (P.normalizedToL2 hμ ⟨i, squaredNorm_pos hpos P i⟩) f •
        P.normalizedToL2 hμ ⟨i, squaredNorm_pos hpos P i⟩ := by
  simp only [partialSum_ofL2, map_sum, map_smul, inner_normalizedToL2_smul hμ hpos]

/-- Projection reproduces every polynomial of degree at most N. -/
theorem IsOrthogonal.partialSum_polynomialToL2 (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) (N : ℕ) (p : ℝ[X]) (hp : p.degree ≤ N) :
    P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P (polynomialToL2 hμ p) N) = p := by
  apply polynomialToL2_injective hμ hpos
  rw [polynomialToL2_partialSum hμ hpos]
  exact (hP.orthonormal_normalizedToL2 hμ hpos).sum_inner_smul_eq_of_mem_span _ _
    (polynomialToL2_mem_span_normalizedToL2 hμ hpos N p hp)

/-- Subtracting data subtracts their projection polynomials. -/
theorem partialSum_sub (hpos : HasPositiveDefiniteIntegral μ)
    (f g : Lp ℝ 2 μ) (N : ℕ) :
    P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P (f - g : Lp ℝ 2 μ) N) = P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P f N) - P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P g N) := by
  apply polynomialToL2_injective hμ hpos
  simp only [map_sub, polynomialToL2_partialSum hμ hpos, inner_sub_right,
    sub_smul, Finset.sum_sub_distrib]

/-- Weighted polynomial projection is a contraction in L². -/
theorem IsOrthogonal.norm_polynomialToL2_partialSum_le (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) (f : Lp ℝ 2 μ) (N : ℕ) :
    ‖polynomialToL2 hμ (P.partialSum μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P f N))‖ ≤ ‖f‖ := by
  rw [polynomialToL2_partialSum hμ hpos]
  exact (hP.orthonormal_normalizedToL2 hμ hpos).norm_sum_inner_smul_le _ _

/-- A factorial upper bound for a genuine Gamma value. -/
theorem Gamma_nat_add_le_factorial_mul_pow (n k : ℕ) {a : ℝ}
    (ha : a ≤ k + 1) (hna : 2 ≤ (n : ℝ) + a) :
    GammaOn (GammaPoint.ofPos ((n : ℝ) + a) (by linarith)) ≤
      n.factorial * ((n : ℝ) + k) ^ k :=
  Gamma_nat_add_le_factorial_mul_pow_raw n k ha hna

/-- A factorial lower bound for a genuine Gamma value. -/
theorem factorial_le_Gamma_nat_add {n : ℕ} (hn : 1 ≤ n) {a : ℝ} (ha : 1 ≤ a) :
    (n.factorial : ℝ) ≤ GammaOn (GammaPoint.ofPos ((n : ℝ) + a)
      (by have := Nat.cast_nonneg (α := ℝ) n; linarith)) :=
  factorial_le_Gamma_nat_add_raw hn ha

/-- A positive rising-factorial ratio has polynomial growth with a classical constant. -/
theorem ascPochhammer_div_factorial_le_polynomial {a : ℝ} (ha : 0 < a)
    (k : ℕ) (hk : a ≤ k + 1) (n : ℕ) :
    (ascPochhammer ℝ n).eval a / (n.factorial : ℝ) ≤
      (1 + |a| + ((k : ℝ) + 1) ^ k / GammaOn (GammaPoint.ofPos a ha)) *
        ((n : ℝ) + 1) ^ k :=
  ascPochhammer_div_factorial_le_polynomial_raw ha k hk n

/-- When alpha+beta >= -1, the maximum absolute value on [-1,1] occurs at
an endpoint. This follows directly from the differential equation's energy. -/
theorem abs_eval_jacobi_le_max_endpoints {α β : ℝ} (hs : 0 ≤ α + β + 1)
    (n : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    |(jacobi α β n).eval x| ≤
      max |(jacobi α β n).eval (-1)| |(jacobi α β n).eval 1| := by
  by_cases hn : n = 0
  · subst n
    simp
  have hn0 : (0 : ℝ) < n := Nat.cast_pos.mpr (Nat.pos_of_ne_zero hn)
  let c : ℝ := (n : ℝ) * (n + α + β + 1)
  have hc : 0 < c := mul_pos hn0 (by linarith)
  let E := fun t : ℝ ↦ c * (jacobi α β n).eval t ^ 2 +
    (1 - t ^ 2) * (jacobi α β n).derivative.eval t ^ 2
  have hd (t : ℝ) : HasDerivAt E
      (2 * (α - β + (α + β + 1) * t) * (jacobi α β n).derivative.eval t ^ 2) t :=
    hasDerivAt_jacobiSoninEnergy α β n t
  have hcont : Continuous E := continuous_iff_continuousAt.mpr (fun t ↦ (hd t).continuousAt)
  have hlow : c * (jacobi α β n).eval x ^ 2 ≤ E x := by
    dsimp [E]
    have hx2 : x ^ 2 ≤ 1 := by nlinarith [hx.1, hx.2]
    exact le_add_of_nonneg_right (mul_nonneg (by linarith) (sq_nonneg _))
  have hupper : E x ≤ max (c * (jacobi α β n).eval (-1) ^ 2)
      (c * (jacobi α β n).eval 1 ^ 2) := by
    by_cases hsign : 0 ≤ α - β + (α + β + 1) * x
    · have hm : MonotoneOn E (Icc x 1) := by
        apply monotoneOn_of_hasDerivWithinAt_nonneg (convex_Icc x 1) hcont.continuousOn
          (fun t _ ↦ (hd t).hasDerivWithinAt)
        intro t ht
        have ht' := interior_subset ht
        have hst : 0 ≤ α - β + (α + β + 1) * t := by
          nlinarith [mul_nonneg hs (sub_nonneg.mpr ht'.1)]
        positivity
      have h := hm (left_mem_Icc.mpr hx.2) (right_mem_Icc.mpr hx.2) hx.2
      apply h.trans
      simp only [E, one_pow, sub_self, zero_mul, add_zero]
      exact le_max_right _ _
    · have hm : AntitoneOn E (Icc (-1) x) := by
        apply antitoneOn_of_hasDerivWithinAt_nonpos (convex_Icc (-1) x) hcont.continuousOn
          (fun t _ ↦ (hd t).hasDerivWithinAt)
        intro t ht
        have ht' := interior_subset ht
        have hst : α - β + (α + β + 1) * t ≤ 0 := by
          nlinarith [mul_nonneg hs (sub_nonneg.mpr ht'.2)]
        exact mul_nonpos_of_nonpos_of_nonneg (mul_nonpos_of_nonneg_of_nonpos (by norm_num) hst)
          (sq_nonneg _)
      have h := hm (left_mem_Icc.mpr hx.1) (right_mem_Icc.mpr hx.1) hx.1
      apply h.trans
      norm_num only [E, neg_one_sq, sub_self, zero_mul, add_zero]
      exact le_max_left _ _
  have hsq : (jacobi α β n).eval x ^ 2 ≤
      max ((jacobi α β n).eval (-1) ^ 2) ((jacobi α β n).eval 1 ^ 2) := by
    rw [← mul_max_of_nonneg _ _ hc.le] at hupper
    exact (mul_le_mul_iff_right₀ hc).mp (hlow.trans hupper)
  rcases le_total ((jacobi α β n).eval (-1) ^ 2) ((jacobi α β n).eval 1 ^ 2) with h | h
  · rw [max_eq_right h] at hsq
    apply (show |(jacobi α β n).eval x| ≤ |(jacobi α β n).eval 1| by
      nlinarith [sq_abs ((jacobi α β n).eval x), sq_abs ((jacobi α β n).eval 1),
        abs_nonneg ((jacobi α β n).eval x), abs_nonneg ((jacobi α β n).eval 1)]).trans
    exact le_max_right _ _
  · rw [max_eq_left h] at hsq
    apply (show |(jacobi α β n).eval x| ≤ |(jacobi α β n).eval (-1)| by
      nlinarith [sq_abs ((jacobi α β n).eval x), sq_abs ((jacobi α β n).eval (-1)),
        abs_nonneg ((jacobi α β n).eval x), abs_nonneg ((jacobi α β n).eval (-1))]).trans
    exact le_max_left _ _

/-- Throughout the full orthogonality range, one parameter shift bounds the
polynomial by endpoint values of two members of the shifted family. -/
theorem abs_eval_jacobi_succ_le_shifted_endpoints {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (n : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    |(jacobi α β (n + 1)).eval x| ≤
      max |(jacobi (α + 1) β (n + 1)).eval (-1)| |(jacobi (α + 1) β (n + 1)).eval 1| +
      max |(jacobi (α + 1) β n).eval (-1)| |(jacobi (α + 1) β n).eval 1| := by
  let A := max |(jacobi (α + 1) β (n + 1)).eval (-1)| |(jacobi (α + 1) β (n + 1)).eval 1|
  let B := max |(jacobi (α + 1) β n).eval (-1)| |(jacobi (α + 1) β n).eval 1|
  have hA : 0 ≤ A := (abs_nonneg _).trans (le_max_left _ _)
  have hB : 0 ≤ B := (abs_nonneg _).trans (le_max_left _ _)
  have hn : (0 : ℝ) ≤ n := Nat.cast_nonneg n
  have hs : 0 ≤ (α + 1) + β + 1 := by linarith
  have hq := abs_eval_jacobi_le_max_endpoints hs (n + 1) hx
  have hr := abs_eval_jacobi_le_max_endpoints hs n hx
  have he := congrArg (fun p : ℝ[X] ↦ p.eval x) (jacobi_raise_alpha α β n)
  simp only [eval_mul, eval_C, eval_sub] at he
  have hd : 0 < 2 * (n : ℝ) + α + β + 3 := by linarith
  have ha : 0 ≤ (n : ℝ) + α + β + 2 := by linarith
  have hb : 0 ≤ (n : ℝ) + β + 1 := by linarith
  have heabs : (2 * (n : ℝ) + α + β + 3) * |(jacobi α β (n + 1)).eval x| ≤
      ((n : ℝ) + α + β + 2) * A + ((n : ℝ) + β + 1) * B := by
    rw [← abs_of_pos hd, ← abs_mul, he]
    apply (abs_sub _ _).trans
    rw [abs_mul, abs_mul, abs_of_nonneg ha, abs_of_nonneg hb]
    exact add_le_add (mul_le_mul_of_nonneg_left hq ha) (mul_le_mul_of_nonneg_left hr hb)
  have hcmp : ((n : ℝ) + α + β + 2) * A + ((n : ℝ) + β + 1) * B ≤
      (2 * (n : ℝ) + α + β + 3) * (A + B) := by
    nlinarith [mul_nonneg (show 0 ≤ (n : ℝ) + 1 by positivity) hA,
      mul_nonneg (show 0 ≤ (n : ℝ) + α + 2 by linarith) hB]
  exact (mul_le_mul_iff_right₀ hd).mp (heabs.trans hcmp)

/-- The negative endpoint has the same bound with beta replacing alpha. -/
theorem abs_eval_jacobi_neg_one_le_polynomial (α : ℝ) {β : ℝ} (hβ : -1 < β)
    (k : ℕ) (hk : β ≤ k) (n : ℕ) :
    |(jacobi α β n).eval (-1)| ≤ jacobiEndpointBound ⟨β, (GammaPoint.ofPos (β + 1) (by linarith)).property⟩ k * ((n : ℝ) + 1) ^ k := by
  rw [eval_jacobi_neg_one, abs_mul, abs_pow, abs_neg, abs_one, one_pow, one_mul]
  simpa only [eval_jacobi_one] using abs_eval_jacobi_one_le_polynomial hβ α k hk n

/-- Both endpoint values are bounded by the sum of their growth constants. -/
theorem max_abs_eval_jacobi_endpoints_le_polynomial {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (k : ℕ) (hkα : α ≤ k) (hkβ : β ≤ k) (n : ℕ) :
    max |(jacobi α β n).eval (-1)| |(jacobi α β n).eval 1| ≤
      (jacobiEndpointBound ⟨α, (GammaPoint.ofPos (α + 1) (by linarith)).property⟩ k + jacobiEndpointBound ⟨β, (GammaPoint.ofPos (β + 1) (by linarith)).property⟩ k) * ((n : ℝ) + 1) ^ k := by
  apply max_le
  · apply (abs_eval_jacobi_neg_one_le_polynomial α hβ k hkβ n).trans
    gcongr
    exact le_add_of_nonneg_left (jacobiEndpointBound_pos hα k).le
  · apply (abs_eval_jacobi_one_le_polynomial hα β k hkα n).trans
    gcongr
    exact le_add_of_nonneg_right (jacobiEndpointBound_pos hβ k).le

/-- A uniform polynomial growth bound on [-1,1], valid for every alpha,beta > -1.
The integer exponent only needs to dominate alpha+1 and beta. -/
theorem abs_eval_jacobi_le_polynomial {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (k : ℕ) (hkα : α + 1 ≤ k) (hkβ : β ≤ k) (n : ℕ)
    {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    |(jacobi α β n).eval x| ≤
      (1 + 2 * (jacobiEndpointBound ⟨α + 1, (GammaPoint.ofPos (α + 1 + 1) (by linarith)).property⟩ k + jacobiEndpointBound ⟨β, (GammaPoint.ofPos (β + 1) (by linarith)).property⟩ k)) *
        ((n : ℝ) + 1) ^ k := by
  have hα' : -1 < α + 1 := by linarith
  have hC : 0 ≤ jacobiEndpointBound ⟨α + 1, (GammaPoint.ofPos (α + 1 + 1) (by linarith)).property⟩ k + jacobiEndpointBound ⟨β, (GammaPoint.ofPos (β + 1) (by linarith)).property⟩ k :=
    add_nonneg (jacobiEndpointBound_pos hα' k).le (jacobiEndpointBound_pos hβ k).le
  rcases n with _ | n
  · simp only [jacobi_zero, eval_one, abs_one, Nat.cast_zero, zero_add, one_pow, mul_one]
    linarith
  have h := abs_eval_jacobi_succ_le_shifted_endpoints hα hβ n hx
  have h1 := max_abs_eval_jacobi_endpoints_le_polynomial hα' hβ k hkα hkβ (n + 1)
  have h0 := max_abs_eval_jacobi_endpoints_le_polynomial hα' hβ k hkα hkβ n
  apply (h.trans (add_le_add h1 h0)).trans
  have hp : ((n : ℝ) + 1) ^ k ≤ (((n + 1 : ℕ) : ℝ) + 1) ^ k := by
    apply pow_le_pow_left₀ (by positivity)
    push_cast
    linarith
  nlinarith [mul_le_mul_of_nonneg_left hp hC,
    pow_nonneg (show 0 ≤ (((n + 1 : ℕ) : ℝ) + 1) by positivity) k]

/-- A factorial comparison bound on reciprocal norms, away from the first
three degrees. The small degrees are handled separately in the global bound. -/
theorem inv_integral_sq_jacobi_succ_le {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (n k : ℕ) (hn : 2 ≤ n) (hk : α + β + 1 ≤ k) :
    (∫ x, (jacobi α β (n + 1)).eval x ^ 2 ∂jacobiMeasure α β)⁻¹ ≤
      ((2 * (n : ℝ) + α + β + 3) * (n + 1) * ((n : ℝ) + k) ^ k) /
        (2 : ℝ) ^ (α + β + 1) := by
  have hn' : (2 : ℝ) ≤ n := by exact_mod_cast hn
  have hA := factorial_le_Gamma_nat_add (show 1 ≤ n by omega)
    (show 1 ≤ α + 2 by linarith)
  have hB := factorial_le_Gamma_nat_add (show 1 ≤ n by omega)
    (show 1 ≤ β + 2 by linarith)
  have hD := Gamma_nat_add_le_factorial_mul_pow n k
    (show α + β + 2 ≤ k + 1 by linarith) (by linarith : 2 ≤ (n : ℝ) + (α + β + 2))
  simp only [GammaOn, GammaPoint.ofPos] at hA hB hD
  have hAg : 0 < Gamma ((n : ℝ) + α + 2) := Gamma_pos_of_pos (by linarith)
  have hBg : 0 < Gamma ((n : ℝ) + β + 2) := Gamma_pos_of_pos (by linarith)
  have hR : 0 < (2 : ℝ) ^ (α + β + 1) := Real.rpow_pos_of_pos (by norm_num) _
  have hd : 0 ≤ 2 * (n : ℝ) + α + β + 3 := by linarith
  have hF : (0 : ℝ) ≤ n.factorial := Nat.cast_nonneg _
  have hprod : Gamma ((n : ℝ) + α + β + 2) * n.factorial ≤
      ((n : ℝ) + k) ^ k * (Gamma ((n : ℝ) + α + 2) * Gamma ((n : ℝ) + β + 2)) := by
    calc
      _ ≤ ((n.factorial : ℝ) * ((n : ℝ) + k) ^ k) * n.factorial := by
        apply mul_le_mul_of_nonneg_right _ hF
        simpa only [add_assoc] using hD
      _ = ((n : ℝ) + k) ^ k * ((n.factorial : ℝ) * n.factorial) := by ring
      _ ≤ _ := by
        apply mul_le_mul_of_nonneg_left _ (by positivity)
        apply mul_le_mul _ _ hF hAg.le
        · simpa only [add_assoc] using hA
        · simpa only [add_assoc] using hB
  rw [integral_sq_jacobi_succ hα hβ, inv_div, Nat.factorial_succ,
    Nat.cast_mul, Nat.cast_add, Nat.cast_one]
  simp only [GammaOn, GammaPoint.ofPos]
  apply (div_le_div_iff₀ (mul_pos (mul_pos hR hAg) hBg) hR).mpr
  convert mul_le_mul_of_nonneg_left hprod
    (mul_nonneg (mul_nonneg hd (show 0 ≤ (n : ℝ) + 1 by positivity)) hR.le) using 1 <;> ring

/-- An explicit fixed-power bound for reciprocal norms in degrees at least three. -/
theorem inv_integral_sq_jacobi_succ_le_polynomial {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (n k : ℕ) (hn : 2 ≤ n)
    (hk : α + β + 1 ≤ k) :
    (∫ x, (jacobi α β (n + 1)).eval x ^ 2 ∂jacobiMeasure α β)⁻¹ ≤
      ((2 + |α + β + 3|) * ((k : ℝ) + 1) ^ k / (2 : ℝ) ^ (α + β + 1)) *
        ((n : ℝ) + 1) ^ (k + 2) := by
  apply (inv_integral_sq_jacobi_succ_le hα hβ n k hn hk).trans
  have hn0 : (0 : ℝ) ≤ n := Nat.cast_nonneg _
  have hd : 0 ≤ 2 * (n : ℝ) + α + β + 3 := by linarith
  calc
    _ ≤ (((2 + |α + β + 3|) * ((n : ℝ) + 1)) * (n + 1) *
        (((k : ℝ) + 1) * (n + 1)) ^ k) / (2 : ℝ) ^ (α + β + 1) := by
      gcongr
      · nlinarith [le_abs_self (α + β + 3), mul_nonneg hn0 (abs_nonneg (α + β + 3))]
      · nlinarith [Nat.cast_nonneg (α := ℝ) k]
    _ = _ := by rw [mul_pow, pow_add, pow_two]; ring

/-- Reciprocal squared norms have a global polynomial bound, including degree zero. -/
theorem inv_integral_sq_jacobi_le_polynomial {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (k : ℕ) (hk : α + β + 1 ≤ k) (n : ℕ) :
    (∫ x, (jacobi α β n).eval x ^ 2 ∂jacobiMeasure α β)⁻¹ ≤
      jacobiReciprocalNormBound ⟨α, hα⟩ ⟨β, hβ⟩ k * ((n : ℝ) + 1) ^ (k + 2) := by
  let B : ℝ := (2 + |α + β + 3|) * ((k : ℝ) + 1) ^ k / (2 : ℝ) ^ (α + β + 1)
  have hB : 0 ≤ B := by dsimp [B]; positivity
  have h0 := inv_nonneg.mpr (integral_sq_jacobi_pos hα hβ 0).le
  have h1 := inv_nonneg.mpr (integral_sq_jacobi_pos hα hβ 1).le
  have h2 := inv_nonneg.mpr (integral_sq_jacobi_pos hα hβ 2).le
  have hC := (jacobiReciprocalNormBound_pos hα hβ k).le
  have hBC : B ≤ jacobiReciprocalNormBound ⟨α, hα⟩ ⟨β, hβ⟩ k := by
    change B ≤ 1 + _ + _ + _ + B
    linarith
  have hsmall (j : ℕ) (hj : j ≤ 2) :
      (∫ x, (jacobi α β j).eval x ^ 2 ∂jacobiMeasure α β)⁻¹ ≤
        jacobiReciprocalNormBound ⟨α, hα⟩ ⟨β, hβ⟩ k := by
    unfold jacobiReciprocalNormBound
    change _ ≤ 1 + _ + _ + _ + B
    interval_cases j <;> linarith
  by_cases hn : n ≤ 2
  · apply (hsmall n hn).trans
    exact le_mul_of_one_le_right hC (one_le_pow₀ (by linarith [Nat.cast_nonneg (α := ℝ) n]))
  obtain ⟨m, rfl⟩ := Nat.exists_eq_succ_of_ne_zero (show n ≠ 0 by omega)
  have hm : 2 ≤ m := by omega
  apply (inv_integral_sq_jacobi_succ_le_polynomial hα hβ m k hm hk).trans
  change B * ((m : ℝ) + 1) ^ (k + 2) ≤ _
  apply mul_le_mul hBC _ (by positivity) hC
  apply pow_le_pow_left₀ (by positivity)
  push_cast
  linarith

/-- Uniform polynomial growth of the diagonal Jacobi Christoffel kernel. -/
theorem jacobi_kernel_le_polynomial {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (k : ℕ) (hkα : α + 1 ≤ k) (hkβ : β ≤ k) (hks : α + β + 1 ≤ k)
    (N : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    (jacobiSequence hα hβ).kernel (jacobiMeasure α β) N x x (Polynomial.Sequence.KernelInput.of _ _ _ (hi := fun _ => (hasFiniteMoments_jacobiMeasure hα hβ).integrable_mul _ _)) ≤
      jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k * ((N : ℝ) + 1) ^ (3 * k + 3) := by
  simp only [Sequence.kernel_of_positive (hasPositiveDefiniteIntegral_jacobiMeasure hα hβ) (jacobiSequence hα hβ), jacobiSequence_apply, squaredNorm_jacobiSequence]
  calc
    _ ≤ ∑ n ∈ Finset.range (N + 1),
        jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k * ((N : ℝ) + 1) ^ (3 * k + 2) := by
      apply Finset.sum_le_sum
      intro n hn
      apply (jacobi_kernel_term_le_polynomial hα hβ k hkα hkβ hks n hx).trans
      apply mul_le_mul_of_nonneg_left _ (jacobiKernelBound_nonneg hα hβ k)
      apply pow_le_pow_left₀ (by positivity)
      have hn' : n ≤ N := Nat.lt_succ_iff.mp (Finset.mem_range.mp hn)
      exact_mod_cast Nat.add_le_add_right hn' 1
    _ = _ := by
      simp only [Finset.sum_const, Finset.card_range, nsmul_eq_mul, Nat.cast_add, Nat.cast_one]
      rw [show 3 * k + 3 = (3 * k + 2) + 1 by omega, pow_succ]
      ring

/-- A polynomial bound for the square root of the kernel, sufficient for
uniform geometric convergence after multiplication by the analytic error. -/
theorem sqrt_jacobi_kernel_le_polynomial {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (k : ℕ) (hkα : α + 1 ≤ k) (hkβ : β ≤ k) (hks : α + β + 1 ≤ k)
    (N : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    √((jacobiSequence hα hβ).kernel (jacobiMeasure α β) N x x (Polynomial.Sequence.KernelInput.of _ _ _ (hi := fun _ => (hasFiniteMoments_jacobiMeasure hα hβ).integrable_mul _ _))) ≤
      √(jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k) * ((N : ℝ) + 1) ^ (3 * k + 3) := by
  apply (Real.sqrt_le_sqrt (jacobi_kernel_le_polynomial hα hβ k hkα hkβ hks N hx)).trans
  rw [Real.sqrt_mul (jacobiKernelBound_nonneg hα hβ k)]
  apply mul_le_mul_of_nonneg_left _ (Real.sqrt_nonneg _)
  exact Real.sqrt_le_self_iff.mpr (Or.inr (one_le_pow₀ (by linarith [Nat.cast_nonneg (α := ℝ) N])))

/-- Scalar extension preserves the L² norm. -/
theorem norm_ofReal_compLp (f : Lp ℝ 2 μ) :
    ‖(@RCLike.ofRealCLM 𝕂 _).compLp f‖ = ‖f‖ := by
  have h := congrArg RCLike.re (inner_ofReal_compLp 𝕂 f f)
  rw [← norm_sq_eq_re_inner, RCLike.ofReal_re, real_inner_self_eq_norm_sq] at h
  nlinarith [norm_nonneg ((@RCLike.ofRealCLM 𝕂 _).compLp f), norm_nonneg f]

/-- The scalar-extended polynomial inner product with arbitrary L² data is
its integral pairing with the real polynomial. -/
theorem inner_polynomialToL2Scalar_left (p : ℝ[X]) (f : Lp 𝕂 2 μ) :
    inner 𝕂 (polynomialToL2Scalar 𝕂 hμ p) f = ∫ x, f x * ((p.eval x : ℝ) : 𝕂) ∂μ := by
  rw [L2.inner_def]
  apply integral_congr_ae
  filter_upwards [coe_polynomialToL2Scalar 𝕂 hμ p] with x hp
  simp [hp, RCLike.inner_apply]

/-- Real orthogonal polynomial normalization remains orthonormal after scalar extension. -/
theorem IsOrthogonal.orthonormal_normalizedToL2Scalar (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) : Orthonormal 𝕂 (fun n ↦ P.normalizedToL2Scalar 𝕂 hμ ⟨n, squaredNorm_pos hpos P n⟩) := by
  rw [orthonormal_iff_ite]
  intro n m
  rw [normalizedToL2Scalar, normalizedToL2Scalar, inner_ofReal_compLp,
    orthonormal_iff_ite.mp (hP.orthonormal_normalizedToL2 hμ hpos)]
  split_ifs <;> simp

/-- Scalar-extended polynomial classes of degree at most N belong to the first
N+1 normalized vectors over the enlarged scalar field. -/
theorem polynomialToL2Scalar_mem_span (hpos : HasPositiveDefiniteIntegral μ)
    (N : ℕ) (p : ℝ[X]) (hp : p.degree ≤ N) :
    polynomialToL2Scalar 𝕂 hμ p ∈ Submodule.span 𝕂
      ((fun n ↦ P.normalizedToL2Scalar 𝕂 hμ ⟨n, squaredNorm_pos hpos P n⟩) '' (Finset.range (N + 1) : Set ℕ)) := by
  have hs := polynomialToL2_mem_span_normalizedToL2 (P := P) hμ hpos N p hp
  clear hp
  change ((@RCLike.ofRealCLM 𝕂 _).compLpₗ 2 μ) (polynomialToL2 hμ p) ∈ _
  generalize polynomialToL2 hμ p = v at hs ⊢
  induction hs using Submodule.span_induction with
  | mem v hv =>
    obtain ⟨i, hi, rfl⟩ := hv
    exact Submodule.subset_span (Set.mem_image_of_mem _ hi)
  | zero =>
    rw [map_zero]
    exact Submodule.zero_mem _
  | add u v _ _ hu hv => simpa only [map_add] using Submodule.add_mem _ hu hv
  | smul c u _ hu =>
    rw [map_smul, RCLike.real_smul_eq_coe_smul (K := 𝕂)]
    exact Submodule.smul_mem _ _ hu

/-- Integral-coefficient terms represent the finite orthonormal projection. -/
theorem inner_normalizedToL2Scalar_smul (hpos : HasPositiveDefiniteIntegral μ)
    (f : Lp 𝕂 2 μ) (n : ℕ) :
    inner 𝕂 (P.normalizedToL2Scalar 𝕂 hμ ⟨n, squaredNorm_pos hpos P n⟩) f • P.normalizedToL2Scalar 𝕂 hμ ⟨n, squaredNorm_pos hpos P n⟩ =
      P.integralCoefficientScalar 𝕂 μ (IntegralCoefficientInput.ofL2 hμ hpos (P n) (P.ne_zero n) f) • polynomialToL2Scalar 𝕂 hμ (P n) := by
  rw [normalizedToL2Scalar_eq 𝕂 hμ P ⟨n, squaredNorm_pos hpos P n⟩, smul_smul]
  congr 1
  simpa only [normalizedToL2Scalar_eq] using
    (inner_normalizedToL2Scalar_mul_inv_sqrt 𝕂 hμ (P := P) hpos f n)

/-- The scalar-extended kernel vector reproduces partial-sum evaluation. -/
theorem inner_kernelVectorScalar (hpos : HasPositiveDefiniteIntegral μ)
    (N : ℕ) (x : ℝ) (f : Lp 𝕂 2 μ) :
    inner 𝕂 (P.kernelVectorScalar 𝕂 hμ ⟨N, fun i _ ↦ squaredNorm_pos hpos P i⟩ x) f = P.partialSumValueScalar 𝕂 μ (IntegralCoefficientFamilyInput.ofL2 hμ hpos P f N) x := by
  rw [kernelVectorScalar, kernelVector_of_positive hμ P hpos]
  change inner 𝕂 (((@RCLike.ofRealCLM 𝕂 _).compLpₗ 2 μ)
    (∑ i ∈ Finset.range (N + 1), ((P i).eval x / √(P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _))) •
      P.normalizedToL2 hμ ⟨i, squaredNorm_pos hpos P i⟩)) f = _
  simp only [map_sum, map_smul, sum_inner,
    RCLike.real_smul_eq_coe_smul (K := 𝕂), inner_smul_left, RCLike.conj_ofReal,
    partialSumValueScalar_ofL2]
  apply Finset.sum_congr rfl
  intro i _
  change (((P i).eval x / √(P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _)) : ℝ) : 𝕂) *
    inner 𝕂 (P.normalizedToL2Scalar 𝕂 hμ ⟨i, squaredNorm_pos hpos P i⟩) f = _
  rw [div_eq_mul_inv, RCLike.ofReal_mul]
  calc
    _ = (inner 𝕂 (P.normalizedToL2Scalar 𝕂 hμ ⟨i, squaredNorm_pos hpos P i⟩) f *
        (((√(P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _)))⁻¹ : ℝ) : 𝕂)) * (((P i).eval x : ℝ) : 𝕂) := by ring
    _ = _ := by rw [inner_normalizedToL2Scalar_mul_inv_sqrt 𝕂 hμ hpos]

/-- Scalar extension leaves the norm of the Christoffel kernel vector unchanged. -/
theorem IsOrthogonal.norm_kernelVectorScalar (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) (N : ℕ) (x : ℝ) :
    ‖P.kernelVectorScalar 𝕂 hμ ⟨N, fun i _ ↦ squaredNorm_pos hpos P i⟩ x‖ = √(P.kernel μ N x x (Polynomial.Sequence.KernelInput.of _ _ _)) := by
  rw [kernelVectorScalar, norm_ofReal_compLp, hP.norm_kernelVector hμ hpos]

/-- The scalar kernel reproduces every real polynomial of degree at most N. -/
theorem IsOrthogonal.inner_kernelVectorScalar_polynomialToL2Scalar (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) (N : ℕ) (x : ℝ)
    (p : ℝ[X]) (hp : p.degree ≤ N) :
    inner 𝕂 (P.kernelVectorScalar 𝕂 hμ ⟨N, fun i _ ↦ squaredNorm_pos hpos P i⟩ x) (polynomialToL2Scalar 𝕂 hμ p) =
      ((p.eval x : ℝ) : 𝕂) := by
  change inner 𝕂 ((@RCLike.ofRealCLM 𝕂 _).compLp _) ((@RCLike.ofRealCLM 𝕂 _).compLp _) = _
  rw [inner_ofReal_compLp, inner_kernelVector hμ hpos,
    hP.partialSum_polynomialToL2 hμ hpos N p hp]

/-- The general scalar coefficient specializes on the full classical input domain. -/
theorem integralCoefficientScalar_jacobiSequence {n : ℕ}
    (f : IntegralCoefficientInput ℂ (jacobiMeasure α β) (jacobi α β n)) :
    (jacobiSequence hα hβ).integralCoefficientScalar ℂ (jacobiMeasure α β) f =
      jacobiComplexCoefficient α β f := by
  simp only [Sequence.integralCoefficientScalar, squaredNorm_jacobiSequence,
    jacobiSequence_apply, jacobiComplexCoefficient]
  rfl

/-- The general scalar coefficient specializes on the full classical input domain. -/
theorem integralCoefficientScalar_gegenbauerSequence {n : ℕ}
    (f : IntegralCoefficientInput ℂ (gegenbauerMeasure lam) (gegenbauer lam n)) :
    (gegenbauerSequence hlam hzero).integralCoefficientScalar ℂ (gegenbauerMeasure lam) f =
      gegenbauerComplexCoefficient lam f := by
  simp only [Sequence.integralCoefficientScalar, squaredNorm_gegenbauerSequence,
    gegenbauerSequence_apply, gegenbauerComplexCoefficient]
  rfl

/-- Modes below the derivative order vanish, including degenerate parameters. -/
theorem iterate_derivative_jacobi_eq_zero (α β : R) {k n : ℕ} (hn : n < k) :
    derivative^[k] (jacobi α β n) = 0 := by
  exact iterate_derivative_eq_zero ((natDegree_jacobi_le α β n).trans_lt hn)

/-- Uniform complex Jacobi approximation, with an explicit nonsharp polynomial factor. -/
theorem norm_sub_sum_jacobiComplexCoefficient_le_uniform_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
    (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
    (hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
    (f : Lp ℂ 2 (jacobiMeasure α β)) (hf : ∀ᵐ x ∂jacobiMeasure α β, f x = F x)
    (k : ℕ) (hkα : α + 1 ≤ k) (hkβ : β ≤ k) (hks : α + β + 1 ≤ k)
    (N : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    ‖F x - ∑ i ∈ Finset.range (N + 1),
        (jacobiComplexCoefficient α β (IntegralCoefficientInput.ofJacobiL2 hα hβ f i)) * (((jacobi α β i).eval x : ℝ) : ℂ)‖ ≤
      (1 + √(jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k) * ((N : ℝ) + 1) ^ (3 * k + 3) *
        √((jacobiMeasure α β).real univ)) * (2 * M / (ρ ^ N * (ρ - 1))) := by
  have h := (isOrthogonal_jacobiSequence hα hβ).norm_sub_partialSumValueScalar_le_kernel_of_analytic
    (hasFiniteMoments_jacobiMeasure hα hβ) (hasPositiveDefiniteIntegral_jacobiMeasure hα hβ)
    ((ae_mem_Ioo_jacobiMeasure α β).mono (fun _ hx ↦ Ioo_subset_Icc_self hx)) hρ hF hM f hf N hx
  simp only [Sequence.partialSumValueScalar_ofL2, Sequence.integralCoefficientScalar, squaredNorm_jacobiSequence,
    jacobiSequence_apply, jacobiComplexCoefficient, IntegralCoefficientInput.ofJacobiL2,
    IntegralCoefficientInput.ofL2, RCLike.ofReal_eq_complex_ofReal] at h ⊢
  have hM0 : 0 ≤ M := (norm_nonneg (F x)).trans (hM x (ofReal_mem_bernsteinEllipse hρ hx))
  apply h.trans
  gcongr
  exact sqrt_jacobi_kernel_le_polynomial hα hβ k hkα hkβ hks N hx

/-- Complex analytic Jacobi expansions converge uniformly on the closed interval. -/
theorem tendstoUniformlyOn_sum_jacobiComplexCoefficient_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
    (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
    (hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
    (f : Lp ℂ 2 (jacobiMeasure α β)) (hf : ∀ᵐ x ∂jacobiMeasure α β, f x = F x) :
    TendstoUniformlyOn (fun N x ↦ ∑ i ∈ Finset.range (N + 1),
      (jacobiComplexCoefficient α β (IntegralCoefficientInput.ofJacobiL2 hα hβ f i)) * (((jacobi α β i).eval x : ℝ) : ℂ))
      (fun x : ℝ ↦ F x) atTop (Icc (-1) 1) := by
  obtain ⟨k, hkα, hkβ, hks⟩ := exists_jacobiKernel_exponent α β
  apply tendstoUniformlyOn_of_polynomial_geometric_bound (3 * k + 3) hρ
    (A := √(jacobiKernelBound ⟨α, hα⟩ ⟨β, hβ⟩ k) * √((jacobiMeasure α β).real univ))
    (B := ρ - 1) (C := 2 * M)
  intro N x hx
  convert norm_sub_sum_jacobiComplexCoefficient_le_uniform_of_analytic hα hβ hρ hF hM f hf
    k hkα hkβ hks N hx using 1
  ring

/-- Geometric weighted L² error for complex analytic Gegenbauer expansions. -/
theorem norm_sub_sum_gegenbauerComplexCoefficient_le_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
    (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
    (hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
    (f : Lp ℂ 2 (gegenbauerMeasure lam)) (hf : ∀ᵐ x ∂gegenbauerMeasure lam, f x = F x) (N : ℕ) :
    ‖f - ∑ i ∈ Finset.range (N + 1), (gegenbauerComplexCoefficient lam (IntegralCoefficientInput.ofGegenbauerL2 hlam hzero f i)) •
      polynomialToL2Scalar ℂ (hasFiniteMoments_gegenbauerMeasure hlam) (gegenbauer lam i)‖ ≤
      √((gegenbauerMeasure lam).real univ) * (2 * M / (ρ ^ N * (ρ - 1))) := by
  have hP := isOrthogonal_gegenbauerSequence hlam hzero
  simpa only [Sequence.integralCoefficientScalar, squaredNorm_gegenbauerSequence, gegenbauerSequence_apply,
    gegenbauerComplexCoefficient, IntegralCoefficientInput.ofGegenbauerL2, IntegralCoefficientInput.ofL2, RCLike.ofReal_eq_complex_ofReal] using
    hP.norm_sub_sum_integralCoefficientScalar_le_of_analytic
      (hasFiniteMoments_gegenbauerMeasure hlam) (hasPositiveDefiniteIntegral_gegenbauerMeasure hlam)
      ((ae_mem_Ioo_gegenbauerMeasure lam).mono (fun _ hx ↦ Ioo_subset_Icc_self hx)) hρ hF hM f hf N

/-- Uniform complex Gegenbauer approximation throughout the nondegenerate range. -/
theorem norm_sub_sum_gegenbauerComplexCoefficient_le_uniform_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
    (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
    (hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
    (f : Lp ℂ 2 (gegenbauerMeasure lam)) (hf : ∀ᵐ x ∂gegenbauerMeasure lam, f x = F x)
    (k : ℕ) (hk : lam + 1 / 2 ≤ k) (hks : 2 * lam ≤ k)
    (N : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    ‖F x - ∑ i ∈ Finset.range (N + 1),
        (gegenbauerComplexCoefficient lam (IntegralCoefficientInput.ofGegenbauerL2 hlam hzero f i)) * (((gegenbauer lam i).eval x : ℝ) : ℂ)‖ ≤
      (1 + √(jacobiKernelBound ⟨lam - 1 / 2, by change -1 < lam - 1 / 2; linarith⟩ ⟨lam - 1 / 2, by change -1 < lam - 1 / 2; linarith⟩ k) *
        ((N : ℝ) + 1) ^ (3 * k + 3) * √((gegenbauerMeasure lam).real univ)) *
        (2 * M / (ρ ^ N * (ρ - 1))) := by
  have hP := isOrthogonal_gegenbauerSequence hlam hzero
  have h := hP.norm_sub_partialSumValueScalar_le_kernel_of_analytic
      (hasFiniteMoments_gegenbauerMeasure hlam) (hasPositiveDefiniteIntegral_gegenbauerMeasure hlam)
      ((ae_mem_Ioo_gegenbauerMeasure lam).mono (fun _ hx ↦ Ioo_subset_Icc_self hx))
      hρ hF hM f hf N hx
  simp only [Sequence.partialSumValueScalar_ofL2, Sequence.integralCoefficientScalar, squaredNorm_gegenbauerSequence,
    gegenbauerSequence_apply, gegenbauerComplexCoefficient, IntegralCoefficientInput.ofGegenbauerL2,
    IntegralCoefficientInput.ofL2, RCLike.ofReal_eq_complex_ofReal, gegenbauer_kernel_eq_jacobi hlam hzero] at h ⊢
  have hM0 : 0 ≤ M := (norm_nonneg (F x)).trans (hM x (ofReal_mem_bernsteinEllipse hρ hx))
  apply h.trans
  gcongr
  exact sqrt_jacobi_kernel_le_polynomial (by linarith) (by linarith) k
    (by linarith) (by linarith) (by linarith) N hx

/-- Complex analytic Gegenbauer expansions converge uniformly on the closed interval. -/
theorem tendstoUniformlyOn_sum_gegenbauerComplexCoefficient_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
    (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
    (hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
    (f : Lp ℂ 2 (gegenbauerMeasure lam)) (hf : ∀ᵐ x ∂gegenbauerMeasure lam, f x = F x) :
    TendstoUniformlyOn (fun N x ↦ ∑ i ∈ Finset.range (N + 1),
      (gegenbauerComplexCoefficient lam (IntegralCoefficientInput.ofGegenbauerL2 hlam hzero f i)) * (((gegenbauer lam i).eval x : ℝ) : ℂ))
      (fun x : ℝ ↦ F x) atTop (Icc (-1) 1) := by
  obtain ⟨k, hk0, hks⟩ := exists_gegenbauerKernel_exponent lam
  apply tendstoUniformlyOn_of_polynomial_geometric_bound (3 * k + 3) hρ
    (A := √(jacobiKernelBound ⟨lam - 1 / 2, by change -1 < lam - 1 / 2; linarith⟩ ⟨lam - 1 / 2, by change -1 < lam - 1 / 2; linarith⟩ k) *
      √((gegenbauerMeasure lam).real univ)) (B := ρ - 1) (C := 2 * M)
  intro N x hx
  convert norm_sub_sum_gegenbauerComplexCoefficient_le_uniform_of_analytic hlam hzero hρ hF hM f hf
    k hk0 hks N hx using 1
  ring