/-- The integral Gram matrix, without installing a measure-dependent global
inner product instance on polynomials. -/
noncomputable def integralGram {ι : Type*} (μ : Measure (σ → ℝ)) (P : ι → MvPolynomial σ ℝ)
    (hP : ∀ i j, Integrable (fun x => eval x (P i) * eval x (P j)) μ := by
      intro i j
      aesop (config := { enableSimp := false }) (rule_sets := [PolynomialIntegrability])) :
    Matrix ι ι ℝ := fun i j => integralInner μ (IntegrablePair.of μ (P i) (P j))

/-- Every ordered finite polynomial basis admits orthonormalization with an
upper triangular, positive-diagonal connection matrix. -/
theorem exists_orthonormal_pos_upperTriangular_basis
    (μ : Measure (σ → ℝ)) (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ)
    {V : Submodule ℝ (MvPolynomial σ ℝ)} (b : Basis (Fin n) ℝ V) :
    ∃ q : Basis (Fin n) ℝ V,
      integralGram μ (fun i => (q i : MvPolynomial σ ℝ)) = 1 ∧
        (q.toMatrix b).IsUpperTriangular ∧ (∀ i, 0 < q.toMatrix b i i) := by
  let core := integralInnerProductCore μ hμ hpos
  let : NormedAddCommGroup (MvPolynomial σ ℝ) :=
    InnerProductSpace.Core.toNormedAddCommGroup (𝕜 := ℝ)
  let : InnerProductSpace ℝ (MvPolynomial σ ℝ) :=
    InnerProductSpace.ofCore (inferInstance : PreInnerProductSpace.Core ℝ (MvPolynomial σ ℝ))
  obtain ⟨q, ht, hp⟩ := InnerProductSpace.exists_orthonormalBasis_pos_upperTriangular b
  exact ⟨q.toBasis, Matrix.gram_eq_one_iff_orthonormal.mpr q.orthonormal, ht, hp⟩

/-- The change of basis to any orthonormal polynomial basis factors the
modified-measure Gram matrix, in any finite number of variables. -/
theorem integralGram_eq_transpose_toMatrix_mul {ι : Type*} [Fintype ι] [DecidableEq ι]
    (μ : Measure (σ → ℝ)) (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ)
    {V : Submodule ℝ (MvPolynomial σ ℝ)} (q b : Basis ι ℝ V)
    (hq : integralGram μ (fun i => (q i : MvPolynomial σ ℝ)) = 1) :
    integralGram μ (fun i => (b i : MvPolynomial σ ℝ)) =
      (q.toMatrix b).transpose * q.toMatrix b := by
  let core := integralInnerProductCore μ hμ hpos
  let : NormedAddCommGroup (MvPolynomial σ ℝ) :=
    InnerProductSpace.Core.toNormedAddCommGroup (𝕜 := ℝ)
  let : InnerProductSpace ℝ (MvPolynomial σ ℝ) :=
    InnerProductSpace.ofCore (inferInstance : PreInnerProductSpace.Core ℝ (MvPolynomial σ ℝ))
  have ho : Orthonormal ℝ q := Matrix.gram_eq_one_iff_orthonormal.mp hq
  change Matrix.gram ℝ b = _
  simpa only [Basis.toBasis_toOrthonormalBasis] using
    InnerProductSpace.gram_eq_transpose_toMatrix_mul (q.toOrthonormalBasis ho) b

/-- A nonzero nonnegative polynomial modification preserves strict positivity
of the polynomial pairing, even if its density vanishes on part of the support. -/
theorem HasPositiveDefiniteIntegral.withDensity_polynomial
    (hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
    (u : MvPolynomial σ ℝ) (hu0 : u ≠ 0) (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u) :
    HasPositiveDefiniteIntegral (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)) := by
  apply HasPositiveDefiniteIntegral.of_positive
  intro p hp
  rw [← integralInner_self _ p
    (h := (hμ.withDensity_polynomial u hu).integrable_mul p p),
    integralInner_withDensity_polynomial u p p hu (hpq := hμ.integrable_mul _ _)]
  have hi := hμ.integrable_mul (u * p) p
  have hn : 0 ≤ᵐ[μ] fun x => eval x (u * p) * eval x p := by
    filter_upwards [hu] with x hx
    change 0 ≤ eval x (u * p) * eval x p
    simpa only [map_mul, mul_assoc, ← pow_two] using mul_nonneg hx (sq_nonneg (eval x p))
  apply lt_of_le_of_ne (integral_nonneg_of_ae hn)
  intro heq
  have hz := (integral_eq_zero_iff_of_nonneg_ae hn hi).mp heq.symm
  have hup : integralInner μ (IntegrablePair.of μ (u * p) (u * p)) = 0 := by
    unfold integralInner
    calc
      _ = ∫ x, (0 : ℝ) ∂μ := integral_congr_ae (by
        filter_upwards [hz] with x hx
        simp only [Pi.zero_apply, map_mul] at hx ⊢
        rcases mul_eq_zero.mp hx with h | h
        · rw [h, zero_mul]
        · simp [h])
      _ = 0 := integral_zero _ _
  exact (mul_ne_zero hu0 hp) (hpos.eq_zero_of_integralInner_self_eq_zero hup)

/-- The modified pairing is polynomial multiplication in the original pairing. -/
theorem integralInner_withDensity_polynomial (u p q : MvPolynomial σ ℝ)
    (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u)
    {hpq : Integrable (fun x => eval x (u * p) * eval x q) μ} :
    integralInner (polynomialDensityMeasure u (μ.comap Subtype.val))
      (IntegrablePair.of _ p q (h := (integrable_mul_withDensity_polynomial_iff u p q hu).2 hpq)) =
      integralInner μ (IntegrablePair.of μ (u * p) q) := by
  obtain ⟨w, hw, hwval, _, he⟩ := exists_polynomialDensityMeasure_extension u
  change (∫ x, eval x p * eval x q ∂polynomialDensityMeasure u (μ.comap Subtype.val)) =
    ∫ x, eval x (u*p) * eval x q ∂μ
  rw [he, integral_withDensity_eq_integral_toReal_smul hw.coe_nnreal_ennreal
    (Filter.Eventually.of_forall (fun _ => ENNReal.coe_lt_top))]
  apply integral_congr_ae
  filter_upwards [hu] with x hx
  rw [hwval x hx]
  simp only [ENNReal.coe_toReal, NNReal.coe_mk, smul_eq_mul, map_mul]
  ring

/-- The modified finite Gram matrix consists of exact multiplication pairings;
this holds on arbitrary multivariate polynomial families. -/
theorem integralGram_withDensity_polynomial {n : ℕ}
    (u : MvPolynomial σ ℝ) (P : Fin n → MvPolynomial σ ℝ)
    (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u)
    (hi : ∀ i j, Integrable (fun x => eval x (u*P i)*eval x (P j)) μ := by
      intro i j
      aesop (config := { enableSimp := false }) (rule_sets := [PolynomialIntegrability])) :
    integralGram (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)) P
      (hP := fun i j => (integrable_mul_withDensity_polynomial_iff u (P i) (P j) hu).2 (hi i j)) =
      fun i j => integralInner μ (IntegrablePair.of μ (u * P i) (P j)) := by
  ext i j
  exact integralInner_withDensity_polynomial u (P i) (P j) hu (hpq := hi i j)

/-- A degree-$d$ modification has no Gram entries between polynomial degree
blocks separated by more than $d$. This does not assert sparsity within blocks. -/
theorem integralInner_withDensity_eq_zero_of_degree
    {ι : Type*} (P : ι → MvPolynomial σ ℝ) (degree : ι → ℕ)
    (hdegree : ∀ i, (P i).totalDegree ≤ degree i)
    (horth : ∀ i q, q.totalDegree < degree i →
      ∃ hi : Integrable (fun x => eval x q * eval x (P i)) μ,
        integralInner μ (IntegrablePair.of μ q (P i) hi) = 0)
    (u : MvPolynomial σ ℝ) (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u) (i j : ι)
    (hij : degree j + u.totalDegree < degree i ∨
      degree i + u.totalDegree < degree j) :
    integralInner (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)) (IntegrablePair.of (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)) (P i) (P j) (h := by
      apply (integrable_mul_withDensity_polynomial_iff u (P i) (P j) hu).2
      rcases hij with hij | hij
      · have hd : (u*P j).totalDegree < degree i := by
          have := totalDegree_mul u (P j)
          have := hdegree j
          omega
        obtain ⟨hi, _⟩ := horth i (u*P j) hd
        simpa only [map_mul, mul_comm, mul_left_comm, mul_assoc] using hi
      · have hd : (u*P i).totalDegree < degree j := by
          have := totalDegree_mul u (P i)
          have := hdegree i
          omega
        exact (horth j (u*P i) hd).choose)) = 0 := by
  rcases hij with hij | hij
  · have hd : (u*P j).totalDegree < degree i := by
      have := totalDegree_mul u (P j)
      have := hdegree j
      omega
    obtain ⟨hi, hz⟩ := horth i (u*P j) hd
    rw [integralInner_comm, integralInner_withDensity_polynomial u (P j) (P i) hu (hpq := hi)]
    exact hz
  · have hd : (u*P i).totalDegree < degree j := by
      have := totalDegree_mul u (P i)
      have := hdegree i
      omega
    obtain ⟨hi, hz⟩ := horth j (u*P i) hd
    rw [integralInner_withDensity_polynomial u (P i) (P j) hu (hpq := hi)]
    exact hz

/-- The positive Cholesky factor of a polynomially modified Gram matrix is
banded by total degree. Additional zeros inside those blocks are a separate issue. -/
theorem gram_factor_eq_zero_of_degree {n : ℕ}
    (P : Fin n → MvPolynomial σ ℝ) (degree : Fin n → ℕ)
    (hmon : Monotone degree) (hdegree : ∀ i, (P i).totalDegree ≤ degree i)
    (horth : ∀ i q, q.totalDegree < degree i →
      ∃ hi : Integrable (fun x => eval x q * eval x (P i)) μ,
        integralInner μ (IntegrablePair.of μ q (P i) hi) = 0)
    (u : MvPolynomial σ ℝ) (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u)
    (R : Matrix (Fin n) (Fin n) ℝ) (hR : R.IsUpperTriangular)
    (hp : ∀ i, 0 < R i i)
    (hi : ∀ i j, Integrable (fun x => eval x (P i)*eval x (P j))
      (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)))
    (hg : integralGram (MvPolynomial.polynomialDensityMeasure (u) ((μ).comap Subtype.val)) P (hP := hi) =
      R.transpose * R) :
    ∀ i j, degree i + u.totalDegree < degree j → R i j = 0 := by
  apply hR.eq_zero_of_transpose_mul_self_degree_band (fun i => ne_of_gt (hp i))
    degree hmon u.totalDegree
  intro i j hij
  rw [← hg]
  exact integralInner_withDensity_eq_zero_of_degree P degree hdegree horth u hu i j (Or.inr hij)

/-- A separated polynomial density gives the product of the modified factor pairings. -/
theorem integralInner_tensorProductModification (μ : σ → Measure ℝ)
    [∀ i, SigmaFinite (μ i)] (h p q : σ → ℝ[X])
    (hh : ∀ i, ∀ᵐ x ∂μ i, 0 ≤ (h i).eval x)
    (hu : ∀ᵐ x ∂Measure.pi μ, 0 ≤ eval x (tensorProductPolynomial h))
    (hi : ∀ i, Integrable (fun x => (h i*p i).eval x*(q i).eval x) (μ i) := by
      intro i
      aesop (config := { enableSimp := false }) (rule_sets := [PolynomialIntegrability])) :
    integralInner (MvPolynomial.polynomialDensityMeasure ((tensorProductPolynomial h)) (((Measure.pi μ)).comap Subtype.val)) (IntegrablePair.of (MvPolynomial.polynomialDensityMeasure ((tensorProductPolynomial h)) (((Measure.pi μ)).comap Subtype.val)) (tensorProductPolynomial p) (tensorProductPolynomial q) (h := (MvPolynomial.integrable_mul_withDensity_polynomial_iff _ _ _ hu).2 (by simpa only [tensorProductPolynomial_mul] using integrable_mul_tensorProductPolynomial μ (fun i => h i*p i) q hi))) =
      ∏ i, Polynomial.integralInner (Polynomial.polynomialDensityMeasure ((h i)) (((μ i)).comap Subtype.val)) (Polynomial.IntegrablePair.of (Polynomial.polynomialDensityMeasure ((h i)) (((μ i)).comap Subtype.val)) (p i) (q i) (h := (Polynomial.integrable_polynomialDensityMeasure_comap_mul_iff (h i) (p i) (q i) (hh i)).2 (hi i))) := by
  have hpq : Integrable (fun x => eval x (tensorProductPolynomial h*tensorProductPolynomial p)*eval x (tensorProductPolynomial q)) (Measure.pi μ) := by
    simpa only [tensorProductPolynomial_mul] using
      integrable_mul_tensorProductPolynomial μ (fun i => h i*p i) q hi
  rw [MvPolynomial.integralInner_withDensity_polynomial _ _ _ hu (hpq := hpq)]
  simp only [← tensorProductPolynomial_mul]
  rw [integralInner_tensorProductPolynomial μ (fun i => h i*p i) q (h := hi)]
  apply Finset.prod_congr rfl
  intro i _
  exact (Polynomial.integralInner_polynomialDensityMeasure_comap (h i) (p i) (q i) (hh i)
    (hpq := hi i)).symm

/-- Leaving the univariate modification band in any coordinate makes the
multivariate connection coefficient exactly zero. -/
theorem integralInner_tensorProductModification_eq_zero_of_gap (μ : σ → Measure ℝ)
    (hμ : ∀ i, Polynomial.HasFiniteMoments (μ i)) (h : σ → ℝ[X])
    (hh : ∀ i, ∀ᵐ x ∂μ i, 0 ≤ (h i).eval x)
    (hu : ∀ᵐ x ∂Measure.pi μ, 0 ≤ eval x (tensorProductPolynomial h))
    (P Q : σ → Polynomial.Sequence ℝ)
    (hP : ∀ i, (P i).IsOrthogonal (μ i))
    (hQ : ∀ i, (Q i).IsOrthogonal (Polynomial.polynomialDensityMeasure ((h i)) (((μ i)).comap Subtype.val)))
    (n m : σ → ℕ) (i : σ) (hgap : n i < m i ∨ m i + (h i).natDegree < n i) :
    integralInner (MvPolynomial.polynomialDensityMeasure ((tensorProductPolynomial h)) (((Measure.pi μ)).comap Subtype.val)) (IntegrablePair.of (MvPolynomial.polynomialDensityMeasure ((tensorProductPolynomial h)) (((Measure.pi μ)).comap Subtype.val)) (tensorProductPolynomial (fun j ↦ P j (n j))) (tensorProductPolynomial (fun j ↦ Q j (m j)))) = 0 := by
  let : ∀ i, IsFiniteMeasure (μ i) := fun i ↦ (hμ i).isFiniteMeasure
  rw [integralInner_tensorProductModification μ h _ _ hh hu]
  apply Finset.prod_eq_zero (Finset.mem_univ i)
  exact (hP i).integralInner_eq_zero_of_polynomial_modification (hμ i) (h i) (hh i) (hQ i) hgap

/-- The complete degree basis rescaled by positive inverse norms. -/
noncomputable def normalizedProriolDegreeBasis {α β γ : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (hγ : -1 < γ) (N : ℕ) :
    Basis {d : Fin 2 →₀ ℕ // d.sum (fun _ e ↦ e) ≤ N} ℝ
      (restrictTotalDegree (Fin 2) ℝ N) :=
  (proriolDegreeBasis hα hβ hγ N).isUnitSMul (fun d ↦ isUnit_iff_ne_zero.mpr
    (inv_ne_zero (proriolNorm_pos hα hβ hγ (d.1 0+d.1 1) (d.1 1)).ne'))

/-- The normalized complete degree basis in degree/inner order. -/
noncomputable def orderedProriolBasis {α β γ : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (hγ : -1 < γ) (N : ℕ) :
    Basis (Fin (Fintype.card (DiskDegreeIndex N))) ℝ (restrictTotalDegree (Fin 2) ℝ N) :=
  (normalizedProriolDegreeBasis hα hβ hγ N).reindex
    ((diskDegreeIndexEquiv N).symm.trans (diskDegreeOrder N).symm.toEquiv)

/-- Both support bounds apply to the actual canonical Cholesky factor. -/
theorem integralGramFactor_proriol_eq_zero {α β γ : ℝ}
    (hα : -1 < α) (hβ : -1 < β) (hγ : -1 < γ) (a b c N : ℕ)
    (hαa : -1 < α+a) (hβb : -1 < β+b) (hγc : -1 < γ+c)
    (i j : Fin (Fintype.card (DiskDegreeIndex N)))
    (hgap : ((diskDegreeOrder N j).1 1 < (diskDegreeOrder N i).1 1 ∨
      (diskDegreeOrder N i).1 1+(b+c) < (diskDegreeOrder N j).1 1) ∨
      ((diskDegreeOrder N j).1 0+(diskDegreeOrder N j).1 1 <
        (diskDegreeOrder N i).1 0+(diskDegreeOrder N i).1 1 ∨
        (diskDegreeOrder N i).1 0+(diskDegreeOrder N i).1 1+(a+b+c) <
          (diskDegreeOrder N j).1 0+(diskDegreeOrder N j).1 1)) :
    integralGramFactor (proriolMeasure (α+a) (β+b) (γ+c))
      (hasFiniteMoments_proriolMeasure hαa hβb hγc)
      (hasPositiveDefiniteIntegral_proriolMeasure hαa hβb hγc)
      (orderedProriolBasis hα hβ hγ N) i j = 0 := by
  rw [integralGramFactor_proriol_eq_connection]
  simp only [orderedProriolBasis, Basis.toMatrix_apply, Basis.reindex_apply,
    Basis.repr_reindex_apply]
  rw [normalizedProriolDegreeBasis_connection]
  apply mul_eq_zero_of_right
  have he := proriol_connection_eq_zero hα hβ hγ a b c N hαa hβb hγc
    ((diskDegreeIndexEquiv N) (diskDegreeOrder N j))
    ((diskDegreeIndexEquiv N) (diskDegreeOrder N i)) hgap
  simpa using he

/-- The finite normalized connection for a unit increase of the disk parameter. -/
noncomputable def dunklXuRaisingMatrix {β : ℝ} (hβ : -1 < β) (N : ℕ) :
    Matrix {d : Fin 2 →₀ ℕ // d.sum (fun _ e => e) ≤ N}
      {d : Fin 2 →₀ ℕ // d.sum (fun _ e => e) ≤ N} ℝ :=
  (normalizedDunklXuDegreeBasis (show -1 < β + 1 by linarith) N).toMatrix
    (normalizedDunklXuDegreeBasis hβ N)

/-- The positive normalized factor has the four-term disk support. -/
theorem dunklXuRaisingMatrix_eq_zero {β : ℝ} (hβ : -1 < β) (N : ℕ)
    (i j : {d : Fin 2 →₀ ℕ // d.sum (fun _ e => e) ≤ N})
    (hout : ¬ ((i.1 0 + i.1 1 = j.1 0 + j.1 1 ∨
      i.1 0 + i.1 1 + 2 = j.1 0 + j.1 1) ∧
      (i.1 1 = j.1 1 ∨ i.1 1 + 2 = j.1 1))) :
    dunklXuRaisingMatrix hβ N i j = 0 := by
  rw [dunklXuRaisingMatrix_apply,
    dunklXuBasis_raise_repr_eq_zero hβ
      ⟨i.1 0 + i.1 1, ⟨i.1 1, by omega⟩⟩
      ⟨j.1 0 + j.1 1, ⟨j.1 1, by omega⟩⟩ hout, mul_zero]

/-- Every diagonal entry of the normalized disk factor is strictly positive. -/
theorem dunklXuRaisingMatrix_diag_pos {β : ℝ} (hβ : -1 < β) (N : ℕ)
    (i : {d : Fin 2 →₀ ℕ // d.sum (fun _ e => e) ≤ N}) :
    0 < dunklXuRaisingMatrix hβ N i i := by
  rw [dunklXuRaisingMatrix_apply]
  exact mul_pos (mul_pos (dunklXuNorm_pos (by linarith) _ _)
    (inv_pos.mpr (dunklXuNorm_pos hβ _ _))) (dunklXuBasis_raise_repr_self_pos hβ ⟨i.1 0 + i.1 1, ⟨i.1 1, by omega⟩⟩)

/-- The factor is upper triangular in degree order, with increasing inner index
within each degree block. -/
theorem dunklXuRaisingMatrix_eq_zero_of_order {β : ℝ} (hβ : -1 < β) (N : ℕ)
    (i j : {d : Fin 2 →₀ ℕ // d.sum (fun _ e => e) ≤ N})
    (hij : j.1 0 + j.1 1 < i.1 0 + i.1 1 ∨
      (j.1 0 + j.1 1 = i.1 0 + i.1 1 ∧ j.1 1 < i.1 1)) :
    dunklXuRaisingMatrix hβ N i j = 0 := by
  apply dunklXuRaisingMatrix_eq_zero hβ N i j
  omega

/-- A polynomial outer modification admits an expansion using only the last
degree(u)+1 outer indices and the same fixed inner lift. -/
theorem IsOrthogonal.separated_modification_expansion {A : Type*} [CommRing A] [Algebra ℝ A]
    {u : ℝ[X]} {μ : Measure {x : ℝ // 0 ≤ u.eval x}}
    (hμ : HasFiniteMoments (μ.map Subtype.val))
    (hpos : HasPositiveDefiniteIntegral (μ.map Subtype.val)) (hu0 : u ≠ 0)
    {P Q : Sequence ℝ} (hP : P.IsOrthogonal (μ.map Subtype.val))
    (hQ : Q.IsOrthogonal (polynomialDensityMeasure u μ))
    (hunits : ∀ i, IsUnit (Q i).leadingCoeff)
    (φ : ℝ[X] →ₐ[ℝ] A) (v : A) (n : ℕ) :
    φ (P n) * v = ∑ m ∈ (Finset.range (n+1)).filter (fun m => n ≤ m+u.natDegree),
      (Q.basis hunits).repr (P n) m • (φ (Q m) * v) := by
  classical
  rw [hQ.separated_expansion (HasFiniteMoments.withDensity_polynomial hμ)
    (HasPositiveDefiniteIntegral.withDensity_polynomial hμ hpos hu0) hunits φ v P n]
  symm
  apply Finset.sum_subset (Finset.filter_subset _ _)
  intro m hm hn
  have hout : m+u.natDegree < n := by
    simp only [Finset.mem_filter, hm, true_and, not_le] at hn
    exact hn
  rw [IsOrthogonal.basis_repr_eq_zero_of_polynomial_modification hμ hpos hu0 hP hQ hunits (Or.inr hout), zero_smul]