/-- Connection denominators commute with ring homomorphisms. -/
@[simp] theorem map_jacobiAlphaConnectionDenominator (f : R →+* S) (α β : R) (n k : ℕ) :
    f (jacobiAlphaConnectionDenominator α β n k) =
      jacobiAlphaConnectionDenominator (f α) (f β) n k := by
  have hp (j : ℕ) (x : R) : f ((ascPochhammer R j).eval x) =
      (ascPochhammer S j).eval (f x) := by
    rw [← ascPochhammer_map f, eval_map]
    exact (eval₂_at_apply f x).symm
  simp only [jacobiAlphaConnectionDenominator, map_mul, hp, map_add, map_natCast,
    map_one, map_ofNat]

/-- The common connection denominator commutes with ring homomorphisms. -/
@[simp] theorem map_jacobiAlphaConnectionScale (f : R →+* S) (α β : R) (n : ℕ) :
    f (jacobiAlphaConnectionScale α β n) = jacobiAlphaConnectionScale (f α) (f β) n := by
  simp [jacobiAlphaConnectionScale]

/-- Cleared connection coefficients commute with ring homomorphisms. -/
@[simp] theorem map_jacobiAlphaConnectionCleared (f : R →+* S) (α β γ : R) (n : ℕ)
    (k : Fin (n + 1)) :
    f (jacobiAlphaConnectionCleared α β γ n k) =
      jacobiAlphaConnectionCleared (f α) (f β) (f γ) n k := by
  simp [jacobiAlphaConnectionCleared]

/-- The explicit one-parameter Jacobi connection formula with denominators
cleared, valid even at degenerate parameters in every commutative rational algebra
(DLMF 18.18.14). -/
theorem jacobiAlphaConnectionScale_smul {R : Type*} [CommRing R] [Algebra ℚ R]
    (α β γ : R) (n : ℕ) :
    jacobiAlphaConnectionScale α β n • jacobi γ β n =
      ∑ k : Fin (n + 1), jacobiAlphaConnectionCleared α β γ n k • jacobi α β k := by
  let a : MvPolynomial (Fin 3) ℚ := MvPolynomial.X 0
  let b : MvPolynomial (Fin 3) ℚ := MvPolynomial.X 1
  let c : MvPolynomial (Fin 3) ℚ := MvPolynomial.X 2
  have hid : jacobiAlphaConnectionScale a b n • jacobi c b n =
      ∑ k : Fin (n + 1), jacobiAlphaConnectionCleared a b c n k • jacobi a b k := by
    apply eq_of_map_eval₂_eq_on (fun _ ↦ Set.Ioi (-1)) (fun _ ↦ Set.Ioi_infinite _)
    intro v hv
    change (jacobiAlphaConnectionScale a b n • jacobi c b n).map
      (MvPolynomial.aeval v).toRingHom =
      (∑ k : Fin (n + 1), jacobiAlphaConnectionCleared a b c n k • jacobi a b k).map
        (MvPolynomial.aeval v).toRingHom
    simp only [Polynomial.map_smul, Polynomial.map_sum, map_jacobi]
    simpa [a, b, c] using
      jacobiAlphaConnectionScale_smul_of_neg_one_lt (hv 0 (Set.mem_univ _))
        (hv 1 (Set.mem_univ _)) (v 2) n
  let f : MvPolynomial (Fin 3) ℚ →ₐ[ℚ] R := MvPolynomial.aeval ![α, β, γ]
  have h := congrArg (Polynomial.map f.toRingHom) hid
  simp only [Polynomial.map_smul, Polynomial.map_sum, map_jacobi] at h
  simpa [f, a, b, c] using h

/-- The second-parameter connection formula follows by reflection from
DLMF 18.18.14, with all denominators cleared. -/
theorem jacobiBetaConnectionScale_smul (α β δ : R) (n : ℕ) :
    jacobiAlphaConnectionScale β α n • jacobi α δ n =
      ∑ k : Fin (n + 1),
        ((-1 : R) ^ (n + k) * jacobiAlphaConnectionCleared β α δ n k) • jacobi α β k := by
  have h := congrArg (fun p : R[X] ↦ p.comp (-X)) (jacobiAlphaConnectionScale_smul β α δ n)
  simp only [smul_comp, sum_comp, jacobi_comp_neg_X, ← smul_eq_C_mul, smul_smul] at h
  have hh := congrArg (fun p : R[X] ↦ (-1 : R) ^ n • p) h
  rw [Finset.smul_sum] at hh
  simp only [smul_smul] at hh
  have he : (-1 : R) ^ n * (-1 : R) ^ n = 1 := by rw [← mul_pow]; simp
  rw [show (-1 : R) ^ n * (jacobiAlphaConnectionScale β α n * (-1 : R) ^ n) =
    jacobiAlphaConnectionScale β α n by rw [mul_left_comm, he, mul_one]] at hh
  convert hh using 1
  apply Finset.sum_congr rfl
  intro k _
  congr 1
  rw [pow_add]
  ring

/-- The classical coefficient recovers its determinant numerator. -/
theorem jacobiFiniteConnectionCoefficient_mul_det (α β γ δ : K) (N : ℕ)
    (j : Fin (N+1)) (d : FiniteConnectionInput (fun i : Fin (N+1) => jacobi α β i)) :
    jacobiFiniteConnectionCoefficient α β γ δ N j d * (jacobiConnectionMatrix α β N).det =
      jacobiConnectionNumerator α β γ δ N j :=
  div_mul_cancel₀ _ d.nonsingular

/-- The determinant of a degree-indexed family's coefficient matrix is the
product of its top coefficients; no nonvanishing assumption is needed. -/
theorem det_connectionMatrix (p : Fin N → R[X]) (hp : ∀ j, (p j).natDegree ≤ j) :
    (connectionMatrix p).det = ∏ j : Fin N, (p j).coeff j :=
  Matrix.det_of_isUpperTriangular (Matrix.matrixOfPolynomials_isUpperTriangular p hp)

/-- The universal cleared finite connection formula (Cramer's rule).
Both the family members and the polynomial have degree less than the matrix size. -/
theorem det_smul_eq_sum_connectionNumerator (p : Fin N → R[X]) (q : R[X])
    (hp : ∀ j, (p j).degree < N) (hq : q.degree < N) :
    (connectionMatrix p).det • q = ∑ j : Fin N, connectionNumerator p q j • p j := by
  ext k
  simp only [coeff_smul, finsetSum_coeff, smul_eq_mul]
  by_cases hk : k < N
  · have h := congrFun (Matrix.mulVec_cramer (connectionMatrix p) (fun i ↦ q.coeff i))
      (⟨k, hk⟩ : Fin N)
    simpa only [Matrix.mulVec, dotProduct, connectionMatrix, Matrix.of_apply,
      Pi.smul_apply, smul_eq_mul, connectionNumerator, mul_comm] using h.symm
  · have hqk : q.coeff k = 0 := coeff_eq_zero_of_degree_lt
      (hq.trans_le (Nat.cast_le.mpr (Nat.le_of_not_gt hk)))
    rw [hqk, mul_zero]
    symm
    apply Finset.sum_eq_zero
    intro j _
    rw [coeff_eq_zero_of_degree_lt ((hp j).trans_le (Nat.cast_le.mpr (Nat.le_of_not_gt hk))),
      mul_zero]

/-- Dividing the universal identity yields the finite polynomial expansion. -/
theorem eq_sum_finiteConnectionCoefficient (p : Fin N → K[X]) (q : K[X])
    (hp : ∀ j, (p j).degree < N) (hq : q.degree < N)
    (hd : (connectionMatrix p).det ≠ 0) :
    q = ∑ j : Fin N, finiteConnectionCoefficient p q j ⟨hd⟩ • p j := by
  have h := congrArg (fun r : K[X] ↦ (connectionMatrix p).det⁻¹ • r)
    (det_smul_eq_sum_connectionNumerator p q hp hq)
  simpa only [smul_smul, inv_mul_cancel₀ hd, one_smul, Finset.smul_sum,
    finiteConnectionCoefficient, div_eq_mul_inv, mul_comm] using h

/-- A universal finite connection identity between arbitrary Jacobi parameters,
with the determinant denominator cleared, including singular specialisations. -/
theorem det_smul_jacobi_eq_sum_connectionNumerator (α β γ δ : R) (N : ℕ) :
    (jacobiConnectionMatrix α β N).det • jacobi γ δ N =
      ∑ j : Fin (N + 1), jacobiConnectionNumerator α β γ δ N j • jacobi α β j := by
  apply det_smul_eq_sum_connectionNumerator
  · intro j
    exact (degree_jacobi_le α β j).trans_lt (Nat.cast_lt.mpr j.isLt)
  · exact (degree_jacobi_le γ δ N).trans_lt (Nat.cast_lt.mpr (Nat.lt_succ_self N))

/-- The finite Jacobi expansion when the target determinant is nonzero.
No nondegeneracy condition is imposed on the source parameters. -/
theorem jacobi_eq_sum_finiteConnectionCoefficient (α β γ δ : K) (N : ℕ)
    (hd : (jacobiConnectionMatrix α β N).det ≠ 0) :
    jacobi γ δ N =
      ∑ j : Fin (N + 1), jacobiFiniteConnectionCoefficient α β γ δ N j ⟨hd⟩ • jacobi α β j := by
  exact eq_sum_finiteConnectionCoefficient _ _
    (fun j ↦ (degree_jacobi_le α β j).trans_lt (Nat.cast_lt.mpr j.isLt))
    ((degree_jacobi_le γ δ N).trans_lt (Nat.cast_lt.mpr (Nat.lt_succ_self N))) hd

/-- Positivity of target leading coefficients makes the finite determinant positive. -/
theorem det_jacobiConnectionMatrix_pos {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (N : ℕ) :
    0 < (jacobiConnectionMatrix α β N).det := by
  rw [jacobiConnectionMatrix, det_connectionMatrix _ (fun j ↦ natDegree_jacobi_le α β j)]
  exact Finset.prod_pos (fun j _ ↦ coeff_jacobi_self_pos hα hβ j)

/-- Arbitrary real Jacobi parameters expand into every target in the full
orthogonality range, with the explicit finite determinant coefficients. -/
theorem jacobi_eq_sum_finiteConnectionCoefficient_of_neg_one_lt {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (γ δ : ℝ) (N : ℕ) :
    jacobi γ δ N =
      ∑ j : Fin (N + 1), jacobiFiniteConnectionCoefficient α β γ δ N j ⟨(det_jacobiConnectionMatrix_pos hα hβ N).ne'⟩ • jacobi α β j :=
  jacobi_eq_sum_finiteConnectionCoefficient α β γ δ N (det_jacobiConnectionMatrix_pos hα hβ N).ne'

/-- Determinant coefficients agree with the existing polynomial-basis coordinates. -/
theorem jacobiFiniteConnectionCoefficient_eq_repr {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (γ δ : ℝ) (N : ℕ) (j : Fin (N + 1)) :
    jacobiFiniteConnectionCoefficient α β γ δ N j ⟨(det_jacobiConnectionMatrix_pos hα hβ N).ne'⟩ =
      (jacobiBasis hα hβ).repr (jacobi γ δ N) j := by
  have he := congrArg (fun p : ℝ[X] ↦ (jacobiBasis hα hβ).repr p j)
    (jacobi_eq_sum_finiteConnectionCoefficient_of_neg_one_lt hα hβ γ δ N)
  have hr (i : Fin (N + 1)) :
      (jacobiBasis hα hβ).repr (jacobi α β i) = Finsupp.single (i : ℕ) 1 := by
    rw [← jacobiBasis_apply hα hβ, Module.Basis.repr_self]
  simp only [map_sum, map_smul, hr, Finsupp.finsetSum_apply,
    Finsupp.smul_apply, smul_eq_mul, Finsupp.single_apply] at he
  simpa only [← Fin.ext_iff, mul_ite, mul_one, mul_zero, Finset.sum_ite_eq',
    Finset.mem_univ, ite_true] using he.symm

/-- A one-step mixed-moment recurrence, with unrestricted source parameter `γ`
(DLMF 18.18.14). -/
theorem integral_jacobi_mixed_succ {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (γ : ℝ) (n : ℕ) :
    2 * (n + 1 : ℝ) * (∫ x, (jacobi γ β (n + 1)).eval x ∂jacobiMeasure α β) =
      (γ - α) * ∫ x, (jacobi (γ + 1) (β + 1) n).eval x ∂jacobiMeasure α (β + 1) := by
  have hi0 := integrable_jacobiMeasure hα hβ (jacobi γ β (n + 1)).continuous.continuousOn
  have hi1 := integrable_jacobiMeasure hα (by linarith : -1 < β + 1)
    (jacobi (γ + 1) (β + 1) n).continuous.continuousOn
  have hw0 := (integrable_jacobiMeasure_iff _ _ _).mp hi0
  have hw1 := (integrable_jacobiMeasure_iff _ _ _).mp hi1
  have hi := integral_eq_zero_of_hasDerivAt_jacobiWeight_mul hα hβ
    (jacobi (γ + 1) (β + 1) n).continuous.continuousOn
    (hasDerivAt_jacobiWeight_mul_jacobi_mixed_onOpen α β γ n)
    ((hw1.const_mul (γ - α)).sub (hw0.const_mul (2 * (n + 1 : ℝ))))
  rw [integral_sub (hw1.const_mul _) (hw0.const_mul _), integral_const_mul,
    integral_const_mul] at hi
  rw [integral_jacobiMeasure α β ⟨fun x ↦ (jacobi γ β (n + 1)).eval (x : ℝ), hi0⟩,
    integral_jacobiMeasure α (β + 1) ⟨fun x ↦ (jacobi (γ + 1) (β + 1) n).eval (x : ℝ), hi1⟩]
  exact (sub_eq_zero.mp hi).symm

/-- Integrating a Jacobi polynomial against a weight with a different first
parameter gives a rising factorial times a shifted mass (DLMF 18.18.14). -/
theorem integral_jacobi_mixed {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (γ : ℝ) (n : ℕ) :
    (2 : ℝ) ^ n * n.factorial * (∫ x, (jacobi γ β n).eval x ∂jacobiMeasure α β) =
      (ascPochhammer ℝ n).eval (γ - α) * ∫ _ : ℝ, (1 : ℝ) ∂jacobiMeasure α (β + n) := by
  induction n generalizing β γ with
  | zero => simp
  | succ n ih =>
    have hm := integral_jacobi_mixed_succ hα hβ γ n
    have hi := ih (by linarith : -1 < β + 1) (γ + 1)
    rw [ascPochhammer_succ_left]
    simp only [eval_mul, eval_X, eval_comp, eval_add, eval_one]
    rw [show γ + 1 - α = γ - α + 1 by ring,
      show β + 1 + (n : ℝ) = β + (n + 1 : ℕ) by push_cast; ring] at hi
    rw [pow_succ, Nat.factorial_succ, Nat.cast_mul, Nat.cast_add, Nat.cast_one]
    push_cast at hi
    linear_combination (2 : ℝ) ^ n * n.factorial * hm + (γ - α) * hi

/-- Integration by parts against a positive-degree Jacobi polynomial
(DLMF 18.18.14). -/
theorem integral_mul_jacobi_succ {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (q : ℝ[X]) (k : ℕ) :
    2 * (k + 1 : ℝ) * (∫ x, q.eval x * (jacobi α β (k + 1)).eval x ∂jacobiMeasure α β) =
      ∫ x, q.derivative.eval x * (jacobi (α + 1) (β + 1) k).eval x
        ∂jacobiMeasure (α + 1) (β + 1) := by
  have hi0 := integrable_jacobiMeasure hα hβ
    (q.continuous.mul (jacobi α β (k + 1)).continuous).continuousOn
  have hi1 := integrable_jacobiMeasure (by linarith : -1 < α + 1)
    (by linarith : -1 < β + 1)
    (q.derivative.continuous.mul (jacobi (α + 1) (β + 1) k).continuous).continuousOn
  have hw0 := (integrable_jacobiMeasure_iff _ _ _).mp hi0
  have hw1 := (integrable_jacobiMeasure_iff _ _ _).mp hi1
  simp only [Pi.mul_apply] at hw0 hw1
  have hi := integral_eq_zero_of_hasDerivAt_jacobiWeight_mul hα hβ
    (q.continuous.mul (jacobi (α + 1) (β + 1) k).continuous).continuousOn
    (hasDerivAt_jacobiWeight_mul_polynomial_mul_jacobi α β q k)
    (hw1.sub (hw0.const_mul (2 * (k + 1 : ℝ))))
  rw [integral_sub hw1 (hw0.const_mul _), integral_const_mul] at hi
  rw [integral_jacobiMeasure α β
    ⟨fun x ↦ q.eval (x : ℝ) * (jacobi α β (k + 1)).eval (x : ℝ), hi0⟩,
    integral_jacobiMeasure (α + 1) (β + 1)
      ⟨fun x ↦ q.derivative.eval (x : ℝ) * (jacobi (α + 1) (β + 1) k).eval (x : ℝ), hi1⟩]
  exact (sub_eq_zero.mp hi).symm

/-- Repeated integration by parts lowers the source and target degrees together.
The source first parameter is unrestricted (DLMF 18.18.14). -/
theorem integral_jacobi_mul_jacobi_shift {α β : ℝ} (hα : -1 < α) (hβ : -1 < β)
    (γ : ℝ) (k m : ℕ) :
    (4 : ℝ) ^ k * k.factorial *
      (∫ x, (jacobi γ β (k + m)).eval x * (jacobi α β k).eval x ∂jacobiMeasure α β) =
      (ascPochhammer ℝ k).eval ((k + m : ℕ) + γ + β + 1) *
        ∫ x, (jacobi (γ + k) (β + k) m).eval x ∂jacobiMeasure (α + k) (β + k) := by
  induction k generalizing α β γ with
  | zero => simp
  | succ k ih =>
    have hb := integral_mul_jacobi_succ hα hβ (jacobi γ β (k + 1 + m)) k
    have hd (x : ℝ) : (jacobi γ β (k + 1 + m)).derivative.eval (x : ℝ) =
        (((k + m : ℕ) + γ + β + 2) / 2) * (jacobi (γ + 1) (β + 1) (k + m)).eval (x : ℝ) := by
      rw [show k + 1 + m = (k + m) + 1 by omega, derivative_jacobi_succ, eval_mul, eval_C]
      norm_num [div_eq_mul_inv, mul_assoc, mul_comm, mul_left_comm]
    simp_rw [hd, mul_assoc] at hb
    rw [integral_const_mul] at hb
    have hi := ih (by linarith : -1 < α + 1) (by linarith : -1 < β + 1) (γ + 1)
    rw [ascPochhammer_succ_left]
    simp only [eval_mul, eval_X, eval_comp, eval_add, eval_one]
    have ha : α + 1 + (k : ℝ) = α + (k + 1 : ℕ) := by push_cast; ring
    have hβ' : β + 1 + (k : ℝ) = β + (k + 1 : ℕ) := by push_cast; ring
    have hγ : γ + 1 + (k : ℝ) = γ + (k + 1 : ℕ) := by push_cast; ring
    rw [ha, hβ', hγ] at hi
    rw [show (k + m : ℕ) + (γ + 1) + (β + 1) + 1 =
      (k + 1 + m : ℕ) + γ + β + 1 + 1 by push_cast; ring] at hi
    rw [pow_succ, Nat.factorial_succ, Nat.cast_mul, Nat.cast_add, Nat.cast_one]
    push_cast at hi hb ⊢
    linear_combination 2 * (4 : ℝ) ^ k * k.factorial * hb +
      ((k : ℝ) + m + γ + β + 2) * hi

/-- Raising the second weight parameter shifts the mass by a Pochhammer ratio,
with the denominator cleared (DLMF 18.18.14). -/
theorem integral_one_jacobiMeasure_beta_add_nat {α β : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (m : ℕ) :
    (ascPochhammer ℝ m).eval (α + β + 2) *
        (∫ _ : ℝ, (1 : ℝ) ∂jacobiMeasure α (β + m)) =
      (2 : ℝ) ^ m * (ascPochhammer ℝ m).eval (β + 1) *
        ∫ _ : ℝ, (1 : ℝ) ∂jacobiMeasure α β := by
  have hβm : -1 < β + (m : ℝ) := by linarith [Nat.cast_nonneg (α := ℝ) m]
  have hβΓ := Gamma_add_nat_eq_ascPochhammer_mul (by linarith : 0 < β + 1) m
  have hαβΓ := Gamma_add_nat_eq_ascPochhammer_mul (by linarith : 0 < α + β + 2) m
  simp only [GammaOn, GammaPoint.ofPos] at hβΓ hαβΓ
  rw [integral_one_jacobiMeasure hα hβm, integral_one_jacobiMeasure hα hβ]
  simp only [GammaOn, GammaPoint.ofPos]
  rw [show α + (β + (m : ℝ)) + 1 = (α + β + 1) + m by ring,
    rpow_add (by norm_num : (0 : ℝ) < 2), rpow_natCast,
    show β + (m : ℝ) + 1 = (β + 1) + m by ring,
    hβΓ,
    show α + (β + (m : ℝ)) + 2 = (α + β + 2) + m by ring,
    hαβΓ]
  have hp : (ascPochhammer ℝ m).eval (α + β + 2) ≠ 0 :=
    (ascPochhammer_pos m (α + β + 2) (by linarith)).ne'
  field_simp