/-- Integer parameter raising is multiplication of classical domain-typed weights. -/
theorem semiclassicalJacobiMeasure_add_nat (t α β γ : ℝ) (ht : 1 < t)
    (k l m : ℕ) :
    semiclassicalJacobiMeasure ⟨t, by assumption⟩ (α + k) (β + l) (γ + m) =
      openIntervalWeightMeasure 0 1
        (semiclassicalJacobiWeight ⟨t, ht⟩ α β γ *
          semiclassicalJacobiWeight ⟨t, ht⟩ k l m) := by
  change openIntervalWeightMeasure 0 1 _ = _
  congr 1
  ext x
  rw [ContinuousMap.mul_apply, NNReal.coe_mul, semiclassicalJacobiWeight_add_nat]
  simp only [semiclassicalJacobiWeight_apply, Real.rpow_natCast]
  ring

/-- The simultaneous unit shift on the classical domain. -/
theorem semiclassicalJacobiWeight_add_one (t : Ioi (1 : ℝ)) (α β γ : ℝ)
    (x : Ioo (0 : ℝ) 1) :
    (semiclassicalJacobiWeight t (α + 1) (β + 1) (γ + 1) x : ℝ) =
      (x : ℝ) * (1 - (x : ℝ)) * ((t : ℝ) - x) *
        (semiclassicalJacobiWeight t α β γ x : ℝ) := by
  simp only [semiclassicalJacobiWeight_apply]
  exact semiclassicalJacobiWeight_add_one_rpow t α β γ t.property x.property

/-- The connection matrix between two positive orthonormal semiclassical bases. -/
noncomputable def semiclassicalJacobiConnection {t α β α' β' : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα' : -1 < α') (hβ' : -1 < β') (γ' : ℝ) : Matrix ℕ ℕ ℝ :=
  (semiclassicalJacobiBasis ht hα' hβ' γ').toMatrix (semiclassicalJacobiBasis ht hα hβ γ)

/-- Every leading connection section is upper triangular in degree order. -/
theorem semiclassicalJacobiConnection_upperTriangular {t α β α' β' : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα' : -1 < α') (hβ' : -1 < β') (γ' : ℝ) (N : ℕ) :
    (Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiConnection ht hα hβ γ hα' hβ' γ' i j)).IsUpperTriangular := by
  simpa only [semiclassicalJacobiConnection, Basis.toMatrix_apply, semiclassicalJacobiBasis_apply,
    semiclassicalJacobiBasis, Sequence.basis_eq_self] using
    (semiclassicalJacobiSequence_spec ht hα' hβ' γ').1.connection_upperTriangular
      (hasFiniteMoments_semiclassicalJacobiMeasure ht hα' hβ')
      (hasPositiveDefiniteIntegral_semiclassicalJacobiMeasure ht hα' hβ')
      (fun n ↦ isUnit_iff_ne_zero.mpr ((semiclassicalJacobiSequence_spec ht hα' hβ' γ').2.1 n).ne')
      (semiclassicalJacobiSequence ht hα hβ γ) N

/-- Positive normalization makes each connection diagonal entry positive. -/
theorem semiclassicalJacobiConnection_diag_pos {t α β α' β' : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα' : -1 < α') (hβ' : -1 < β') (γ' : ℝ) (n : ℕ) :
    0 < semiclassicalJacobiConnection ht hα hβ γ hα' hβ' γ' n n := by
  simpa only [semiclassicalJacobiConnection, Basis.toMatrix_apply, semiclassicalJacobiBasis_apply,
    semiclassicalJacobiBasis, Sequence.basis_eq_self] using
    (semiclassicalJacobiSequence_spec ht hα' hβ' γ').1.connection_diag_pos
      (hasFiniteMoments_semiclassicalJacobiMeasure ht hα' hβ')
      (hasPositiveDefiniteIntegral_semiclassicalJacobiMeasure ht hα' hβ')
      (fun n ↦ isUnit_iff_ne_zero.mpr ((semiclassicalJacobiSequence_spec ht hα' hβ' γ').2.1 n).ne')
      (semiclassicalJacobiSequence ht hα hβ γ)
      (semiclassicalJacobiSequence_spec ht hα hβ γ).2.1
      (semiclassicalJacobiSequence_spec ht hα' hβ' γ').2.1 n

/-- Every leading section is the transpose-product factor of the target Gram matrix. -/
theorem semiclassicalJacobiConnection_gram {t α β α' β' : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα' : -1 < α') (hβ' : -1 < β') (γ' : ℝ) (N : ℕ) :
    Matrix.of (fun i j : Fin N ↦ integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ α' β' γ') (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ α' β' γ') (semiclassicalJacobiSequence ht hα hβ γ i) (semiclassicalJacobiSequence ht hα hβ γ j))) =
      (Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiConnection ht hα hβ γ hα' hβ' γ' i j)).transpose *
        Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiConnection ht hα hβ γ hα' hβ' γ' i j) := by
  simpa only [semiclassicalJacobiConnection, Basis.toMatrix_apply, semiclassicalJacobiBasis_apply,
    semiclassicalJacobiBasis, Sequence.basis_eq_self] using
    (semiclassicalJacobiSequence_spec ht hα' hβ' γ').1.connection_gram
      (hasFiniteMoments_semiclassicalJacobiMeasure ht hα' hβ')
      (hasPositiveDefiniteIntegral_semiclassicalJacobiMeasure ht hα' hβ')
      (fun n ↦ isUnit_iff_ne_zero.mpr ((semiclassicalJacobiSequence_spec ht hα' hβ' γ').2.1 n).ne')
      (semiclassicalJacobiSequence ht hα hβ γ)
      (fun n ↦ by simpa only [integralInner_self] using
        (semiclassicalJacobiSequence_spec ht hα' hβ' γ').2.2 n) N

/-- The cubic multiplier in the semiclassical Jacobi Pearson equation. -/
noncomputable def semiclassicalJacobiSigma (t : ℝ) : ℝ[X] := X * (1 - X) * (C t - X)

/-- The quadratic coefficient in the semiclassical Jacobi Pearson equation. -/
noncomputable def semiclassicalJacobiTau (t α β γ : ℝ) : ℝ[X] :=
  C (α + 1) * (1 - X) * (C t - X) - C (β + 1) * X * (C t - X) -
    C (γ + 1) * X * (1 - X)

/-- In the integrable endpoint range the shifted weight has a continuous
closed-interval extension with both endpoint values zero. -/
theorem exists_continuous_semiclassicalJacobiWeight_add_one (t : Ioi (1 : ℝ))
    {α β γ : ℝ} (hα : -1 < α) (hβ : -1 < β) :
    ∃ F : C(Icc (0 : ℝ) 1, ℝ),
      (∀ x : Ioo (0 : ℝ) 1, F ⟨x, x.property.1.le, x.property.2.le⟩ =
        (semiclassicalJacobiWeight t (α + 1) (β + 1) (γ + 1) x : ℝ)) ∧
      F ⟨0, by constructor <;> norm_num⟩ = 0 ∧ F ⟨1, by constructor <;> norm_num⟩ = 0 := by
  let F : C(Icc (0 : ℝ) 1, ℝ) := ⟨fun x ↦
    (x : ℝ) ^ (α + 1) * (1 - (x : ℝ)) ^ (β + 1) * ((t : ℝ) - x) ^ (γ + 1),
    ((continuous_subtype_val.rpow_const (fun _ ↦ Or.inr (by linarith : 0 ≤ α + 1))).mul
      ((continuous_const.sub continuous_subtype_val).rpow_const
        (fun _ ↦ Or.inr (by linarith : 0 ≤ β + 1)))).mul
      ((continuous_const.sub continuous_subtype_val).rpow_const
        (fun x ↦ Or.inl (by change (t : ℝ) - x ≠ 0; have ht : 1 < (t : ℝ) := t.property; linarith [x.property.2])))⟩
  refine ⟨F, fun _ ↦ rfl, ?_, ?_⟩
  · change (0 : ℝ)^(α+1) * (1-0)^(β+1) * ((t : ℝ)-0)^(γ+1) = 0
    simp [Real.zero_rpow (by linarith : α+1 ≠ 0)]
  · change (1 : ℝ)^(α+1) * (1-1)^(β+1) * ((t : ℝ)-1)^(γ+1) = 0
    simp [Real.zero_rpow (by linarith : β+1 ≠ 0)]

/-- The left endpoint is a limit; no restriction on the opposite exponent is needed. -/
theorem tendsto_semiclassicalJacobiWeight_add_one_zero (t : Ioi (1 : ℝ)) (β γ : ℝ)
    {α : ℝ} (hα : -1 < α) :
    Tendsto (fun x : Ioo (0 : ℝ) 1 ↦
      (semiclassicalJacobiWeight t (α + 1) (β + 1) (γ + 1) x : ℝ))
      (Filter.comap Subtype.val (𝓝 (0 : ℝ))) (𝓝 0) := by
  have hval : Tendsto (Subtype.val : Ioo (0 : ℝ) 1 → ℝ)
      (Filter.comap Subtype.val (𝓝 0)) (𝓝 0) := Filter.tendsto_comap
  let u : Ioo (0 : ℝ) 1 → ℝ≥0 := fun x ↦ ⟨x, x.property.1.le⟩
  have hu : Tendsto u (Filter.comap Subtype.val (𝓝 0)) (𝓝 0) :=
    NNReal.tendsto_coe.mp hval
  have h1 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ (x : ℝ) ^ (α + 1))
      (Filter.comap Subtype.val (𝓝 0)) (𝓝 0) := by
    have h := ((NNReal.continuous_rpow_const (by linarith : (0 : ℝ) ≤ α + 1)).tendsto 0).comp hu
    rw [NNReal.zero_rpow (by linarith : α + 1 ≠ 0)] at h
    exact NNReal.tendsto_coe.mpr h
  have h2 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ (1 - (x : ℝ)) ^ (β + 1))
      (Filter.comap Subtype.val (𝓝 0)) (𝓝 ((1 - 0 : ℝ) ^ (β + 1))) :=
    ((continuousAt_const.sub continuousAt_id).rpow_const (Or.inl (by norm_num))).tendsto.comp hval
  have h3 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ ((t : ℝ) - x) ^ (γ + 1))
      (Filter.comap Subtype.val (𝓝 0)) (𝓝 (((t : ℝ) - 0) ^ (γ + 1))) :=
    ((continuousAt_const.sub continuousAt_id).rpow_const
      (Or.inl (by have ht : 1 < (t : ℝ) := t.property; dsimp; linarith))).tendsto.comp hval
  simpa only [semiclassicalJacobiWeight_apply, zero_mul] using (h1.mul h2).mul h3

/-- The right endpoint is a limit; the left exponent can be arbitrary. -/
theorem tendsto_semiclassicalJacobiWeight_add_one_one (t : Ioi (1 : ℝ)) (α γ : ℝ)
    {β : ℝ} (hβ : -1 < β) :
    Tendsto (fun x : Ioo (0 : ℝ) 1 ↦
      (semiclassicalJacobiWeight t (α + 1) (β + 1) (γ + 1) x : ℝ))
      (Filter.comap Subtype.val (𝓝 (1 : ℝ))) (𝓝 0) := by
  have hval : Tendsto (Subtype.val : Ioo (0 : ℝ) 1 → ℝ)
      (Filter.comap Subtype.val (𝓝 1)) (𝓝 1) := Filter.tendsto_comap
  let u : Ioo (0 : ℝ) 1 → ℝ≥0 := fun x ↦ ⟨1 - (x : ℝ), by linarith [x.property.2]⟩
  have hu : Tendsto u (Filter.comap Subtype.val (𝓝 1)) (𝓝 0) := by
    rw [← NNReal.tendsto_coe]
    have h := (tendsto_const_nhds (x := (1 : ℝ))).sub hval
    rw [sub_self] at h
    exact h
  have h1 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ (x : ℝ) ^ (α + 1))
      (Filter.comap Subtype.val (𝓝 1)) (𝓝 ((1 : ℝ) ^ (α + 1))) :=
    (continuousAt_id.rpow_const (Or.inl (by norm_num))).tendsto.comp hval
  have h2 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ (1 - (x : ℝ)) ^ (β + 1))
      (Filter.comap Subtype.val (𝓝 1)) (𝓝 0) := by
    have h := ((NNReal.continuous_rpow_const (by linarith : (0 : ℝ) ≤ β + 1)).tendsto 0).comp hu
    rw [NNReal.zero_rpow (by linarith : β + 1 ≠ 0)] at h
    exact NNReal.tendsto_coe.mpr h
  have h3 : Tendsto (fun x : Ioo (0 : ℝ) 1 ↦ ((t : ℝ) - x) ^ (γ + 1))
      (Filter.comap Subtype.val (𝓝 1)) (𝓝 (((t : ℝ) - 1) ^ (γ + 1))) :=
    ((continuousAt_const.sub continuousAt_id).rpow_const
      (Or.inl (by have ht : 1 < (t : ℝ) := t.property; dsimp; linarith))).tendsto.comp hval
  simpa only [semiclassicalJacobiWeight_apply, mul_zero, zero_mul] using (h1.mul h2).mul h3

/-- The Pearson identity tested against any polynomial, with zero endpoint terms. -/
theorem integral_semiclassicalJacobi_pearson {t α β γ : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (p : ℝ[X]) :
    ∫ x, (semiclassicalJacobiSigma t * p.derivative + semiclassicalJacobiTau t α β γ * p).eval x
      ∂semiclassicalJacobiMeasure ⟨t, by assumption⟩ α β γ = 0 := by
  let q := semiclassicalJacobiSigma t * p.derivative + semiclassicalJacobiTau t α β γ * p
  obtain ⟨F, hF, h0, h1⟩ := exists_continuous_semiclassicalJacobiWeight_add_one
    (γ := γ) ⟨t, ht⟩ hα hβ
  let P : C(Icc (0 : ℝ) 1, ℝ) := ⟨fun x ↦ F x * p.eval (x : ℝ),
    F.continuous.mul (p.continuous.comp continuous_subtype_val)⟩
  have hq := integrable_semiclassicalJacobiMeasure (γ := γ) ht hα hβ q.continuous.continuousOn
  have hg := (integrable_semiclassicalJacobiMeasure_iff ⟨t, ht⟩ α β γ _).mp hq
  have hd (x : Ioo (0 : ℝ) 1) :
      HasDerivAtOnOpen ⟨Ioo (0 : ℝ) 1, isOpen_Ioo⟩
        (fun y ↦ P ⟨y, ⟨y.property.1.le, y.property.2.le⟩⟩)
        ((semiclassicalJacobiWeight ⟨t, ht⟩ α β γ x : ℝ) * q.eval (x : ℝ)) x := by
    have he : (fun y : Ioo (0 : ℝ) 1 ↦ P ⟨y, ⟨y.property.1.le, y.property.2.le⟩⟩) =
        fun y ↦ (semiclassicalJacobiWeight ⟨t, ht⟩ (α+1) (β+1) (γ+1) y : ℝ) * p.eval (y : ℝ) := by
      funext y
      exact congrArg (fun r ↦ r * p.eval (y : ℝ)) (hF y)
    erw [he]
    convert ((hasDerivAt_semiclassicalJacobiWeight_add_one ⟨t, ht⟩ α β γ x).mul
      ⟨_, fun _ ↦ rfl, p.hasDerivAt (x : ℝ)⟩) using 1
    rw [semiclassicalJacobiWeight_add_one]
    simp only [q, semiclassicalJacobiSigma, semiclassicalJacobiTau, eval_add, eval_mul,
      eval_sub, eval_one, eval_X, eval_C]
    ring
  have hi := intervalIntegral.integral_eq_sub_of_hasDerivAtOnOpen
    (show (0 : ℝ) < 1 by norm_num) P ⟨_, hg⟩ hd
  change (∫ x : Ioo (0 : ℝ) 1,
    (semiclassicalJacobiWeight ⟨t, ht⟩ α β γ x : ℝ) * q.eval (x : ℝ)) =
    F ⟨1, by norm_num⟩ * p.eval 1 - F ⟨0, by norm_num⟩ * p.eval 0 at hi
  rw [h0, h1, zero_mul, zero_mul, sub_self] at hi
  exact (integral_semiclassicalJacobiMeasure ⟨t, ht⟩ α β γ ⟨_, hq⟩).trans hi

/-- The derivative matrix from the original to the simultaneously raised family. -/
noncomputable def semiclassicalJacobiDerivative {t α β : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα1 : -1 < α+1) (hβ1 : -1 < β+1) : Matrix ℕ ℕ ℝ :=
  (semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)).toMatrix
    (fun n ↦ (semiclassicalJacobiSequence ht hα hβ γ n).derivative)

/-- Column n of the actual derivative matrix is supported only at n-1 and n-2. -/
theorem semiclassicalJacobiDerivative_eq_zero {t α β : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα1 : -1 < α+1) (hβ1 : -1 < β+1) (n j : ℕ)
    (hgap : n ≤ j ∨ j+2 < n) :
    semiclassicalJacobiDerivative ht hα hβ γ hα1 hβ1 j n = 0 := by
  rw [semiclassicalJacobiDerivative, Basis.toMatrix_apply, semiclassicalJacobiBasis_repr]
  exact integralInner_semiclassicalJacobi_derivative_eq_zero ht hα hβ hα1 hβ1 n j hgap

/-- The coefficient operator differentiates every finite polynomial expansion exactly. -/
theorem semiclassicalJacobi_derivative_coefficientMap_exact {t α β : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα1 : -1 < α+1) (hβ1 : -1 < β+1) (v : ℕ →₀ ℝ) :
    (semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)).repr.symm
      ((semiclassicalJacobiBasis ht hα hβ γ).coefficientMap
        (semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)) derivative v) =
      ((semiclassicalJacobiBasis ht hα hβ γ).repr.symm v).derivative :=
  Basis.repr_symm_coefficientMap _ _ _ _

/-- The sparse matrix acts by a finite sum over input coefficients. -/
theorem semiclassicalJacobi_derivative_coefficientMap_apply {t α β : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ)
    (hα1 : -1 < α+1) (hβ1 : -1 < β+1) (v : ℕ →₀ ℝ) (i : ℕ) :
    (semiclassicalJacobiBasis ht hα hβ γ).coefficientMap
        (semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)) derivative v i =
      v.sum (fun j a ↦ semiclassicalJacobiDerivative ht hα hβ γ hα1 hβ1 i j * a) := by
  rw [Basis.coefficientMap_apply]
  simp only [semiclassicalJacobiDerivative, semiclassicalJacobiBasis_apply]

/-- The exact polynomial action of weighted differentiation. -/
noncomputable def semiclassicalJacobiWeightedDerivativeMap (t α β γ : ℝ) : ℝ[X] →ₗ[ℝ] ℝ[X] :=
  (LinearMap.mulLeft ℝ (semiclassicalJacobiSigma t)).comp derivative +
    LinearMap.mulLeft ℝ (semiclassicalJacobiTau t α β γ)

/-- The weighted derivative matrix from the raised family back to the original family. -/
noncomputable def semiclassicalJacobiWeightedDerivative {t α β : ℝ} (ht : 1<t)
    (hα : -1<α) (hβ : -1<β) (γ : ℝ) (hα1 : -1<α+1) (hβ1 : -1<β+1) : Matrix ℕ ℕ ℝ :=
  (semiclassicalJacobiBasis ht hα hβ γ).toMatrix (fun j ↦
    semiclassicalJacobiWeightedDerivativeMap t α β γ
      (semiclassicalJacobiSequence ht hα1 hβ1 (γ+1) j))

/-- Column j of weighted differentiation has support only at j+1 and j+2. -/
theorem semiclassicalJacobiWeightedDerivative_eq_zero {t α β : ℝ} (ht : 1<t)
    (hα : -1<α) (hβ : -1<β) (γ : ℝ) (hα1 : -1<α+1) (hβ1 : -1<β+1)
    (n j : ℕ) (hgap : n≤j ∨ j+2<n) :
    semiclassicalJacobiWeightedDerivative ht hα hβ γ hα1 hβ1 n j=0 := by
  rw [semiclassicalJacobiWeightedDerivative_eq_neg_transpose]
  simp only [Matrix.neg_apply, Matrix.transpose_apply,
    semiclassicalJacobiDerivative_eq_zero ht hα hβ γ hα1 hβ1 n j hgap, neg_zero]

/-- Weighted differentiation acts exactly on every finite expansion. -/
theorem semiclassicalJacobi_weightedDerivative_coefficientMap_exact {t α β : ℝ} (ht : 1<t)
    (hα : -1<α) (hβ : -1<β) (γ : ℝ) (hα1 : -1<α+1) (hβ1 : -1<β+1) (v : ℕ →₀ ℝ) :
    (semiclassicalJacobiBasis ht hα hβ γ).repr.symm
      ((semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)).coefficientMap
        (semiclassicalJacobiBasis ht hα hβ γ) (semiclassicalJacobiWeightedDerivativeMap t α β γ) v) =
      semiclassicalJacobiWeightedDerivativeMap t α β γ
        ((semiclassicalJacobiBasis ht hα1 hβ1 (γ+1)).repr.symm v) :=
  Basis.repr_symm_coefficientMap _ _ _ _

/-- The unit-shift connection is the unique positive upper triangular Cholesky factor. -/
theorem semiclassicalJacobiRaisingSection_cholesky_unique {t α β : ℝ} (ht : 1 < t)
    (hα : -1 < α) (hβ : -1 < β) (γ : ℝ) (N : ℕ)
    (R : Matrix (Fin N) (Fin N) ℝ) (hR : R.IsUpperTriangular)
    (hdiag : ∀ i, 0 < R i i)
    (hgram : R.transpose * R = t • (1 : Matrix (Fin N) (Fin N) ℝ) -
      Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiMatrix ht hα hβ γ i j)) :
    R = semiclassicalJacobiRaisingSection ht hα hβ γ 1 N :=
  Matrix.eq_of_pos_upperTriangular_gram_eq hR
    (semiclassicalJacobiRaisingSection_triangular ht hα hβ γ 1 N).1 hdiag
    (semiclassicalJacobiRaisingSection_triangular ht hα hβ γ 1 N).2
    (hgram.trans (semiclassicalJacobiRaisingSection_cholesky ht hα hβ γ N))

/-- The padded intertwining identity gives the Jacobi update using a genuine
finite inverse. The inverse is specified by both inverse equations. -/
theorem semiclassicalJacobi_padded_update (γ : ℝ) (c N : ℕ) :
    ∃ S : Matrix (Fin N) (Fin N) ℝ,
      semiclassicalJacobiRaisingSection ht hα hβ γ c N * S = 1 ∧
      S * semiclassicalJacobiRaisingSection ht hα hβ γ c N = 1 ∧
      Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiMatrix ht hα hβ (γ+c) i j) =
        (Matrix.of (fun i : Fin N ↦ fun k : Fin (N+1) ↦
          semiclassicalJacobiConnection ht hα hβ γ hα hβ (γ+c) i k)) *
        (Matrix.of (fun k : Fin (N+1) ↦ fun j : Fin N ↦
          semiclassicalJacobiMatrix ht hα hβ γ k j)) * S := by
  obtain ⟨u, hu⟩ := semiclassicalJacobiRaisingSection_isUnit ht hα hβ γ c N
  refine ⟨↑u⁻¹, ?_, ?_, ?_⟩
  · rw [← hu]; exact u.mul_inv
  · rw [← hu]; exact u.inv_mul
  · rw [semiclassicalJacobi_padded_intertwining]
    change _ = _ * semiclassicalJacobiRaisingSection ht hα hβ γ c N * _
    rw [← hu, Units.mul_inv_cancel_right]

/-- The rectangular Q has orthonormal columns. -/
theorem semiclassicalJacobiQRSection_orthogonal (γ : ℝ) (N : ℕ) :
    (semiclassicalJacobiQRSection ht hα hβ γ N).transpose *
      semiclassicalJacobiQRSection ht hα hβ γ N = 1 := by
  ext i j
  simp only [Matrix.mul_apply, Matrix.transpose_apply, semiclassicalJacobiQRSection]
  rw [← semiclassicalJacobi_integralInner_eq_sum ht hα hβ γ _ _ (N+1)
    (by have hd := natDegree_mul_le (p := C t-X) (q := semiclassicalJacobiSequence ht hα hβ (γ+2) i)
        have htdeg : (C t-X).natDegree ≤ 1 := by compute_degree
        simp only [Sequence.natDegree_eq] at hd
        omega), semiclassicalJacobi_integralInner_shift_two ht hα hβ]
  rw [← semiclassicalJacobiBasis_repr ht hα hβ, ← semiclassicalJacobiBasis_apply ht hα hβ (γ+2) i,
    Basis.repr_self]
  simp [Matrix.one_apply, Finsupp.single_apply, Fin.ext_iff]

/-- The twice-raised Jacobi section is the compression by the rectangular Q.
All entries are exact; the old Jacobi section has one extra row and column. -/
theorem semiclassicalJacobiQRSection_update (γ : ℝ) (N : ℕ) :
    (semiclassicalJacobiQRSection ht hα hβ γ N).transpose *
      (Matrix.of (fun i j : Fin (N+1) ↦ semiclassicalJacobiMatrix ht hα hβ γ i j)) *
      semiclassicalJacobiQRSection ht hα hβ γ N =
    Matrix.of (fun i j : Fin N ↦ semiclassicalJacobiMatrix ht hα hβ (γ+2) i j) := by
  ext i j
  simp only [Matrix.mul_apply, Matrix.transpose_apply, Matrix.of_apply,
    semiclassicalJacobiQRSection, semiclassicalJacobiMatrix, Basis.toMatrix_apply]
  simp_rw [Finset.sum_mul]
  rw [Finset.sum_comm]
  have hi : ((C t-X)*semiclassicalJacobiSequence ht hα hβ (γ+2) i).natDegree < N+1 := by
    have hd := natDegree_mul_le (p := C t-X) (q := semiclassicalJacobiSequence ht hα hβ (γ+2) i)
    have htdeg : (C t-X).natDegree ≤ 1 := by compute_degree
    simp only [Sequence.natDegree_eq] at hd
    omega
  have hj : ((C t-X)*semiclassicalJacobiSequence ht hα hβ (γ+2) j).natDegree < N+1 := by
    have hd := natDegree_mul_le (p := C t-X) (q := semiclassicalJacobiSequence ht hα hβ (γ+2) j)
    have htdeg : (C t-X).natDegree ≤ 1 := by compute_degree
    simp only [Sequence.natDegree_eq] at hd
    omega
  rw [← semiclassicalJacobi_integralInner_mul_eq_sum ht hα hβ γ _ _ X (N+1) hi hj]
  have he : X * ((C t-X)*semiclassicalJacobiSequence ht hα hβ (γ+2) j) =
      (C t-X)*(X*semiclassicalJacobiSequence ht hα hβ (γ+2) j) := by ring
  have htransfer := semiclassicalJacobi_integralInner_shift_two ht hα hβ γ
    (semiclassicalJacobiSequence ht hα hβ (γ+2) i)
    (X*semiclassicalJacobiSequence ht hα hβ (γ+2) j)
  rw [semiclassicalJacobiBasis_repr ht hα hβ]
  simp only [integralInner] at htransfer ⊢
  rw [he, htransfer]
  apply integral_congr_ae
  filter_upwards [] with x
  exact mul_comm _ _