/-- The unique positive upper triangular factor for an ordered polynomial basis. -/
noncomputable def integralGramFactor : Matrix (Fin n) (Fin n) ℝ :=
  (existsUnique_pos_upperTriangular_integralGram_factor μ hμ hpos b).choose

/-- Any concrete positive upper triangular factor identifies the canonical factor. -/
theorem integralGramFactor_eq (R : Matrix (Fin n) (Fin n) ℝ)
    (hR : R.IsUpperTriangular) (hp : ∀ i, 0 < R i i)
    (hg : integralGram μ (fun i => (b i : MvPolynomial σ ℝ)) = R.transpose*R) :
    integralGramFactor μ hμ hpos b = R :=
  ((existsUnique_pos_upperTriangular_integralGram_factor μ hμ hpos b).choose_spec.2
    R ⟨hR,hp,hg⟩).symm

/-- Support of the canonical factor can be verified by any positive triangular Gram factor. -/
theorem integralGramFactor_hasOffsetSupport_iff {G : Type*} [AddCommGroup G]
    (label : Fin n → G) (S : Finset G) :
    (integralGramFactor μ hμ hpos b).HasOffsetSupport label label S ↔
      ∃ R : Matrix (Fin n) (Fin n) ℝ, R.IsUpperTriangular ∧ (∀ i, 0 < R i i) ∧
        integralGram μ (fun i => (b i : MvPolynomial σ ℝ)) = R.transpose*R ∧
        R.HasOffsetSupport label label S := by
  constructor
  · intro h
    have hs := integralGramFactor_spec μ hμ hpos b
    exact ⟨_,hs.1,hs.2.1,hs.2.2,h⟩
  · rintro ⟨R,hR,hp,hg,hs⟩
    rw [integralGramFactor_eq μ hμ hpos b R hR hp hg]
    exact hs

/-- A function of two variables is a product of functions of the individual variables. -/
def IsProductSeparated (f : α → β → K) : Prop :=
  ∃ a : α → K, ∃ b : β → K, ∀ x y, f x y = a x * b y

/-- A zero selected upper-factor row forces the corresponding bordered minor to vanish. -/
theorem det_submatrix_eq_zero_of_lower_factorization
    (W L U : Matrix (Fin (n+m)) (Fin (n+m)) ℝ)
    (hW : W = L*U) (hL : L.transpose.IsUpperTriangular)
    (col : Fin n → Fin (n+m)) (i : Fin n)
    (hz : ∀ j, U (Fin.castAdd m i) (col j) = 0) :
    (W.submatrix (Fin.castAdd m) col).det = 0 := by
  rw [hW, lowerTriangular_mul_submatrix L U hL col, det_mul,
    det_eq_zero_of_row_eq_zero (A := U.submatrix (Fin.castAdd m) col) i hz, mul_zero]

/-- Degree two has only the top and middle terms of the double conversion. -/
theorem gegenbauer_double_conversion_two (lam : R) :
    C ((lam+2)*(lam+3))*gegenbauer lam 2 =
      C (lam*(lam+1))*gegenbauer (lam+2) 2 - C (2*lam*(lam+2)) := by
  have h₀ := gegenbauer_conversion lam 0
  have h₁ := gegenbauer_conversion (lam+1) 0
  have he : lam+1+1=lam+2 := by ring
  simp only [Nat.cast_zero, zero_add, he, gegenbauer_zero,
    map_mul, map_add, map_ofNat, map_one] at h₀ h₁ ⊢
  linear_combination (C lam+3)*h₀+C lam*h₁

/-- Degree three has only the top and middle terms of the double conversion. -/
theorem gegenbauer_double_conversion_three (lam : R) :
    C ((lam+3)*(lam+4)*(lam+2))*gegenbauer lam 3 =
      C (lam*(lam+1)*(lam+2))*gegenbauer (lam+2) 3 -
      C (2*lam*(lam+1)*(lam+3))*gegenbauer (lam+2) 1 := by
  have h₀ := gegenbauer_conversion lam 1
  have h₁ := gegenbauer_conversion (lam+1) 1
  have h₂ := gegenbauer_conversion_one (lam+1)
  have he : lam+1+1=lam+2 := by ring
  simp only [Nat.cast_one, he, map_mul, map_add, map_ofNat, map_one] at h₀ h₁ h₂ ⊢
  linear_combination
    (C lam+4)*(C lam+2)*h₀+C lam*(C lam+2)*h₁-C lam*(C lam+4)*h₂

/-- Retain only entries that raise the column level by exactly one. -/
def maximalStep (A : Matrix ι ι ℝ) (level : ι → ℤ) : Matrix ι ι ℝ :=
  fun i j => if level j = level i+1 then A i j else 0

/-- At the maximal possible level gap, only maximal-step entries contribute. -/
theorem pow_apply_eq_maximalStep (A : Matrix ι ι ℝ) (level : ι → ℤ)
    (hA : ∀ i j, level i+1 < level j → A i j = 0)
    (n : ℕ) (i j : ι) (hgap : level j = level i+n) :
    (A^n) i j = ((maximalStep A level)^n) i j := by
  induction n generalizing i j with
  | zero => rfl
  | succ n ih =>
    rw [pow_succ', pow_succ', mul_apply, mul_apply]
    apply Finset.sum_congr rfl
    intro k _
    by_cases hk : level k = level i+1
    · rw [maximalStep, ite_eq_left hk, ih k j (by omega)]
    · rw [maximalStep, ite_eq_right hk, zero_mul]
      by_cases hbig : level i+1 < level k
      · rw [hA i k hbig, zero_mul]
      · rw [pow_apply_eq_zero_of_level_gap A level hA n k j (by omega), mul_zero]

/-- One positive path suffices for a strictly positive entry of a nonnegative matrix power. -/
theorem pow_apply_pos_of_positive_path (A : Matrix ι ι ℝ)
    (hA : ∀ i j, 0 ≤ A i j) (n : ℕ) (p : ℕ → ι)
    (hp : ∀ k < n, 0 < A (p k) (p (k+1))) : 0 < (A^n) (p 0) (p n) := by
  induction n generalizing p with
  | zero => simp
  | succ n ih =>
    rw [pow_succ', mul_apply]
    have ht := ih (fun k => p (k+1)) (fun k hk => hp (k+1) (by omega))
    have hv : 0 < A (p 0) (p 1) * (A^n) (p 1) (p (n+1)) :=
      mul_pos (hp 0 (by omega)) ht
    apply hv.trans_le
    exact Finset.single_le_sum
      (fun k _ => mul_nonneg (hA (p 0) k) (pow_apply_nonneg hA n k (p (n+1))))
      (Finset.mem_univ (p 1))

/-- Diagonal sign changes telescope along every matrix-power path. -/
theorem pow_apply_sign_conjugate (A : Matrix ι ι ℝ) (s : ι → ℝ)
    (hs : ∀ i, s i*s i=1) (n : ℕ) (i j : ι) :
    ((Matrix.of (fun i j => s i*A i j*s j))^n) i j = s i*(A^n) i j*s j := by
  induction n generalizing i j with
  | zero =>
    by_cases h : i=j
    · subst j; simp [hs]
    · simp [h]
  | succ n ih =>
    rw [pow_succ', pow_succ', mul_apply, mul_apply, Finset.mul_sum, Finset.sum_mul]
    apply Finset.sum_congr rfl
    intro k _
    rw [ih]
    change (s i*A i k*s k)*(s k*(A^n) k j*s j) = s i*(A i k*(A^n) k j)*s j
    calc
      _ = (s i*(A i k*(A^n) k j)*s j)*(s k*s k) := by ring
      _ = _ := by rw [hs, mul_one]

/-- Exact factorization after all lower inverse-path coefficients vanish. -/
theorem apply_rightInverse_eq_of_powers_vanish (f : R →ₗ[ℝ] ℝ)
    (K Y : R) (t : ℝ) (hy : (1+t • K)*Y=1) (r : ℕ)
    (hz : ∀ q < r, f (K^q)=0) :
    f Y = (-t)^r * f (K^r*Y) := by
  have he := eq_geom_sum_add_pow_mul_of_one_sub_mul_eq_one (-t • K) Y
    (by simpa using hy) r
  conv_lhs => rw [he, map_add, map_sum]
  have hs : (∑ q ∈ Finset.range r, f ((-t • K)^q))=0 := by
    apply Finset.sum_eq_zero
    intro q hq
    rw [smul_pow, map_smul, hz q (Finset.mem_range.mp hq), smul_zero]
  rw [hs, zero_add, smul_pow, smul_mul_assoc, map_smul, smul_eq_mul]

/-- A nonzero surviving inverse-path remainder gives a nonzero pairing. -/
theorem apply_rightInverse_ne_zero_of_powers_vanish (f : R →ₗ[ℝ] ℝ)
    (K Y : R) (t : ℝ) (ht : t ≠ 0) (hy : (1+t • K)*Y=1) (r : ℕ)
    (hz : ∀ q < r, f (K^q)=0) (hr : f (K^r*Y) ≠ 0) : f Y ≠ 0 := by
  rw [f.apply_rightInverse_eq_of_powers_vanish K Y t hy r hz]
  exact mul_ne_zero (pow_ne_zero r (neg_ne_zero.mpr ht)) hr

/-- The first Gram row of an upper triangular matrix is its first row times the first pivot. -/
theorem IsUpperTriangular.transpose_mul_self_first_row {n : ℕ}
    {D : Matrix (Fin (n+1)) (Fin (n+1)) ℝ} (hD : D.IsUpperTriangular) (j : Fin (n+1)) :
    (D.transpose*D) 0 j = D 0 0 * D 0 j := by
  rw [mul_apply]
  apply Finset.sum_eq_single 0
  · intro k _ hk
    change D k 0 * D k j = 0
    rw [hD (Fin.pos_iff_ne_zero.mpr hk), zero_mul]
  · simp

/-- For affine leading and off-diagonal blocks, a zero original entry leaves the inverse sandwich. -/
theorem schur_gram_first_row_eq_of_pencil {n : ℕ}
    (A : Matrix ι ι ℝ) (B : Matrix ι (Fin (n+1)) ℝ)
    (D : Matrix (Fin (n+1)) (Fin (n+1)) ℝ)
    (K : Matrix ι ι ℝ) (E : Matrix ι (Fin (n+1)) ℝ)
    (F : Matrix (Fin (n+1)) ι ℝ) (t : ℝ)
    (hA : IsUnit A.det) (hD : D.IsUpperTriangular)
    (hAA : A.transpose*A = 1+t • K) (hAB : A.transpose*B = t • E)
    (hBA : B.transpose*A = t • F) (j : Fin (n+1))
    (hij : (B.transpose*B+D.transpose*D) 0 j = 0) :
    D 0 0 * D 0 j = -t^2 * (F*(1+t • K)⁻¹*E) 0 j := by
  have he := congrArg (fun M : Matrix (Fin (n+1)) (Fin (n+1)) ℝ => M 0 j)
    (schur_gram_upper_blocks A B D hA)
  rw [hAA, hAB, hBA, sub_apply, hij, hD.transpose_mul_self_first_row] at he
  simp only [Matrix.smul_mul, Matrix.mul_smul, smul_apply, smul_eq_mul] at he
  nlinarith [he]

/-- A selected entry of rectangular left and right multiplication, as a linear functional. -/
noncomputable def sandwichEntry (F : Matrix (Fin (n+1)) ι ℝ) (E : Matrix ι (Fin (n+1)) ℝ)
    (i j : Fin (n+1)) : Matrix ι ι ℝ →ₗ[ℝ] ℝ where
  toFun M := (F*M*E) i j
  map_add' M N := by simp [Matrix.mul_add, Matrix.add_mul]
  map_smul' c M := by simp [Matrix.mul_smul, Matrix.smul_mul]

/-- Vanishing lower powers and one surviving path pairing force local factor nonvanishing. -/
theorem eventually_cholesky_first_row_ne_zero
    (K : Matrix ι ι ℝ) (E : Matrix ι (Fin (n+1)) ℝ)
    (F : Matrix (Fin (n+1)) ι ℝ) (j : Fin (n+1)) (r : ℕ)
    (hz : ∀ q < r, (F*K^q*E) 0 j = 0) (hr : (F*K^r*E) 0 j ≠ 0) :
    ∀ᶠ t in 𝓝[≠] (0 : ℝ),
      ∀ (A : Matrix ι ι ℝ) (B : Matrix ι (Fin (n+1)) ℝ)
        (D : Matrix (Fin (n+1)) (Fin (n+1)) ℝ),
        IsUnit A.det → D.IsUpperTriangular →
        A.transpose*A = 1+t • K → A.transpose*B = t • E → B.transpose*A = t • F →
        (B.transpose*B+D.transpose*D) 0 j = 0 → D 0 j ≠ 0 := by
  have hh := eventually_apply_inv_one_add_smul_ne_zero (sandwichEntry F E 0 j) K r hz hr
  filter_upwards [hh, self_mem_nhdsWithin] with t ht ht0
  intro A B D hA hD hAA hAB hBA hij
  exact schur_gram_first_row_ne_zero_of_pencil A B D K E F t (by simpa using ht0)
    hA hD hAA hAB hBA j hij ht

/-- Squaring recovers the affine squared coordinate everywhere in the domain. -/
theorem diskDeformation_sq (δ : ℝ) (x : {x : ℝ // 0 ≤ 1-δ+δ*x^2}) :
    diskDeformation δ x ^ 2 = 1-δ+δ*(x : ℝ)^2 := Real.sq_sqrt x.property

/-- Parameters in the closed unit interval admit every real source coordinate. -/
theorem diskDeformation_radical_nonneg {δ : ℝ} (hδ : 0 ≤ δ) (hδ1 : δ ≤ 1) (x : ℝ) :
    0 ≤ 1-δ+δ*x^2 := add_nonneg (sub_nonneg.mpr hδ1) (mul_nonneg hδ (sq_nonneg x))

/-- The branch maps unit-interval source points into the unit interval. -/
theorem diskDeformation_mem_Icc {δ : ℝ} (hδ : 0 ≤ δ) (_hδ1 : δ ≤ 1)
    (x : {x : ℝ // 0 ≤ 1-δ+δ*x^2}) (hx : (x : ℝ) ∈ Icc (-1 : ℝ) 1) :
    diskDeformation δ x ∈ Icc (-1 : ℝ) 1 := by
  have hn : 0 ≤ diskDeformation δ x := Real.sqrt_nonneg _
  have he := diskDeformation_sq δ x
  have hx2 : (x : ℝ)^2 ≤ 1 := by nlinarith [hx.1,hx.2]
  constructor <;> nlinarith [mul_le_mul_of_nonneg_left hx2 hδ]

/-- Even test functions have a single classical branch integral. -/
theorem integral_diskDeformedMeasure_even (δ : ℝ)
    {μ : Measure {x : ℝ // 0 ≤ 1-δ+δ*x^2}} [IsFiniteMeasure μ]
    (hx : ∀ᵐ x : {x : ℝ // 0 ≤ 1-δ+δ*x^2} ∂μ, (x : ℝ) ∈ Icc (-1 : ℝ) 1)
    {f : ℝ → ℝ} (hf : Continuous f) (he : Function.Even f) :
    ∫ x, f x ∂diskDeformedMeasure δ μ = ∫ t, f (diskDeformation δ t) ∂μ := by
  rw [integral_diskDeformedMeasure δ hx hf]
  have heq : (fun t => f (-diskDeformation δ t)) = (fun t => f (diskDeformation δ t)) := _root_.funext (fun t => he (diskDeformation δ t))
  rw [heq]
  ring

/-- Nonzero measure mass survives the symmetric pushforward. -/
theorem diskDeformedMeasure_ne_zero (δ : ℝ)
    {μ : Measure {x : ℝ // 0 ≤ 1-δ+δ*x^2}} (hμ : μ ≠ 0) :
    diskDeformedMeasure δ μ ≠ 0 := by
  intro h
  have hmap : μ.map (diskDeformation δ) = 0 := by
    have hh : μ.map (diskDeformation δ) + μ.map (fun x => -diskDeformation δ x) = 0 :=
      (smul_eq_zero.mp h).resolve_left (by norm_num)
    apply le_antisymm _ bot_le
    calc
      _ ≤ μ.map (diskDeformation δ) + μ.map (fun x => -diskDeformation δ x) :=
        Measure.le_add_right le_rfl
      _ = 0 := hh
  exact hμ ((Measure.map_eq_zero_iff (continuous_diskDeformation δ).measurable.aemeasurable).mp hmap)

/-- Every nonzero deformation parameter has finite fibers on its classical domain. -/
theorem finite_fiber_diskDeformation (δ : ℝ) (hδ : δ ≠ 0) (y : ℝ) :
    {x | diskDeformation δ x = y}.Finite := by
  let p : ℝ[X] := C δ * X^2 + C (1-δ-y^2)
  have hp : p ≠ 0 := by
    intro h
    have he := congrArg (fun q : ℝ[X] => q.coeff 2) h
    have hc : p.coeff 2 = δ := by
      change (C δ * X^2 + C (1-δ-y^2) : ℝ[X]).coeff 2 = δ
      simp only [coeff_add, coeff_C_mul, coeff_X_pow, coeff_C]
      norm_num
    rw [hc, coeff_zero] at he
    exact hδ he
  have hf := Set.Finite.preimage (f := (Subtype.val : {x : ℝ // 0 ≤ 1-δ+δ*x^2} → ℝ))
    Subtype.val_injective.injOn (finite_setOfPred_isRoot hp)
  apply hf.subset
  intro x hx
  change p.eval (x : ℝ) = 0
  have he := diskDeformation_sq δ x
  rw [hx] at he
  simp only [p, eval_add, eval_mul, eval_C, eval_pow, eval_X]
  linarith

/-- Finite supported real source measures retain all polynomial moments. -/
@[aesop (rule_sets := [PolynomialIntegrability]) unsafe 80% apply]
theorem hasFiniteMoments_diskDeformedMeasure {μ : Measure ℝ} [IsFiniteMeasure μ]
    {δ : ℝ} (hδ : 0 ≤ δ) (hδ1 : δ ≤ 1) (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1) :
    HasFiniteMoments (diskDeformedMeasure δ (μ.comap Subtype.val)) :=
  hasFiniteMoments_diskDeformedMeasure_on_domain hδ hδ1 (ae_mem_Icc_comap_diskDeformation δ hx)

/-- Atomlessness survives supported real-measure transport and nondegenerate deformation. -/
theorem nullSingletonClass_diskDeformedMeasure {μ : Measure ℝ} [NullSingletonClass μ]
    {δ : ℝ} (hδ : 0 < δ) (_hδ1 : δ ≤ 1) :
    NullSingletonClass (diskDeformedMeasure δ (μ.comap Subtype.val)) :=
  nullSingletonClass_diskDeformedMeasure_on_domain hδ.ne' _

/-- A nonzero finite atomless real source gives a positive deformed polynomial pairing. -/
@[aesop (rule_sets := [PolynomialIntegrability]) unsafe 80% apply]
theorem hasPositiveDefiniteIntegral_diskDeformedMeasure {μ : Measure ℝ}
    [IsFiniteMeasure μ] [NullSingletonClass μ] (hμ : μ ≠ 0)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1) (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1) :
    HasPositiveDefiniteIntegral (diskDeformedMeasure δ (μ.comap Subtype.val)) :=
  hasPositiveDefiniteIntegral_diskDeformedMeasure_on_domain
    (comap_diskDeformation_ne_zero hδ.le hδ1 hμ) hδ hδ1 (ae_mem_Icc_comap_diskDeformation δ hx)

/-- The degree-shifted outer integral acquires exactly the factor `δ^k`. -/
theorem integral_diskOuterMeasure_diskDeformedMeasure {μ : Measure ℝ} [IsFiniteMeasure μ]
    (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    {δ : ℝ} (hδ : 0 ≤ δ) (hδ1 : δ ≤ 1) (k : ℕ)
    {f : ℝ → ℝ} (hf : Continuous f) :
    ∫ x, f x ∂diskOuterMeasure k ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val) =
      δ^k * ((1/2 : ℝ) * (∫ t, f (diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ hδ1 t⟩) ∂diskOuterMeasure k (μ.comap Subtype.val)) +
        (1/2 : ℝ) * (∫ t, f (-diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ hδ1 t⟩) ∂diskOuterMeasure k (μ.comap Subtype.val))) := by
  have hc : Continuous (fun t : ℝ =>
      diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ hδ1 t⟩) :=
    (continuous_diskDeformation δ).comp (continuous_id.subtype_mk _)
  have hx' := ae_mem_Icc_diskDeformedMeasure hδ hδ1 hx
  erw [integral_diskOuterMeasure _ hx' k
      ⟨f, integrable_diskOuterMeasure_of_continuous hx' k hf⟩,
    integral_diskDeformedMeasure_comap hx hδ hδ1 (by
      change Continuous (fun x : ℝ => (1-x^2)^k*f x)
      fun_prop),
    integral_diskOuterMeasure μ hx k
      ⟨fun t => f (diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ hδ1 t⟩),
        integrable_diskOuterMeasure_of_continuous hx k (hf.comp hc)⟩,
    integral_diskOuterMeasure μ hx k
      ⟨fun t => f (-diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ hδ1 t⟩),
        integrable_diskOuterMeasure_of_continuous hx k (hf.comp hc.neg)⟩]
  dsimp only
  simp_rw [neg_sq, one_sub_diskDeformation_sq hδ hδ1, mul_pow, mul_assoc,
    integral_const_mul]
  ring

/-- The transported polynomial is even. -/
theorem even_eval_evenDiskTransport (δ : {δ : ℝ // δ ≠ 0}) (p : ℝ[X]) :
    Function.Even (evenDiskTransport δ p).eval := by
  intro x
  simp [evenDiskTransport, eval_comp]

/-- The transported polynomial is explicitly a polynomial in `X²`. -/
theorem evenDiskTransport_eq_comp_sq (δ : {δ : ℝ // δ ≠ 0}) (p : ℝ[X]) :
    evenDiskTransport δ p =
      ((p.contract 2).comp (C (δ : ℝ)⁻¹ * (X-C (1-(δ : ℝ))))).comp (X^2) := by
  rw [evenDiskTransport, comp_assoc]
  congr 1
  simp only [mul_comp, C_comp, sub_comp, X_comp]

/-- The inverse substitution recovers an even polynomial on either branch. -/
theorem eval_evenDiskTransport_sq (δ : ℝ) (hδ : δ ≠ 0) (p : ℝ[X])
    (hp : ∀ i, i % 2 ≠ 0 → p.coeff i = 0) (x t : ℝ)
    (hx : x^2 = 1-δ+δ*t^2) :
    (evenDiskTransport ⟨δ,hδ⟩ p).eval x = p.eval t := by
  have he := congrArg (fun q : ℝ[X] => q.eval t) (comp_sq_contract_of_even_coeff p hp)
  simp only [eval_comp, eval_pow, eval_X] at he
  rw [evenDiskTransport, eval_comp]
  have ha : (C δ⁻¹ * (X^2-C (1-δ)) : ℝ[X]).eval x = t^2 := by
    simp only [eval_mul, eval_C, eval_sub, eval_pow, eval_X, hx]
    field_simp
    ring
  rw [ha, he]

/-- Transport preserves the degree of every even polynomial. -/
theorem natDegree_evenDiskTransport (δ : ℝ) (hδ : δ ≠ 0) (p : ℝ[X])
    (hp : ∀ i, i % 2 ≠ 0 → p.coeff i = 0) :
    (evenDiskTransport ⟨δ,hδ⟩ p).natDegree = p.natDegree := by
  have he := natDegree_evenDiskTransport_comp_sq δ hδ (p.contract 2)
  rwa [comp_sq_contract_of_even_coeff p hp] at he

/-- A positive parameter preserves positivity of the leading coefficient. -/
theorem leadingCoeff_evenDiskTransport_pos (δ : ℝ) (hδ : 0 < δ) (p : ℝ[X])
    (hp : ∀ i, i % 2 ≠ 0 → p.coeff i = 0) (hl : 0 < p.leadingCoeff) :
    0 < (evenDiskTransport ⟨δ,hδ.ne'⟩ p).leadingCoeff := by
  have he := leadingCoeff_evenDiskTransport_comp_sq δ hδ.ne' (p.contract 2)
  rw [comp_sq_contract_of_even_coeff p hp] at he
  rw [he]
  exact mul_pos hl (pow_pos (inv_pos.mpr hδ) _)

/-- Transport by the reciprocal parameter is the inverse on even polynomials. -/
theorem evenDiskTransport_inv (δ : ℝ) (hδ : δ ≠ 0) (p : ℝ[X])
    (hp : ∀ i, i % 2 ≠ 0 → p.coeff i = 0) :
    evenDiskTransport ⟨δ⁻¹,inv_ne_zero hδ⟩ (evenDiskTransport ⟨δ,hδ⟩ p) = p := by
  rw [evenDiskTransport_eq_comp_sq ⟨δ,hδ⟩ p, evenDiskTransport_comp_sq, comp_assoc]
  have he : (C δ⁻¹ * (X-C (1-δ))).comp
      (C (δ⁻¹)⁻¹ * (X^2-C (1-δ⁻¹))) = (X^2 : ℝ[X]) := by
    simp only [mul_comp, C_comp, sub_comp, X_comp, inv_inv]
    have hc : C δ⁻¹ * C δ = (1 : ℝ[X]) := by rw [← C_mul, inv_mul_cancel₀ hδ, C_1]
    simp only [map_sub, map_one]
    linear_combination (X^2+C δ⁻¹)*hc
  rw [he, comp_sq_contract_of_even_coeff p hp]

/-- The factor `δ⁻¹^b` makes transport an isometry at inner degree `2*b`. -/
theorem integralInner_normalized_evenDiskTransport {μ : Measure ℝ} [IsFiniteMeasure μ]
    (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1) (b : ℕ) (p q : ℝ[X])
    (hp : ∀ i, i % 2 ≠ 0 → p.coeff i = 0)
    (hq : ∀ i, i % 2 ≠ 0 → q.coeff i = 0) :
    integralInner (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)) (IntegrablePair.of (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)) (δ⁻¹^b • evenDiskTransport ⟨δ,hδ.ne'⟩ p) (δ⁻¹^b • evenDiskTransport ⟨δ,hδ.ne'⟩ q) (h := ((hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 hx).diskOuterMeasure (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 hx) _).integrable_mul _ _)) =
      integralInner (diskOuterMeasure (2*b) (μ.comap Subtype.val)) (IntegrablePair.of (diskOuterMeasure (2*b) (μ.comap Subtype.val)) p q (h := ((hasFiniteMoments_of_ae_mem_Icc hx).diskOuterMeasure hx _).integrable_mul _ _)) := by
  have hμδ := (hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 hx).diskOuterMeasure
    (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 hx) (2*b)
  rw [integralInner_smul_left (h := hμδ.integrable_mul _ _), integralInner_comm _ (evenDiskTransport ⟨δ,hδ.ne'⟩ p),
    integralInner_smul_left (h := hμδ.integrable_mul _ _), integralInner_comm _ (evenDiskTransport ⟨δ,hδ.ne'⟩ q),
    integralInner_evenDiskTransport hx hδ hδ1 (2*b) p q hp hq]
  have he : δ⁻¹^b * (δ⁻¹^b * δ^(2*b)) = 1 := by
    rw [pow_mul]
    calc
      _ = (δ⁻¹ * δ)^(2*b) := by rw [mul_pow, pow_mul]; ring
      _ = 1 := by rw [inv_mul_cancel₀ hδ.ne', one_pow]
  calc
    _ = (δ⁻¹^b * (δ⁻¹^b * δ^(2*b))) * integralInner (diskOuterMeasure (2*b) (μ.comap Subtype.val)) (IntegrablePair.of (diskOuterMeasure (2*b) (μ.comap Subtype.val)) p q (h := ((hasFiniteMoments_of_ae_mem_Icc hx).diskOuterMeasure hx _).integrable_mul _ _)) := by ring
    _ = _ := by rw [he, one_mul]

/-- Positive-leading unit normalization makes the even transport equal the target polynomial. -/
theorem normalized_evenDiskTransport_eq
    {μ : Measure ℝ} [IsFiniteMeasure μ]
    (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1) (hμ : HasFiniteMoments μ)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1) (b n : ℕ)
    (hpos : HasPositiveDefiniteIntegral (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)))
    {P Q : Sequence ℝ} (hP : P.IsOrthogonal (diskOuterMeasure (2*b) (μ.comap Subtype.val)))
    (hQ : Q.IsOrthogonal (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)))
    (hp : ∀ i, i % 2 ≠ 0 → (P n).coeff i = 0)
    (hq : ∀ i, i % 2 ≠ 0 → (Q n).coeff i = 0)
    (hpl : 0 < (P n).leadingCoeff) (hql : 0 < (Q n).leadingCoeff)
    (hpn : integralInner (diskOuterMeasure (2*b) (μ.comap Subtype.val)) (IntegrablePair.of (diskOuterMeasure (2*b) (μ.comap Subtype.val)) (P n) (P n) (h := ((hasFiniteMoments_of_ae_mem_Icc hx).diskOuterMeasure hx _).integrable_mul _ _)) = 1)
    (hqn : integralInner (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)) (IntegrablePair.of (diskOuterMeasure (2*b) ((diskDeformedMeasure δ (μ.comap Subtype.val)).comap Subtype.val)) (Q n) (Q n) (h := ((hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 hx).diskOuterMeasure (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 hx) _).integrable_mul _ _)) = 1) :
    δ⁻¹^b • evenDiskTransport ⟨δ,hδ.ne'⟩ (P n) = Q n := by
  let c := (evenDiskTransport ⟨δ,hδ.ne'⟩ (P n)).leadingCoeff / (Q n).leadingCoeff
  have hc : 0 < c := div_pos (leadingCoeff_evenDiskTransport_pos δ hδ (P n) hp hpl) hql
  have he : evenDiskTransport ⟨δ,hδ.ne'⟩ (P n) = c • Q n :=
    evenDiskTransport_eq_leadingCoeff_div_smul hx hμ hδ hδ1 (2*b) n hpos hP hQ hp hq hpl
  have hn := integralInner_evenDiskTransport hx hδ hδ1 (2*b) (P n) (P n) hp hp
  have hμδ := (hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 hx).diskOuterMeasure
    (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 hx) (2*b)
  have heq := congrArg (fun q => integralInner _ (hμδ.pair q q)) he
  rw [heq, integralInner_smul_left (h := hμδ.integrable_mul _ _), integralInner_comm _ (Q n),
    integralInner_smul_left (h := hμδ.integrable_mul _ _), hqn, hpn] at hn
  have hec : c = δ^b := by
    have hpow : δ^(2*b) = (δ^b)^2 := by rw [Nat.mul_comm, pow_mul]
    rw [hpow] at hn
    nlinarith [pow_pos hδ b]
  rw [he, hec, smul_smul, ← mul_pow, inv_mul_cancel₀ hδ.ne', one_pow, one_smul]

/-- Even outer polynomials in the two systems agree after normalized transport. -/
theorem outer_even_diskDeformation_eq
    {μ ν : Measure ℝ} [IsFiniteMeasure μ] [μ.IsNegInvariant]
    (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    (hμ : Polynomial.HasFiniteMoments μ) (hμpos : Polynomial.HasPositiveDefiniteIntegral μ)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1)
    (hδpos : Polynomial.HasPositiveDefiniteIntegral (diskDeformedMeasure δ (μ.comap Subtype.val)))
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ (μ.comap Subtype.val))) ν)
    (a b : ℕ) :
    δ⁻¹^b • evenDiskTransport ⟨δ,hδ.ne'⟩ (S.outer (2*b) (2*a)) = T.outer (2*b) (2*a) := by
  let := isNegInvariant_diskOuterMeasure μ (2*b)
  let := isNegInvariant_diskOuterMeasure (diskDeformedMeasure δ (μ.comap Subtype.val)) (2*b)
  have hδmom := hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 hx
  have hδsupp := ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 hx
  apply normalized_evenDiskTransport_eq hx hμ hδ hδ1 b (2*a)
    (hδpos.diskOuterMeasure hδmom hδsupp (2*b))
    (S.outer_orthogonal_real (2*b)) (T.outer_orthogonal_real (2*b))
  · intro i hi
    exact (S.outer_orthogonal_real (2*b)).coeff_eq_zero_of_parity
      (hμ.diskOuterMeasure hx (2*b)) (hμpos.diskOuterMeasure hμ hx (2*b))
      (2*a) i (by simpa using hi)
  · intro i hi
    exact (T.outer_orthogonal_real (2*b)).coeff_eq_zero_of_parity
      (hδmom.diskOuterMeasure hδsupp (2*b)) (hδpos.diskOuterMeasure hδmom hδsupp (2*b))
      (2*a) i (by simpa using hi)
  · exact S.outer_leadingCoeff_pos _ _
  · exact T.outer_leadingCoeff_pos _ _
  · simpa only [Polynomial.integralInner_self] using S.outer_norm_real (2*b) (2*a)
  · simpa only [Polynomial.integralInner_self] using T.outer_norm_real (2*b) (2*a)

/-- The normalized inner factors depend only on their common measure. -/
theorem inner_eq {μ μ' ν : Measure ℝ}
    (hν : Polynomial.HasFiniteMoments ν) (hνpos : Polynomial.HasPositiveDefiniteIntegral ν)
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ') ν) (n : ℕ) :
    S.inner n = T.inner n := by
  let c := (S.inner n).leadingCoeff / (T.inner n).leadingCoeff
  have hc : 0 < c := div_pos (S.inner_leadingCoeff_pos n) (T.inner_leadingCoeff_pos n)
  have he : S.inner n = c • T.inner n :=
    T.inner_orthogonal.eq_leadingCoeff_div_smul hν hνpos S.inner_orthogonal n
  have hn := S.inner_norm n
  rw [he] at hn
  simp only [Polynomial.eval_smul, smul_eq_mul, mul_pow, integral_const_mul, T.inner_norm] at hn
  have hc1 : c = 1 := by nlinarith
  simpa [hc1] using he

/-- Canonical even-even disk polynomials agree on corresponding deformed slices. -/
theorem eval_polynomial_diskDeformation
    {μ ν : Measure ℝ} [IsFiniteMeasure μ] [μ.IsNegInvariant]
    (hxμ : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    (hμ : Polynomial.HasFiniteMoments μ) (hμpos : Polynomial.HasPositiveDefiniteIntegral μ)
    (hν : Polynomial.HasFiniteMoments ν) (hνpos : Polynomial.HasPositiveDefiniteIntegral ν)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1)
    (hδpos : Polynomial.HasPositiveDefiniteIntegral (diskDeformedMeasure δ (μ.comap Subtype.val)))
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ (μ.comap Subtype.val))) ν)
    (a b : ℕ) (x t z : ℝ) (hx : x ∈ Icc (-1 : ℝ) 1) (ht : t ∈ Icc (-1 : ℝ) 1)
    (hxt : x^2 = 1-δ+δ*t^2) :
    eval ((fun v : ℝ × ℝ ↦ (![v.1, Real.sqrt (1-v.1^2)*v.2] : Fin 2 → ℝ)) (x,z)) (T.polynomial (2*a+2*b) (2*b)) =
      eval ((fun v : ℝ × ℝ ↦ (![v.1, Real.sqrt (1-v.1^2)*v.2] : Fin 2 → ℝ)) (t,z)) (S.polynomial (2*a+2*b) (2*b)) := by
  simp only [polynomial, Nat.add_sub_cancel]
  rw [← S.outer_even_diskDeformation_eq hxμ hμ hμpos hδ hδ1 hδpos T a b,
    ← S.inner_eq hν hνpos T (2*b)]
  apply eval_cartesianDiskPolynomial_evenDiskTransport hδ b _ _ _ (by simp)
    (S.inner_parity (2*b)) x t z hx ht hxt
  intro i hi
  let := isNegInvariant_diskOuterMeasure μ (2*b)
  exact (S.outer_orthogonal_real (2*b)).coeff_eq_zero_of_parity
    (hμ.diskOuterMeasure hxμ (2*b)) (hμpos.diskOuterMeasure hμ hxμ (2*b))
    (2*a) i (by simpa using hi)

/-- The modified Gram entries of the canonical even-even sector form an affine pencil. -/
theorem integral_quadratic_polynomial_diskDeformation
    {μ ν : Measure ℝ} [IsFiniteMeasure μ] [IsFiniteMeasure ν] [μ.IsNegInvariant]
    (hxμ : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    (hzν : ∀ᵐ z ∂ν, z ∈ Icc (-1 : ℝ) 1)
    (hμ : Polynomial.HasFiniteMoments μ) (hμpos : Polynomial.HasPositiveDefiniteIntegral μ)
    (hν : Polynomial.HasFiniteMoments ν) (hνpos : Polynomial.HasPositiveDefiniteIntegral ν)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1)
    (hδpos : Polynomial.HasPositiveDefiniteIntegral (diskDeformedMeasure δ (μ.comap Subtype.val)))
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ (μ.comap Subtype.val))) ν)
    (a b c d : ℕ) (A B C : ℝ) :
    (∫ v, (A+B*(v 0)^2+C*(v 1)^2)*
      (eval v (T.polynomial (2*a+2*b) (2*b))*eval v (T.polynomial (2*c+2*d) (2*d)))
      ∂cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ (μ.comap Subtype.val))) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ (μ.comap Subtype.val))) ν)) =
    ∫ v, ((A+B)+δ*(C*(v 1)^2-B*(1-(v 0)^2)))*
      (eval v (S.polynomial (2*a+2*b) (2*b))*eval v (S.polynomial (2*c+2*d) (2*d)))
      ∂cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) := by
  apply integral_quadratic_diskDeformedMeasure hxμ hzν hδ.le hδ1 A B C
    ((T.polynomial (2*a+2*b) (2*b)).continuous_eval.mul
      (T.polynomial (2*c+2*d) (2*d)).continuous_eval)
    ((S.polynomial (2*a+2*b) (2*b)).continuous_eval.mul
      (S.polynomial (2*c+2*d) (2*d)).continuous_eval)
  · intro t ht z _
    exact congrArg₂ (fun u v : ℝ => u*v)
      (S.eval_polynomial_diskDeformation hxμ hμ hμpos hν hνpos hδ hδ1 hδpos T
        a b _ t z (diskDeformation_mem_Icc_of_unitInterval hδ.le hδ1 ht) ht (diskDeformation_sq_of_unitInterval hδ.le hδ1 t))
      (S.eval_polynomial_diskDeformation hxμ hμ hμpos hν hνpos hδ hδ1 hδpos T
        c d _ t z (diskDeformation_mem_Icc_of_unitInterval hδ.le hδ1 ht) ht (diskDeformation_sq_of_unitInterval hδ.le hδ1 t))
  · intro t ht z _
    have hn : -diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ.le hδ1 t⟩ ∈ Icc (-1 : ℝ) 1 := by
      have hh := diskDeformation_mem_Icc_of_unitInterval hδ.le hδ1 ht
      constructor <;> linarith [hh.1,hh.2]
    have hs : (-diskDeformation δ ⟨t,diskDeformation_radical_nonneg hδ.le hδ1 t⟩)^2 = 1-δ+δ*t^2 := by
      simpa using diskDeformation_sq_of_unitInterval hδ.le hδ1 t
    exact congrArg₂ (fun u v : ℝ => u*v)
      (S.eval_polynomial_diskDeformation hxμ hμ hμpos hν hνpos hδ hδ1 hδpos T
        a b _ t z hn ht hs)
      (S.eval_polynomial_diskDeformation hxμ hμ hμpos hν hνpos hδ hδ1 hδpos T
        c d _ t z hn ht hs)

/-- Positive leading coefficients force the top-degree pairing to be positive. -/
theorem IsOrthogonal.integralInner_pos_of_degree_eq {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
    {P : Sequence ℝ} (hP : P.IsOrthogonal μ) {p : ℝ[X]} {n : ℕ}
    (hp : p.degree = n) (hlp : 0 < p.leadingCoeff) (hlP : 0 < (P n).leadingCoeff) :
    0 < integralInner μ (IntegrablePair.of μ p (P n)) := by
  rw [hP.integralInner_eq_leadingCoeff_div_mul hμ hp]
  apply mul_pos (div_pos hlp hlP)
  simpa [integralInner_self] using hpos (P n) (P.ne_zero n)

/-- Exact double raising, with the absent low-degree terms having zero coefficient. -/
theorem jacobi_symmetric_double_raise {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    jacobi β β n =
      jacobiSymmetricDoubleRaiseA β n • jacobi (β + 2) (β + 2) n +
      jacobiSymmetricDoubleRaiseB β n • jacobi (β + 2) (β + 2) (n - 2) +
      jacobiSymmetricDoubleRaiseC β n • jacobi (β + 2) (β + 2) (n - 4) := by
  have hb : -1 < β + 1 := by linarith
  rw [(jacobi_symmetric_raise_all hβ n).1,
    (jacobi_symmetric_raise_all hb n).1,
    (jacobi_symmetric_raise_all hb (n - 2)).1]
  simp only [jacobiSymmetricDoubleRaiseA, jacobiSymmetricDoubleRaiseB,
    jacobiSymmetricDoubleRaiseC, smul_add, add_smul, mul_smul, Nat.sub_sub,
    show β + 1 + 1 = β + 2 by ring]
  abel

/-- The highest remaining degree has a strictly positive coefficient. -/
theorem jacobiSymmetricDoubleRaiseA_pos {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    0 < jacobiSymmetricDoubleRaiseA β n :=
  mul_pos (jacobiSymmetricRaiseA_pos hβ n)
    (jacobiSymmetricRaiseA_pos (by linarith) n)

/-- The middle degree has a strictly negative coefficient whenever it exists. -/
theorem jacobiSymmetricDoubleRaiseB_neg {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    jacobiSymmetricDoubleRaiseB β (n + 2) < 0 := by
  unfold jacobiSymmetricDoubleRaiseB
  exact add_neg (mul_neg_of_pos_of_neg (jacobiSymmetricRaiseA_pos hβ _)
    (jacobiSymmetricRaiseB_neg (by linarith) n))
    (mul_neg_of_neg_of_pos (jacobiSymmetricRaiseB_neg hβ n)
      (jacobiSymmetricRaiseA_pos (by linarith) _))

/-- The lowest remaining degree has a strictly positive coefficient whenever it exists. -/
theorem jacobiSymmetricDoubleRaiseC_pos {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    0 < jacobiSymmetricDoubleRaiseC β (n + 4) := by
  unfold jacobiSymmetricDoubleRaiseC
  simp only [show n + 4 - 2 = n + 2 by omega]
  exact mul_pos_of_neg_of_neg (jacobiSymmetricRaiseB_neg hβ (n + 2))
    (jacobiSymmetricRaiseB_neg (by linarith) n)

/-- Exact three-term support of the shifted-weight pairing. -/
theorem integralInner_jacobi_double_raise {β : ℝ} (hβ : -1 < β) (n m : ℕ) :
    integralInner (jacobiMeasure (β + 2) (β + 2)) (IntegrablePair.of (jacobiMeasure (β + 2) (β + 2)) (jacobi β β n) (jacobi (β + 2) (β + 2) m) (h := (hasFiniteMoments_jacobiMeasure (by linarith) (by linarith)).integrable_mul _ _)) =
      (jacobiSymmetricDoubleRaiseA β n * (if n = m then 1 else 0) +
       jacobiSymmetricDoubleRaiseB β n * (if n - 2 = m then 1 else 0) +
       jacobiSymmetricDoubleRaiseC β n * (if n - 4 = m then 1 else 0)) *
        integralInner (jacobiMeasure (β + 2) (β + 2)) (IntegrablePair.of (jacobiMeasure (β + 2) (β + 2)) (jacobi (β + 2) (β + 2) m) (jacobi (β + 2) (β + 2) m) (h := (hasFiniteMoments_jacobiMeasure (by linarith) (by linarith)).integrable_mul _ _)) := by
  have hb : -1 < β + 2 := by linarith
  have hM := hasFiniteMoments_jacobiMeasure hb hb
  have he := congrArg (fun p : ℝ[X] => integralInner (jacobiMeasure (β+2) (β+2))
    (hM.pair p (jacobi (β+2) (β+2) m))) (jacobi_symmetric_double_raise hβ n)
  rw [he]
  simp only [integralInner_add_left (hasFiniteMoments_jacobiMeasure hb hb),
    integralInner_smul_left (h := hM.integrable_mul _ _)]
  rw [integralInner_jacobi_diagonal hb n m, integralInner_jacobi_diagonal hb (n - 2) m,
    integralInner_jacobi_diagonal hb (n - 4) m]
  split_ifs <;> ring

/-- Equal-degree shifted-weight pairings are strictly positive, including degrees zero and one. -/
theorem integralInner_jacobi_double_raise_top_pos {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    0 < integralInner (jacobiMeasure (β + 2) (β + 2)) (IntegrablePair.of (jacobiMeasure (β + 2) (β + 2)) (jacobi β β n) (jacobi (β + 2) (β + 2) n) (h := (hasFiniteMoments_jacobiMeasure (by linarith) (by linarith)).integrable_mul _ _)) := by
  rw [integralInner_jacobi_double_raise hβ]
  have hb : jacobiSymmetricDoubleRaiseB β n * (if n - 2 = n then 1 else 0) = 0 := by
    by_cases hn : n < 2
    · rw [jacobiSymmetricDoubleRaiseB_eq_zero hn, zero_mul]
    · rw [ite_eq_right (by omega), mul_zero]
  have hc : jacobiSymmetricDoubleRaiseC β n * (if n - 4 = n then 1 else 0) = 0 := by
    by_cases hn : n < 4
    · rw [jacobiSymmetricDoubleRaiseC_eq_zero hn, zero_mul]
    · rw [ite_eq_right (by omega), mul_zero]
  simp only [ite_true, mul_one, hb, hc, add_zero]
  exact mul_pos (jacobiSymmetricDoubleRaiseA_pos hβ n) (integralInner_jacobi_self_pos (by linarith) n)

/-- Dropping two degrees gives a strictly negative pairing. -/
theorem integralInner_jacobi_double_raise_middle_neg {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    integralInner (jacobiMeasure (β + 2) (β + 2)) (IntegrablePair.of (jacobiMeasure (β + 2) (β + 2)) (jacobi β β (n + 2)) (jacobi (β + 2) (β + 2) n) (h := (hasFiniteMoments_jacobiMeasure (by linarith) (by linarith)).integrable_mul _ _)) < 0 := by
  rw [integralInner_jacobi_double_raise hβ]
  have hc : jacobiSymmetricDoubleRaiseC β (n + 2) *
      (if n + 2 - 4 = n then 1 else 0) = 0 := by
    by_cases hn : n + 2 < 4
    · rw [jacobiSymmetricDoubleRaiseC_eq_zero hn, zero_mul]
    · rw [ite_eq_right (by omega), mul_zero]
  rw [ite_eq_right (by omega : n + 2 ≠ n)]
  simp only [mul_zero, zero_add, Nat.add_sub_cancel, ite_true, mul_one, hc, add_zero]
  exact mul_neg_of_neg_of_pos (jacobiSymmetricDoubleRaiseB_neg hβ n)
    (integralInner_jacobi_self_pos (by linarith) n)

/-- Dropping four degrees gives a strictly positive pairing. -/
theorem integralInner_jacobi_double_raise_bottom_pos {β : ℝ} (hβ : -1 < β) (n : ℕ) :
    0 < integralInner (jacobiMeasure (β + 2) (β + 2)) (IntegrablePair.of (jacobiMeasure (β + 2) (β + 2)) (jacobi β β (n + 4)) (jacobi (β + 2) (β + 2) n) (h := (hasFiniteMoments_jacobiMeasure (by linarith) (by linarith)).integrable_mul _ _)) := by
  rw [integralInner_jacobi_double_raise hβ, ite_eq_right (by omega : n + 4 ≠ n),
    ite_eq_right (by omega : n + 4 - 2 ≠ n)]
  simp only [mul_zero, zero_add, Nat.add_sub_cancel, ite_true, mul_one]
  exact mul_pos (jacobiSymmetricDoubleRaiseC_pos hβ n)
    (integralInner_jacobi_self_pos (by linarith) n)

/-- In the maximal inner-index step the radius power is exactly the target
outer weight, with no residual polynomial factor. -/
theorem integralInner_y_sq_dunklXuDisk_inner_add_two {β : ℝ} (hβ : -1 < β)
    (n k m : ℕ) :
    integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β n k) (dunklXuDisk β m (k + 2))) =
    Polynomial.integralInner (jacobiMeasure (β + (k + 2 : ℕ) + 1 / 2) (β + (k + 2 : ℕ) + 1 / 2)) (Polynomial.IntegrablePair.of (jacobiMeasure (β + (k + 2 : ℕ) + 1 / 2) (β + (k + 2 : ℕ) + 1 / 2)) (jacobi (β + k + 1 / 2) (β + k + 1 / 2) (n - k)) (jacobi (β + (k + 2 : ℕ) + 1 / 2) (β + (k + 2 : ℕ) + 1 / 2) (m - (k + 2))) (h := (hasFiniteMoments_jacobiMeasure (by linarith [show (0 : ℝ) ≤ k from Nat.cast_nonneg k]) (by linarith [show (0 : ℝ) ≤ k from Nat.cast_nonneg k])).integrable_mul _ _)) *
    Polynomial.integralInner (jacobiMeasure β β) (Polynomial.IntegrablePair.of (jacobiMeasure β β) (Polynomial.X ^ 2 * jacobi β β k) (jacobi β β (k + 2))) := by
  rw [integralInner_y_sq_dunklXuDisk_separated hβ]
  congr 1
  let a := β + ((k + 2 : ℕ) : ℝ) + 1/2
  have ha : -1 < a := by dsimp [a]; have := Nat.cast_nonneg (α := ℝ) (k+2); linarith
  rw [Polynomial.integralInner, integral_jacobiMeasure a a
    ⟨_, (hasFiniteMoments_jacobiMeasure ha ha).integrable_mul _ _⟩]
  apply integral_congr_ae
  filter_upwards [] with x
  have hs : 0 ≤ 1 - (x : ℝ)^2 := by nlinarith [x.property.1, x.property.2]
  rw [show a = (β + 1/2) + (k+2 : ℕ) by dsimp [a]; ring,
    jacobiWeight_add_nat,
    show k + (k + 2) + 2 = 2 * (k + 2) by omega, pow_mul, Real.sq_sqrt hs]
  ring

/-- The inner maximal-degree pairing is strictly positive. -/
theorem integralInner_X_sq_jacobi_add_two_pos {β : ℝ} (hβ : -1 < β) (k : ℕ) :
    0 < Polynomial.integralInner (jacobiMeasure β β) (Polynomial.IntegrablePair.of (jacobiMeasure β β) (Polynomial.X ^ 2 * jacobi β β k) (jacobi β β (k + 2))) := by
  have h := (isOrthogonal_jacobiSequence hβ hβ).integralInner_pos_of_degree_eq
    (hasFiniteMoments_jacobiMeasure hβ hβ)
    (hasPositiveDefiniteIntegral_jacobiMeasure hβ hβ)
    (p := Polynomial.X ^ 2 * jacobi β β k) (n := k + 2)
  apply h
  · rw [degree_mul, degree_pow, degree_X,
      degree_jacobi_of_neg_one_lt hβ hβ]
    push_cast
    ring
  · rw [leadingCoeff_mul, leadingCoeff_pow, leadingCoeff_X, one_pow, one_mul]
    rw [leadingCoeff, show (jacobi β β k).natDegree = k from (jacobiSequence hβ hβ).natDegree_eq k]
    exact coeff_jacobi_self_pos hβ hβ k
  · change 0 < (jacobi β β (k + 2)).leadingCoeff
    rw [leadingCoeff, show (jacobi β β (k + 2)).natDegree = k + 2 from
      (jacobiSequence hβ hβ).natDegree_eq (k + 2)]
    exact coeff_jacobi_self_pos hβ hβ (k + 2)

/-- Simultaneously increasing total degree and inner index by two has positive pairing. -/
theorem integralInner_y_sq_dunklXuDisk_up_pos {β : ℝ} (hβ : -1 < β) (n k : ℕ) :
    0 < integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β (n + k) k) (dunklXuDisk β (n + (k + 2)) (k + 2))) := by
  rw [integralInner_y_sq_dunklXuDisk_inner_add_two hβ]
  simp only [Polynomial.integralInner, Nat.add_sub_cancel, shifted_outerParameter]
  simpa only [Polynomial.integralInner] using mul_pos (integralInner_jacobi_double_raise_top_pos (outerParameter_pos hβ k) n)
    (integralInner_X_sq_jacobi_add_two_pos hβ k)

/-- Increasing the inner index while keeping the total degree gives a negative pairing. -/
theorem integralInner_y_sq_dunklXuDisk_same_degree_neg {β : ℝ} (hβ : -1 < β) (n k : ℕ) :
    integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β (n + 2 + k) k) (dunklXuDisk β (n + (k + 2)) (k + 2))) < 0 := by
  rw [integralInner_y_sq_dunklXuDisk_inner_add_two hβ]
  simp only [Polynomial.integralInner, Nat.add_sub_cancel, shifted_outerParameter]
  simpa only [Polynomial.integralInner] using mul_neg_of_neg_of_pos
    (integralInner_jacobi_double_raise_middle_neg (outerParameter_pos hβ k) n)
    (integralInner_X_sq_jacobi_add_two_pos hβ k)

/-- Increasing the inner index and decreasing total degree by two has positive pairing. -/
theorem integralInner_y_sq_dunklXuDisk_down_pos {β : ℝ} (hβ : -1 < β) (n k : ℕ) :
    0 < integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β (n + 4 + k) k) (dunklXuDisk β (n + (k + 2)) (k + 2))) := by
  rw [integralInner_y_sq_dunklXuDisk_inner_add_two hβ]
  simp only [Polynomial.integralInner, Nat.add_sub_cancel, shifted_outerParameter]
  simpa only [Polynomial.integralInner] using mul_pos (integralInner_jacobi_double_raise_bottom_pos (outerParameter_pos hβ k) n)
    (integralInner_X_sq_jacobi_add_two_pos hβ k)

/-- No other total degrees occur at the maximal inner-index step. -/
theorem integralInner_y_sq_dunklXuDisk_inner_add_two_eq_zero {β : ℝ} (hβ : -1 < β)
    (n m k : ℕ) (h₀ : n ≠ m) (h₂ : n - 2 ≠ m) (h₄ : n - 4 ≠ m) :
    integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β (n + k) k) (dunklXuDisk β (m + (k + 2)) (k + 2))) = 0 := by
  rw [integralInner_y_sq_dunklXuDisk_inner_add_two hβ]
  simp only [Polynomial.integralInner, Nat.add_sub_cancel, shifted_outerParameter]
  have he := integralInner_jacobi_double_raise (outerParameter_pos hβ k) n m
  simp only [Polynomial.integralInner] at he
  rw [he]
  simp [h₀, h₂, h₄]

/-- The normalized quadratic perturbation matrix; its entries are symmetric pairings. -/
noncomputable def dunklXuQuadraticMatrix {β : ℝ} (hβ : -1 < β) (b c : ℝ) :
    Matrix (Σ n : ℕ, Fin (n + 1)) (Σ n : ℕ, Fin (n + 1)) ℝ := fun i j ↦
  integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) ((C c * X 1 ^ 2 - C b * (1 - X 0 ^ 2)) * normalizedDunklXuBasis hβ i) (normalizedDunklXuBasis hβ j))

/-- Positive norm scaling relates a normalized multiplication entry to its raw pairing. -/
theorem integralInner_mul_normalizedDunklXuBasis {β : ℝ} (hβ : -1 < β)
    (u : MvPolynomial (Fin 2) ℝ) (i j : Σ n : ℕ, Fin (n + 1)) :
    integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (u * normalizedDunklXuBasis hβ i) (normalizedDunklXuBasis hβ j)) =
      (dunklXuNorm β ⟨(i.1, i.2), by exact hasFiniteMoments_dunklXuDiskMeasure (by linarith) _⟩)⁻¹ * (dunklXuNorm β ⟨(j.1, j.2), by exact hasFiniteMoments_dunklXuDiskMeasure (by linarith) _⟩)⁻¹ *
        integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (u * dunklXuDisk β i.1 i.2) (dunklXuDisk β j.1 j.2)) := by
  simp only [normalizedDunklXuBasis_apply, mul_smul_comm]
  rw [integralInner_smul_left (h := (hasFiniteMoments_dunklXuDiskMeasure hβ).integrable_mul _ _), integralInner_comm,
    integralInner_smul_left (h := (hasFiniteMoments_dunklXuDiskMeasure hβ).integrable_mul _ _), integralInner_comm]
  ring

/-- The matrix is symmetric because multiplication is symmetric for the integral pairing. -/
theorem dunklXuQuadraticMatrix_symm {β : ℝ} (hβ : -1 < β) (b c : ℝ)
    (i j : Σ n : ℕ, Fin (n + 1)) :
    dunklXuQuadraticMatrix hβ b c i j = dunklXuQuadraticMatrix hβ b c j i := by
  unfold dunklXuQuadraticMatrix
  rw [integralInner_mul_left_right, integralInner_comm]

/-- Across different inner indices the outer-coordinate part vanishes exactly. -/
theorem dunklXuQuadraticMatrix_of_inner_ne {β : ℝ} (hβ : -1 < β) (b c : ℝ)
    (i j : Σ n : ℕ, Fin (n + 1)) (hij : (i.2 : ℕ) ≠ (j.2 : ℕ)) :
    dunklXuQuadraticMatrix hβ b c i j =
      c * ((dunklXuNorm β ⟨(i.1, i.2), by exact hasFiniteMoments_dunklXuDiskMeasure (by linarith) _⟩)⁻¹ * (dunklXuNorm β ⟨(j.1, j.2), by exact hasFiniteMoments_dunklXuDiskMeasure (by linarith) _⟩)⁻¹) *
        integralInner (dunklXuDiskMeasure β) (IntegrablePair.of (dunklXuDiskMeasure β) (X 1 ^ 2 * dunklXuDisk β i.1 i.2) (dunklXuDisk β j.1 j.2)) := by
  rw [dunklXuQuadraticMatrix, integralInner_mul_normalizedDunklXuBasis]
  simp only [sub_mul]
  rw [integralInner_sub_left (hasFiniteMoments_dunklXuDiskMeasure hβ)]
  have he : C b * (1 - X 0 ^ 2) =
      (Polynomial.C b * (1 - Polynomial.X ^ 2)).eval₂ C (X 0 : MvPolynomial (Fin 2) ℝ) := by
    simp
  simp only [he]
  rw [integralInner_outer_mul_dunklXuDisk_eq_zero hβ _ _ _ _ _ hij, sub_zero]
  have hs : C c * X 1 ^ 2 * dunklXuDisk β i.1 i.2 =
    c • (X 1 ^ 2 * dunklXuDisk β i.1 i.2) := by rw [smul_eq_C_mul, mul_assoc]
  simp only [hs]
  rw [integralInner_smul_left (h := (hasFiniteMoments_dunklXuDiskMeasure hβ).integrable_mul _ _)]
  ring

/-- The upward inner bandwidth is two, independently of the outer-coordinate coefficient. -/
theorem dunklXuQuadraticMatrix_eq_zero_of_inner_gap {β : ℝ} (hβ : -1 < β) (b c : ℝ)
    (i j : Σ n : ℕ, Fin (n + 1)) (hij : (i.2 : ℕ) + 2 < (j.2 : ℕ)) :
    dunklXuQuadraticMatrix hβ b c i j = 0 := by
  rw [dunklXuQuadraticMatrix_of_inner_ne hβ b c i j (by omega),
    integralInner_y_sq_dunklXuDisk_eq_zero_of_inner_gap hβ _ _ _ _ hij, mul_zero]

/-- Opposite inner parities do not couple in the quadratic modification. -/
theorem dunklXuQuadraticMatrix_eq_zero_of_inner_parity {β : ℝ} (hβ : -1 < β) (b c : ℝ)
    (i j : Σ n : ℕ, Fin (n + 1)) (hij : (i.2 : ℕ) % 2 ≠ (j.2 : ℕ) % 2) :
    dunklXuQuadraticMatrix hβ b c i j = 0 := by
  rw [dunklXuQuadraticMatrix_of_inner_ne hβ b c i j (by omega),
    integralInner_y_sq_dunklXuDisk_eq_zero_of_inner_parity hβ _ _ _ _ hij, mul_zero]

/-- At the first possible level gap, retaining only maximal boundary and
interior steps leaves the sandwich unchanged. -/
theorem sum_mul_pow_mul_eq_maximalStep
    (K : Matrix ι ι ℝ) (F G : ι → ℝ) (level : ι → ℤ) (lo hi : ℤ)
    (hK : ∀ i j, level i + 1 < level j → K i j = 0)
    (hF : ∀ i, lo + 1 < level i → F i = 0)
    (hG : ∀ i, level i + 1 < hi → G i = 0)
    (q : ℕ) (hgap : hi = lo + q + 2) :
    (∑ i, ∑ j, F i * (K ^ q) i j * G j) =
      ∑ i, ∑ j, (if level i = lo + 1 then F i else 0) *
        ((maximalStep K level) ^ q) i j * (if level j + 1 = hi then G j else 0) := by
  apply Finset.sum_congr rfl
  intro i _
  apply Finset.sum_congr rfl
  intro j _
  by_cases hf : F i = 0
  · simp [hf]
  by_cases hg : G j = 0
  · simp [hg]
  have hfi : level i ≤ lo + 1 := by
    by_contra h
    exact hf (hF i (by omega))
  have hgj : hi ≤ level j + 1 := by
    by_contra h
    exact hg (hG j (by omega))
  by_cases hi' : level i = lo + 1
  · by_cases hj' : level j + 1 = hi
    · rw [ite_eq_left hi', ite_eq_left hj',
        pow_apply_eq_maximalStep K level hK q i j (by omega)]
    · rw [ite_eq_right hj', pow_apply_eq_zero_of_level_gap K level hK q i j (by omega)]
      simp
  · rw [ite_eq_right hi', pow_apply_eq_zero_of_level_gap K level hK q i j (by omega)]
    simp

/-- Diagonal signs telescope through both boundary edges and the interior power. -/
theorem sum_mul_pow_mul_sign_conjugate
    (K : Matrix ι ι ℝ) (F G s : ι → ℝ) (slo shi : ℝ)
    (hs : ∀ i, s i * s i = 1) (q : ℕ) :
    (∑ i, ∑ j, (slo * F i * s i) *
      ((Matrix.of (fun i j ↦ s i * K i j * s j)) ^ q) i j * (s j * G j * shi)) =
      slo * (∑ i, ∑ j, F i * (K ^ q) i j * G j) * shi := by
  simp only [Finset.mul_sum, Finset.sum_mul]
  apply Finset.sum_congr rfl
  intro i _
  apply Finset.sum_congr rfl
  intro j _
  rw [pow_apply_sign_conjugate K s hs]
  calc
    _ = (slo * (F i * (K ^ q) i j * G j) * shi) * (s i * s i) * (s j * s j) := by ring
    _ = _ := by rw [hs, hs, mul_one, mul_one]

/-- Nonnegative factors and one strictly positive path give a positive sandwich. -/
theorem sum_mul_pow_mul_pos_of_positive_path
    (K : Matrix ι ι ℝ) (F G : ι → ℝ)
    (hK : ∀ i j, 0 ≤ K i j) (hF : ∀ i, 0 ≤ F i) (hG : ∀ i, 0 ≤ G i)
    (q : ℕ) (p : ℕ → ι) (hfirst : 0 < F (p 0)) (hlast : 0 < G (p q))
    (hp : ∀ k < q, 0 < K (p k) (p (k + 1))) :
    0 < ∑ i, ∑ j, F i * (K ^ q) i j * G j := by
  have ht := mul_pos (mul_pos hfirst (pow_apply_pos_of_positive_path K hK q p hp)) hlast
  apply ht.trans_le
  calc
    _ ≤ ∑ j, F (p 0) * (K ^ q) (p 0) j * G j :=
      Finset.single_le_sum (fun j _ ↦ mul_nonneg
        (mul_nonneg (hF _) (pow_apply_nonneg hK q _ j)) (hG j)) (Finset.mem_univ _)
    _ ≤ ∑ i, ∑ j, F i * (K ^ q) i j * G j :=
      Finset.single_le_sum (fun i _ ↦ Finset.sum_nonneg (fun j _ ↦ mul_nonneg
        (mul_nonneg (hF i) (pow_apply_nonneg hK q i j)) (hG j))) (Finset.mem_univ _)

/-- The boundary-to-boundary pairing of an interior matrix. -/
noncomputable def dunklXuInteriorFunctional {β : ℝ} (hβ : -1 < β) (u v : ℝ) (r : ℕ) :
    Matrix (DunklXuEvenInterior r) (DunklXuEvenInterior r) ℝ →ₗ[ℝ] ℝ where
  toFun M := ∑ i, ∑ j, dunklXuInteriorLeft hβ u v r i * M i j * dunklXuInteriorRight hβ u v r j
  map_add' M N := by simp [mul_add, add_mul, Finset.sum_add_distrib]
  map_smul' a M := by simp [Finset.mul_sum, mul_assoc, mul_left_comm]

/-- Index embedding for the even-even Cartesian disk sector. -/
def dunklXuEvenIndex (a b : ℕ) : Σ n : ℕ, Fin (n + 1) :=
  ⟨2 * a + 2 * b, ⟨2 * b, by omega⟩⟩

/-- Even-even indices strictly below total degree 2r. -/
abbrev DunklXuEvenInterior (r : ℕ) :=
  {i : Fin r × Fin r // (i.1 : ℕ) + (i.2 : ℕ) < r}

/-- The explicit path stays in degree 2r-2 between its two boundary edges. -/
def dunklXuInteriorPath (r : ℕ) (hr : 2 ≤ r) (k : ℕ) : DunklXuEvenInterior r :=
  ⟨(⟨r - 2 - min k (r - 2), by omega⟩, ⟨min k (r - 2) + 1, by omega⟩), by dsimp; omega⟩

/-- Before the final interior vertex, each edge lowers the outer degree by two
and raises the inner degree by two. -/
theorem dunklXuInteriorPath_step (r : ℕ) (hr : 2 ≤ r) (k : ℕ) (hk : k < r - 2) :
    ((dunklXuInteriorPath r hr k).1.1 : ℕ) =
        ((dunklXuInteriorPath r hr (k + 1)).1.1 : ℕ) + 1 ∧
      ((dunklXuInteriorPath r hr (k + 1)).1.2 : ℕ) =
        ((dunklXuInteriorPath r hr k).1.2 : ℕ) + 1 := by
  dsimp [dunklXuInteriorPath]
  omega

/-- Diagonal sign changes remove every sign of a maximal inner step. -/
theorem dunklXuEvenY2Entry_signed_nonneg {β : ℝ} (hβ : -1 < β) (a b c : ℕ) :
    0 ≤ (-1 : ℝ) ^ a * dunklXuEvenY2Entry hβ a b c (b + 1) * (-1 : ℝ) ^ c := by
  have hs : (-1 : ℝ) ^ c * (-1 : ℝ) ^ c = 1 := by rw [← mul_pow]; norm_num
  by_cases h₀ : a = c
  · subst a
    have hp := (dunklXuEvenY2Entry_up_pos hβ c b).le
    nlinarith [hs]
  by_cases h₁ : a = c + 1
  · subst a
    rw [pow_add, pow_one]
    have hp := (dunklXuEvenY2Entry_middle_neg hβ c b).le
    nlinarith [hs]
  by_cases h₂ : a = c + 2
  · subst a
    rw [pow_add]
    norm_num only [neg_one_sq, mul_one]
    have hp := (dunklXuEvenY2Entry_down_pos hβ c b).le
    nlinarith [hs]
  rw [dunklXuEvenY2Entry_eq_zero hβ a b c h₀ h₁ h₂]
  simp

/-- Every elimination power before the first possible inner gap vanishes. -/
theorem dunklXuInterior_sandwich_eq_zero {β : ℝ} (hβ : -1 < β) (u v : ℝ) (r q : ℕ)
    (hq : q + 2 < r) :
    ∑ i, ∑ j, dunklXuInteriorLeft hβ u v r i *
      (dunklXuInteriorMatrix hβ u v r ^ q) i j * dunklXuInteriorRight hβ u v r j = 0 := by
  apply sum_mul_pow_mul_eq_zero_of_level_gap _ _ _ dunklXuInteriorLevel 0 r
    (dunklXuInteriorMatrix_level_band hβ u v r)
    (dunklXuInteriorLeft_level_band hβ u v r)
    (dunklXuInteriorRight_level_band hβ u v r)
  omega

/-- The first allowed elimination coefficient is nonzero in every degree. -/
theorem dunklXuInterior_first_coefficient_ne_zero {β : ℝ} (hβ : -1 < β) (u : ℝ)
    (r : ℕ) (hr : 2 ≤ r) :
    dunklXuInteriorFunctional hβ u 1 r (dunklXuInteriorMatrix hβ u 1 r ^ (r - 2)) ≠ 0 := by
  let K := maximalStep (dunklXuInteriorMatrix hβ u 1 r) dunklXuInteriorLevel
  let F : DunklXuEvenInterior r → ℝ := fun i ↦
    if dunklXuInteriorLevel i = 1 then dunklXuInteriorLeft hβ u 1 r i else 0
  let G : DunklXuEvenInterior r → ℝ := fun i ↦
    if dunklXuInteriorLevel i + 1 = r then dunklXuInteriorRight hβ u 1 r i else 0
  have hp := sum_mul_pow_mul_pos_of_positive_path
    (Matrix.of (fun i j ↦ dunklXuInteriorSign i * K i j * dunklXuInteriorSign j))
    (fun i ↦ (-1 : ℝ) ^ r * F i * dunklXuInteriorSign i)
    (fun j ↦ dunklXuInteriorSign j * G j * 1)
    (dunklXuInterior_signed_maximal_nonneg hβ u r)
    (dunklXuInterior_signed_left_nonneg hβ u r)
    (fun i ↦ by simpa only [mul_one] using dunklXuInterior_signed_right_nonneg hβ u r i)
    (r - 2) (dunklXuInteriorPath r hr)
    (dunklXuInterior_signed_first_pos hβ u r hr)
    (by simpa only [mul_one] using dunklXuInterior_signed_last_pos hβ u r hr)
    (dunklXuInterior_signed_path_pos hβ u r hr)
  rw [sum_mul_pow_mul_sign_conjugate K F G dunklXuInteriorSign ((-1 : ℝ) ^ r) 1
    dunklXuInteriorSign_sq, mul_one] at hp
  have he := sum_mul_pow_mul_eq_maximalStep (dunklXuInteriorMatrix hβ u 1 r)
    (dunklXuInteriorLeft hβ u 1 r) (dunklXuInteriorRight hβ u 1 r)
    dunklXuInteriorLevel 0 r (dunklXuInteriorMatrix_level_band hβ u 1 r)
    (dunklXuInteriorLeft_level_band hβ u 1 r) (dunklXuInteriorRight_level_band hβ u 1 r)
    (r - 2) (by omega)
  change (∑ i, ∑ j, dunklXuInteriorLeft hβ u 1 r i *
    (dunklXuInteriorMatrix hβ u 1 r ^ (r - 2)) i j * dunklXuInteriorRight hβ u 1 r j) =
      ∑ i, ∑ j, F i * (K ^ (r - 2)) i j * G j at he
  rw [← he] at hp
  intro hz
  change (∑ i, ∑ j, dunklXuInteriorLeft hβ u 1 r i *
    (dunklXuInteriorMatrix hβ u 1 r ^ (r - 2)) i j * dunklXuInteriorRight hβ u 1 r j) = 0 at hz
  rw [hz, mul_zero] at hp
  exact lt_irrefl 0 hp

/-- The first possible power sandwich scales by the coefficient once per edge. -/
theorem dunklXuInteriorFunctional_pow_scale {β : ℝ} (hβ : -1 < β)
    (u v : ℝ) (hv : v ≠ 0) (r q : ℕ) :
    dunklXuInteriorFunctional hβ u v r (dunklXuInteriorMatrix hβ u v r ^ q) =
      v ^ (q + 2) * dunklXuInteriorFunctional hβ (u / v) 1 r
        (dunklXuInteriorMatrix hβ (u / v) 1 r ^ q) := by
  have hK : dunklXuInteriorMatrix hβ u v r = v • dunklXuInteriorMatrix hβ (u / v) 1 r := by
    ext i j
    exact dunklXuQuadraticMatrix_scale hβ u v hv _ _
  have hF (i : DunklXuEvenInterior r) : dunklXuInteriorLeft hβ u v r i =
      v * dunklXuInteriorLeft hβ (u / v) 1 r i := dunklXuQuadraticMatrix_scale hβ u v hv _ _
  have hG (i : DunklXuEvenInterior r) : dunklXuInteriorRight hβ u v r i =
      v * dunklXuInteriorRight hβ (u / v) 1 r i := dunklXuQuadraticMatrix_scale hβ u v hv _ _
  rw [hK, smul_pow]
  change (∑ i, ∑ j, dunklXuInteriorLeft hβ u v r i *
    (v ^ q * (dunklXuInteriorMatrix hβ (u / v) 1 r ^ q) i j) * dunklXuInteriorRight hβ u v r j) = _
  simp_rw [hF, hG]
  change _ = v ^ (q + 2) * (∑ i, ∑ j, dunklXuInteriorLeft hβ (u / v) 1 r i *
    (dunklXuInteriorMatrix hβ (u / v) 1 r ^ q) i j * dunklXuInteriorRight hβ (u / v) 1 r j)
  simp only [Finset.mul_sum]
  apply Finset.sum_congr rfl
  intro i _
  apply Finset.sum_congr rfl
  intro j _
  rw [pow_add]
  ring

/-- Every nonzero y² coefficient gives a surviving first path coefficient,
regardless of its sign or of the x-only quadratic coefficient. -/
theorem dunklXuInterior_first_coefficient_ne_zero_of_vertical_ne_zero {β : ℝ}
    (hβ : -1 < β) (u v : ℝ) (hv : v ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    dunklXuInteriorFunctional hβ u v r (dunklXuInteriorMatrix hβ u v r ^ (r - 2)) ≠ 0 := by
  rw [dunklXuInteriorFunctional_pow_scale hβ u v hv]
  exact mul_ne_zero (pow_ne_zero _ hv) (dunklXuInterior_first_coefficient_ne_zero hβ (u / v) r hr)

/-- The finite inverse sandwich remains nonzero locally for any nonzero
vertical coefficient; no continuity of a Cholesky factor is assumed. -/
theorem eventually_dunklXuInterior_inverse_ne_zero_of_vertical_ne_zero {β : ℝ}
    (hβ : -1 < β) (u v : ℝ) (hv : v ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ t in 𝓝[≠] (0 : ℝ),
      dunklXuInteriorFunctional hβ u v r ((1 + t • dunklXuInteriorMatrix hβ u v r)⁻¹) ≠ 0 := by
  apply eventually_apply_inv_one_add_smul_ne_zero _ _ (r - 2)
  · intro q hq
    exact dunklXuInterior_sandwich_eq_zero hβ u v r q (by omega)
  · exact dunklXuInterior_first_coefficient_ne_zero_of_vertical_ne_zero hβ u v hv r hr

/-- Perturbation block from the top even degree to the lower degree section. -/
noncomputable def dunklXuQuadraticLeftBlock {β : ℝ} (hβ : -1 < β) (u v : ℝ) (r : ℕ) :
    Matrix (Fin (r + 1)) (DunklXuEvenInterior r) ℝ := fun i j ↦
  dunklXuQuadraticMatrix hβ u v (dunklXuEvenBoundary r i) (dunklXuEvenIndex j.1.1 j.1.2)

/-- Perturbation block from the lower degree section to the top even degree. -/
noncomputable def dunklXuQuadraticRightBlock {β : ℝ} (hβ : -1 < β) (u v : ℝ) (r : ℕ) :
    Matrix (DunklXuEvenInterior r) (Fin (r + 1)) ℝ := fun i j ↦
  dunklXuQuadraticMatrix hβ u v (dunklXuEvenIndex i.1.1 i.1.2) (dunklXuEvenBoundary r j)

/-- Perturbation within the top even degree. -/
noncomputable def dunklXuQuadraticBoundaryBlock {β : ℝ} (hβ : -1 < β) (u v : ℝ) (r : ℕ) :
    Matrix (Fin (r + 1)) (Fin (r + 1)) ℝ := fun i j ↦
  dunklXuQuadraticMatrix hβ u v (dunklXuEvenBoundary r i) (dunklXuEvenBoundary r j)

/-- The uneliminated first-to-last top-degree entry is zero for r>=2. -/
theorem dunklXuQuadraticBoundaryBlock_first_last {β : ℝ} (hβ : -1 < β)
    (u v : ℝ) (r : ℕ) (hr : 2 ≤ r) :
    dunklXuQuadraticBoundaryBlock hβ u v r 0 (Fin.last r) = 0 := by
  apply dunklXuQuadraticMatrix_eq_zero_of_inner_gap
  simp only [dunklXuEvenBoundary, dunklXuEvenIndex, Fin.val_zero, Fin.val_last]
  omega

/-- Every compatible triangular Gram factor has a nonzero first-to-last entry
in its top even-degree block near zero. The actual disk matrices instantiate
all path hypotheses; only the displayed factor identities remain assumptions. -/
theorem eventually_dunklXuQuadratic_cholesky_first_last_ne_zero {β : ℝ}
    (hβ : -1 < β) (u v : ℝ) (hv : v ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ t in 𝓝[≠] (0 : ℝ),
      ∀ (A : Matrix (DunklXuEvenInterior r) (DunklXuEvenInterior r) ℝ)
        (B : Matrix (DunklXuEvenInterior r) (Fin (r + 1)) ℝ)
        (D : Matrix (Fin (r + 1)) (Fin (r + 1)) ℝ),
        IsUnit A.det → D.IsUpperTriangular →
        A.transpose * A = 1 + t • dunklXuInteriorMatrix hβ u v r →
        A.transpose * B = t • dunklXuQuadraticRightBlock hβ u v r →
        B.transpose * A = t • dunklXuQuadraticLeftBlock hβ u v r →
        B.transpose * B + D.transpose * D = 1 + t • dunklXuQuadraticBoundaryBlock hβ u v r →
        D 0 (Fin.last r) ≠ 0 := by
  have hn := eventually_cholesky_first_row_ne_zero (dunklXuInteriorMatrix hβ u v r)
    (dunklXuQuadraticRightBlock hβ u v r) (dunklXuQuadraticLeftBlock hβ u v r)
    (Fin.last r) (r - 2)
    (fun q hq ↦ by
      rw [← dunklXuInteriorFunctional_eq_sandwich]
      exact dunklXuInterior_sandwich_eq_zero hβ u v r q (by omega))
    (by rw [← dunklXuInteriorFunctional_eq_sandwich]
        exact dunklXuInterior_first_coefficient_ne_zero_of_vertical_ne_zero hβ u v hv r hr)
  filter_upwards [hn] with t ht A B D hA hD hAA hAB hBA hDD
  apply ht A B D hA hD hAA hAB hBA
  rw [hDD, Matrix.add_apply, Matrix.smul_apply,
    dunklXuQuadraticBoundaryBlock_first_last hβ u v r hr]
  have he : (0 : Fin (r + 1)) ≠ Fin.last r := by
    intro h
    have := congrArg Fin.val h
    simp only [Fin.val_zero, Fin.val_last] at this
    omega
  simp [he]

/-- Restricting to a complete sector commutes with taking the Gram matrix,
even when the sector's indices are interleaved with other sectors. -/
theorem IsUpperTriangular.transpose_mul_self_sector {n : ℕ} {ι γ : Type*} [Fintype ι]
    {R : Matrix (Fin n) (Fin n) ℝ} (hR : R.IsUpperTriangular)
    (hd : ∀ i, R i i ≠ 0) (sector : Fin n → γ)
    (hG : ∀ i j, sector i ≠ sector j → (R.transpose * R) i j = 0)
    (e : ι → Fin n) (he : Function.Injective e) (c : γ)
    (hc : ∀ i, sector (e i) = c) (hcover : ∀ k, sector k = c → k ∈ Set.range e) :
    (R.transpose * R).submatrix e e = (R.submatrix e e).transpose * R.submatrix e e := by
  ext i j
  simp only [submatrix_apply, mul_apply, transpose_apply]
  symm
  apply Fintype.sum_of_injective e he
  · intro k hk
    have hki : sector k ≠ sector (e i) := by
      rw [hc]
      exact fun h ↦ hk (hcover k h)
    rw [hR.apply_eq_zero_of_gram_sectors hd sector hG k (e i) hki, zero_mul]
  · intro k
    rfl

/-- A complete orthogonal sector split into earlier and later indices has the
usual upper block Gram formula, independently of interleaving outside the sector. -/
theorem IsUpperTriangular.gram_sector_blocks {n : ℕ} {ι κ γ : Type*}
    [Fintype ι] [Fintype κ] {R : Matrix (Fin n) (Fin n) ℝ}
    (hR : R.IsUpperTriangular) (hd : ∀ i, R i i ≠ 0) (sector : Fin n → γ)
    (hG : ∀ i j, sector i ≠ sector j → (R.transpose * R) i j = 0)
    (e : ι ⊕ κ → Fin n) (he : Function.Injective e) (c : γ)
    (hc : ∀ i, sector (e i) = c) (hcover : ∀ k, sector k = c → k ∈ Set.range e)
    (hbefore : ∀ i j, e (Sum.inl i) < e (Sum.inr j)) :
    let A := R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inl)
    let B := R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inr)
    let D := R.submatrix (e ∘ Sum.inr) (e ∘ Sum.inr)
    (R.transpose * R).submatrix e e =
      fromBlocks (A.transpose * A) (A.transpose * B) (B.transpose * A)
        (B.transpose * B + D.transpose * D) := by
  classical
  dsimp only
  rw [hR.transpose_mul_self_sector hd sector hG e he c hc hcover]
  have hs : R.submatrix e e = fromBlocks
      (R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inl))
      (R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inr)) 0
      (R.submatrix (e ∘ Sum.inr) (e ∘ Sum.inr)) := by
    ext i j
    cases i with
    | inl i => cases j <;> rfl
    | inr i =>
      cases j with
      | inl j => exact hR (hbefore j i)
      | inr j => rfl
  rw [hs, fromBlocks_transpose, fromBlocks_multiply]
  simp

/-- The modified Gram entry of actual deformed Cartesian polynomials equals
the block pencil when the constant boundary value A+B is normalized to one. -/
theorem integralInner_deformed_quadratic_eq_evenGramPencil {β : ℝ} (hβ : -1 < β)
    {δ : ℝ} (hδ : 0 < δ) (hδ1 : δ ≤ 1)
    (hδpos : Polynomial.HasPositiveDefiniteIntegral
      (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val)))
    (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β))
    (A B C : ℝ) (hAB : A + B = 1)
    (hu : ∀ᵐ x ∂cartesianDiskMeasure
      (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β)),
      0 ≤ eval x (coordinateQuadratic A B C))
    (r : ℕ) (i j : DunklXuEvenInterior r ⊕ Fin (r + 1)) :
    let ci := dunklXuEvenBlockCoordinates r i
    let cj := dunklXuEvenBlockCoordinates r j
    integralInner (MvPolynomial.polynomialDensityMeasure ((coordinateQuadratic A B C)) (((cartesianDiskMeasure
      (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val)))
      (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β)))).comap Subtype.val)) (IntegrablePair.of (MvPolynomial.polynomialDensityMeasure ((coordinateQuadratic A B C)) (((cartesianDiskMeasure
      (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val)))
      (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β)))).comap Subtype.val)) (T.polynomial (2 * ci.1 + 2 * ci.2) (2 * ci.2)) (T.polynomial (2 * cj.1 + 2 * cj.2) (2 * cj.2)) (h := ((T.hasFiniteMoments (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 (ae_mem_Icc_jacobiMeasure _ _))).withDensity_polynomial _ hu).integrable_mul _ _)) =
        dunklXuEvenGramPencil hβ B C δ r i j := by
  dsimp only
  rw [MvPolynomial.integralInner_withDensity_polynomial _ _ _ hu (hpq := (T.hasFiniteMoments (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 (ae_mem_Icc_jacobiMeasure _ _))).integrable_mul _ _)]
  have he := integral_quadratic_deformed_eq_dunklXu_pencil hβ hδ hδ1 hδpos T
    (dunklXuEvenBlockCoordinates r i).1 (dunklXuEvenBlockCoordinates r i).2
    (dunklXuEvenBlockCoordinates r j).1 (dunklXuEvenBlockCoordinates r j).2 A B C
  rw [hAB, one_mul] at he
  rw [dunklXuEvenGramPencil_apply]
  have hcoord : dunklXuEvenBlockCoordinates r i = dunklXuEvenBlockCoordinates r j ↔ i = j :=
    (dunklXuEvenBlockCoordinates_injective r).eq_iff
  simpa only [integralInner, coordinateQuadratic, map_add, map_mul, map_pow, eval_C, eval_X,
    mul_assoc, Prod.mk.eta, hcoord, dunklXuEvenIndex] using he

/-- The distant factor entry survives extraction from a full factor with
interleaved orthogonal sectors. Exact sector coverage and ordering are explicit. -/
theorem eventually_dunklXu_full_factor_entry_ne_zero {β : ℝ} (hβ : -1 < β)
    (u v : ℝ) (hv : v ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ t in 𝓝[≠] (0 : ℝ), ∀ {N : ℕ} {γ : Type}
      (R : Matrix (Fin N) (Fin N) ℝ) (sector : Fin N → γ)
      (e : DunklXuEvenInterior r ⊕ Fin (r + 1) → Fin N) (c : γ),
      R.IsUpperTriangular → (∀ i, R i i ≠ 0) →
      (∀ i j, sector i ≠ sector j → (R.transpose * R) i j = 0) →
      Function.Injective e → (∀ i, sector (e i) = c) →
      (∀ k, sector k = c → k ∈ Set.range e) →
      (∀ i j, e (Sum.inl i) < e (Sum.inr j)) → StrictMono (e ∘ Sum.inr) →
      (R.transpose * R).submatrix e e = dunklXuEvenGramPencil hβ u v t r →
      R (e (Sum.inr 0)) (e (Sum.inr (Fin.last r))) ≠ 0 := by
  filter_upwards [eventually_dunklXuQuadratic_cholesky_first_last_ne_zero hβ u v hv r hr]
    with t ht
  intro N γ R sector e c hR hd hG he hc hcover hbefore hmono hP
  let A := R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inl)
  let B := R.submatrix (e ∘ Sum.inl) (e ∘ Sum.inr)
  let D := R.submatrix (e ∘ Sum.inr) (e ∘ Sum.inr)
  have hblocks := hR.gram_sector_blocks hd sector hG e he c hc hcover hbefore
  change (R.transpose * R).submatrix e e =
    fromBlocks (A.transpose * A) (A.transpose * B) (B.transpose * A)
      (B.transpose * B + D.transpose * D) at hblocks
  have h := hblocks.symm.trans hP
  have hAA : A.transpose * A = 1 + t • dunklXuInteriorMatrix hβ u v r := by
    ext i j
    exact congrArg (fun M ↦ M (Sum.inl i) (Sum.inl j)) h
  have hAB : A.transpose * B = t • dunklXuQuadraticRightBlock hβ u v r := by
    ext i j
    exact congrArg (fun M ↦ M (Sum.inl i) (Sum.inr j)) h
  have hBA : B.transpose * A = t • dunklXuQuadraticLeftBlock hβ u v r := by
    ext i j
    exact congrArg (fun M ↦ M (Sum.inr i) (Sum.inl j)) h
  have hDD : B.transpose * B + D.transpose * D =
      1 + t • dunklXuQuadraticBoundaryBlock hβ u v r := by
    ext i j
    exact congrArg (fun M ↦ M (Sum.inr i) (Sum.inr j)) h
  exact ht A B D (hR.isUnit_det_submatrix hd _ (he.comp Sum.inl_injective))
    (hR.submatrix_strictMono _ hmono) hAA hAB hBA hDD

/-- The actual canonical Cartesian system with the deformed Jacobi outer factor. -/
noncomputable def dunklXuDeformedCartesianSystem {β δ : ℝ} (hβ : -1 < β)
    (hδ : 0 < δ) (hδ1 : δ ≤ 1) :
    CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val)))
      (jacobiMeasure β β) := by
  have hβ' : -1 < β + 1 / 2 := by linarith
  let := isFiniteMeasure_jacobiMeasure hβ' hβ'
  let := isNegInvariant_jacobiMeasure β
  exact cartesianDiskSystem _ _
    (hasFiniteMoments_diskDeformedMeasure hδ.le hδ1 (ae_mem_Icc_jacobiMeasure _ _))
    (hasPositiveDefiniteIntegral_diskDeformed_jacobiMeasure hβ' hβ' hδ hδ1)
    (hasFiniteMoments_jacobiMeasure hβ hβ) (hasPositiveDefiniteIntegral_jacobiMeasure hβ hβ)
    (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 (ae_mem_Icc_jacobiMeasure _ _))

/-- The elementary three scalar inequalities give a nonnegative quadratic density. -/
theorem ae_coordinateQuadratic_nonneg_cartesianDiskMeasure {μ ν : Measure ℝ}
    [SFinite ν] (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    (hz : ∀ᵐ z ∂ν, z ∈ Icc (-1 : ℝ) 1) (A B C : ℝ)
    (hA : 0 ≤ A) (hAB : 0 ≤ A + B) (hAC : 0 ≤ A + C) :
    ∀ᵐ v ∂cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν), 0 ≤ eval v (coordinateQuadratic A B C) := by
  filter_upwards [ae_sq_add_sq_le_one_cartesianDiskMeasure hx hz] with v hv
  simpa [coordinateQuadratic] using
    (quadraticDisk_nonneg_iff A B C).2 ⟨hA, hAB, hAC⟩ (v 0) (v 1) hv

/-- Every admissible nonzero quadratic modification remains positive on nonzero polynomials. -/
theorem CartesianDiskSystem.hasPositiveDefiniteIntegral_coordinateQuadratic
    {μ ν : Measure ℝ} [IsFiniteMeasure μ] [IsFiniteMeasure ν]
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν) (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1)
    (hz : ∀ᵐ z ∂ν, z ∈ Icc (-1 : ℝ) 1) (A B C : ℝ)
    (hA : 0 ≤ A) (hAB : 0 < A + B) (hAC : 0 ≤ A + C) :
    MvPolynomial.HasPositiveDefiniteIntegral (MvPolynomial.polynomialDensityMeasure ((coordinateQuadratic A B C)) (((cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν))).comap Subtype.val)) :=
  (S.hasPositiveDefiniteIntegral hx hz).withDensity_polynomial
    (hasFiniteMoments_cartesianDiskMeasure hx hz) (coordinateQuadratic A B C)
    (coordinateQuadratic_ne_zero_of_boundary_ne_zero A B C hAB.ne')
    (ae_coordinateQuadratic_nonneg_cartesianDiskMeasure hx hz A B C hA hAB.le hAC)

/-- Scaling a nonnegative polynomial density scales every finite Gram matrix. -/
theorem integralGram_withDensity_smul (μ : Measure (σ → ℝ))
    (u : MvPolynomial σ ℝ) (P : Fin n → MvPolynomial σ ℝ)
    (hu : ∀ᵐ x ∂μ, 0 ≤ eval x u) {c : ℝ} (hc : 0 ≤ c)
    (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 ((c • u)) ((μ).comap Subtype.val)) P
      (hP := fun i j => (integrable_mul_withDensity_polynomial_iff (c • u) (P i) (P j)
        (hu.mono (fun x hx => by simpa using mul_nonneg hc hx))).2 (by
          simpa only [smul_eq_C_mul, map_mul, eval_C, mul_assoc] using (hi i j).const_mul c)) =
      c • 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)) := by
  have hcu : ∀ᵐ x ∂μ, 0 ≤ eval x (c • u) := by
    filter_upwards [hu] with x hx
    simpa using mul_nonneg hc hx
  have hci (i j) : Integrable (fun x => eval x ((c • u)*P i)*eval x (P j)) μ := by
    simpa only [smul_eq_C_mul, map_mul, eval_C, mul_assoc] using (hi i j).const_mul c
  rw [integralGram_withDensity_polynomial _ _ hcu (hi := hci),
    integralGram_withDensity_polynomial _ _ hu (hi := hi)]
  ext i j
  change integralInner μ (IntegrablePair.of μ ((c • u) * P i) (P j) (h := hci i j)) =
    c * integralInner μ (IntegrablePair.of μ (u * P i) (P j) (h := hi i j))
  simp only [integralInner, smul_eq_C_mul, map_mul, eval_C, mul_assoc, integral_const_mul]

/-- Positive scalar Gram scaling gives exactly positive square-root factor scaling. -/
theorem integralGramFactor_eq_sqrt_smul_of_gram_eq
    (μ ν : Measure (σ → ℝ)) (hμ : HasFiniteMoments μ) (hν : HasFiniteMoments ν)
    (hpμ : HasPositiveDefiniteIntegral μ) (hpν : HasPositiveDefiniteIntegral ν)
    {V : Submodule ℝ (MvPolynomial σ ℝ)} (b : Basis (Fin n) ℝ V)
    {c : ℝ} (hc : 0 < c)
    (hG : integralGram μ (fun i ↦ (b i : MvPolynomial σ ℝ)) =
      c • integralGram ν (fun i ↦ (b i : MvPolynomial σ ℝ))) :
    integralGramFactor μ hμ hpμ b = Real.sqrt c • integralGramFactor ν hν hpν b := by
  obtain ⟨hR, hd, hg⟩ := integralGramFactor_spec ν hν hpν b
  apply integralGramFactor_eq
  · intro i j hij
    change Real.sqrt c * integralGramFactor ν hν hpν b i j = 0
    rw [hR hij, mul_zero]
  · intro i
    exact mul_pos (Real.sqrt_pos.mpr hc) (hd i)
  · rw [hG, hg]
    rw [Matrix.transpose_smul, Matrix.smul_mul, Matrix.mul_smul, smul_smul,
      ← pow_two, Real.sq_sqrt hc.le]

/-- Positive Gram scaling does not change which factor entries vanish. -/
theorem integralGramFactor_eq_zero_iff_of_gram_scale
    (μ ν : Measure (σ → ℝ)) (hμ : HasFiniteMoments μ) (hν : HasFiniteMoments ν)
    (hpμ : HasPositiveDefiniteIntegral μ) (hpν : HasPositiveDefiniteIntegral ν)
    {V : Submodule ℝ (MvPolynomial σ ℝ)} (b : Basis (Fin n) ℝ V)
    {c : ℝ} (hc : 0 < c)
    (hG : integralGram μ (fun i ↦ (b i : MvPolynomial σ ℝ)) =
      c • integralGram ν (fun i ↦ (b i : MvPolynomial σ ℝ))) (i j : Fin n) :
    integralGramFactor μ hμ hpμ b i j = 0 ↔ integralGramFactor ν hν hpν b i j = 0 := by
  rw [integralGramFactor_eq_sqrt_smul_of_gram_eq μ ν hμ hν hpμ hpν b hc hG]
  change Real.sqrt c * integralGramFactor ν hν hpν b i j = 0 ↔ _
  simp [(Real.sqrt_pos.mpr hc).ne']

/-- Increasing enumeration of the complete degree/inner ordering. -/
noncomputable def diskDegreeOrder (N : ℕ) :
    Fin (Fintype.card (DiskDegreeIndex N)) ≃o DiskDegreeIndex N :=
  Fintype.orderIsoFinOfCardEq _ rfl

/-- The existing complete degree basis, reindexed in degree/inner order. -/
noncomputable def orderedDegreeBasis (N : ℕ) :
    Basis (Fin (Fintype.card (DiskDegreeIndex N))) ℝ (restrictTotalDegree (Fin 2) ℝ N) :=
  (S.degreeBasis hx hz N).reindex
    ((diskDegreeIndexEquiv N).symm.trans (diskDegreeOrder N).symm.toEquiv)

/-- Every even-even index occurs in this sector enumeration. -/
theorem dunklXuEvenSectorEmbedding_covers (r : ℕ)
    (i : Fin (Fintype.card (DiskDegreeIndex (2 * r))))
    (h0 : ((diskDegreeOrder (2 * r) i).1 0) % 2 = 0)
    (h1 : ((diskDegreeOrder (2 * r) i).1 1) % 2 = 0) :
    i ∈ Set.range (dunklXuEvenSectorEmbedding r) := by
  let d := diskDegreeOrder (2 * r) i
  have hd : d.1 0 + d.1 1 ≤ 2 * r := by
    simpa only [sum_exponents_fin_two] using d.2
  have h0' : d.1 0 % 2 = 0 := h0
  have h1' : d.1 1 % 2 = 0 := h1
  have ha : 2 * (d.1 0 / 2) = d.1 0 := by omega
  have hb : 2 * (d.1 1 / 2) = d.1 1 := by omega
  have hs : d.1 0 / 2 + d.1 1 / 2 ≤ r := by omega
  suffices ∃ j, dunklXuEvenBlockCoordinates r j = (d.1 0 / 2, d.1 1 / 2) by
    obtain ⟨j, hj⟩ := this
    refine ⟨j, (diskDegreeOrder (2 * r)).injective ?_⟩
    apply DiskDegreeIndex.ext
    rw [diskDegreeOrder_evenSectorEmbedding, hj]
    ext t
    fin_cases t
    · simpa [d] using ha
    · simpa [d] using hb
  by_cases he : d.1 0 / 2 + d.1 1 / 2 < r
  · exact ⟨Sum.inl ⟨(⟨d.1 0 / 2, by omega⟩, ⟨d.1 1 / 2, by omega⟩), he⟩, rfl⟩
  · refine ⟨Sum.inr ⟨d.1 1 / 2, by omega⟩, ?_⟩
    dsimp [dunklXuEvenBlockCoordinates]
    congr 1
    omega

/-- All interior indices precede every top-degree even index. -/
theorem dunklXuEvenSectorEmbedding_before (r : ℕ) (i : DunklXuEvenInterior r) (j : Fin (r + 1)) :
    dunklXuEvenSectorEmbedding r (Sum.inl i) < dunklXuEvenSectorEmbedding r (Sum.inr j) := by
  apply (diskDegreeOrder (2 * r)).lt_iff_lt.mp
  rw [diskDegreeIndex_lt_iff]
  left
  simp only [diskDegreeOrder_evenSectorEmbedding, Finsupp.add_apply]
  simp only [dunklXuEvenBlockCoordinates, Finsupp.single_eq_same,
    Finsupp.single_eq_of_ne (by decide : (0 : Fin 2) ≠ 1),
    Finsupp.single_eq_of_ne (by decide : (1 : Fin 2) ≠ 0), add_zero, zero_add]
  have := i.2
  omega

/-- The top block is ordered by increasing inner index. -/
theorem dunklXuEvenSectorEmbedding_strictMono (r : ℕ) :
    StrictMono (dunklXuEvenSectorEmbedding r ∘ Sum.inr) := by
  intro i j hij
  apply (diskDegreeOrder (2 * r)).lt_iff_lt.mp
  rw [diskDegreeIndex_lt_iff]
  right
  simp only [Function.comp_apply, diskDegreeOrder_evenSectorEmbedding, Finsupp.add_apply]
  simp only [dunklXuEvenBlockCoordinates, Finsupp.single_eq_same,
    Finsupp.single_eq_of_ne (by decide : (0 : Fin 2) ≠ 1),
    Finsupp.single_eq_of_ne (by decide : (1 : Fin 2) ≠ 0), add_zero, zero_add]
  constructor <;> omega

/-- Nonseparated admissible quadratics force distant canonical factor entries
for sufficiently small positive outer-factor deformations. -/
theorem eventually_integralGramFactor_admissible_quadratic_ne_zero {β : ℝ}
    (hβ : -1 < β) (A B C : ℝ) (hA : 0 ≤ A) (hAB : 0 < A + B) (hAC : 0 ≤ A + C) (hC : C ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ δ in 𝓝[≠] (0 : ℝ), ∀ (_hδ : 0 < δ) (_hδ1 : δ ≤ 1)
      (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β)),
      let μ := cartesianDiskMeasure
        (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β))
      let ω := MvPolynomial.polynomialDensityMeasure ((coordinateQuadratic A B C)) ((μ).comap Subtype.val)
      ∀ (hM : MvPolynomial.HasFiniteMoments ω) (hP : MvPolynomial.HasPositiveDefiniteIntegral ω)
        {N : ℕ} {V : Submodule ℝ (MvPolynomial (Fin 2) ℝ)} (b : Basis (Fin N) ℝ V)
        (n k : Fin N → ℕ)
        (_hb : ∀ i, (b i : MvPolynomial (Fin 2) ℝ) = T.polynomial (n i) (k i))
        (e : DunklXuEvenInterior r ⊕ Fin (r + 1) → Fin N),
        Function.Injective e →
        (∀ i, n (e i) = 2 * (dunklXuEvenBlockCoordinates r i).1 +
          2 * (dunklXuEvenBlockCoordinates r i).2 ∧
          k (e i) = 2 * (dunklXuEvenBlockCoordinates r i).2) →
        (∀ i, (n i - k i) % 2 = 0 → k i % 2 = 0 → i ∈ Set.range e) →
        (∀ i j, e (Sum.inl i) < e (Sum.inr j)) → StrictMono (e ∘ Sum.inr) →
        integralGramFactor ω hM hP b (e (Sum.inr 0)) (e (Sum.inr (Fin.last r))) ≠ 0 := by
  have hnorm : A / (A + B) + B / (A + B) = 1 := by
    rw [← add_div, div_self hAB.ne']
  filter_upwards [eventually_integralGramFactor_deformed_quadratic_ne_zero hβ
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hnorm
    (div_ne_zero hC hAB.ne') r hr] with δ ht
  intro hδ hδ1 T
  dsimp only
  intro hM hP N V b n k hb e he hlabels hcover hbefore hmono
  have hβ' : -1 < β + 1 / 2 := by linarith
  let := isFiniteMeasure_jacobiMeasure hβ' hβ'
  let := isFiniteMeasure_jacobiMeasure hβ hβ
  have hx := ae_mem_Icc_diskDeformedMeasure hδ.le hδ1
    (ae_mem_Icc_jacobiMeasure (β + 1 / 2) (β + 1 / 2))
  have hz := ae_mem_Icc_jacobiMeasure β β
  have hA' := div_nonneg hA hAB.le
  have hAC' : 0 ≤ A / (A + B) + C / (A + B) := by
    rw [← add_div]
    exact div_nonneg hAC hAB.le
  have hu := ae_coordinateQuadratic_nonneg_cartesianDiskMeasure hx hz
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hA' (by rw [hnorm]; norm_num) hAC'
  have hm := (hasFiniteMoments_cartesianDiskMeasure hx hz).withDensity_polynomial _ hu
  have hp := T.hasPositiveDefiniteIntegral_coordinateQuadratic hx hz
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hA' (by rw [hnorm]; norm_num) hAC'
  have hn := ht hδ hδ1 (hasPositiveDefiniteIntegral_diskDeformed_jacobiMeasure
    hβ' hβ' hδ hδ1) T hu hm hp b n k hb e he hlabels hcover hbefore hmono
  have hscale : coordinateQuadratic A B C =
      (A + B) • coordinateQuadratic (A / (A + B)) (B / (A + B)) (C / (A + B)) := by
    simp only [coordinateQuadratic, smul_eq_C_mul, mul_add, ← mul_assoc, ← map_mul]
    field_simp
  intro hzR
  apply hn
  apply (integralGramFactor_eq_zero_iff_of_gram_scale _ _ hM hm hP hp b hAB ?_
    (e (Sum.inr 0)) (e (Sum.inr (Fin.last r)))).mp hzR
  · have hs := integralGram_withDensity_smul _ _
      (fun i => (b i : MvPolynomial (Fin 2) ℝ)) hu hAB.le
      (hi := fun _ _ => (T.hasFiniteMoments hx).integrable_mul _ _)
    ext i j
    have hij := congrFun (congrFun hs i) j
    dsimp only [integralGram, integralInner, IntegrablePair.of] at hij ⊢
    rw [hscale]
    exact hij

/-- The certified measure is determined by its nonnegative density on measurable sets. -/
theorem dunklXuDeformedQuadraticMeasure_apply (β δ A B C : ℝ)
    (d : DunklXuDeformedQuadraticInput β δ A B C)
    {s : Set (Fin 2 → ℝ)} (hs : MeasurableSet s) :
    dunklXuDeformedQuadraticMeasure β δ A B C d s =
      ∫⁻ x : {x : Fin 2 → ℝ // 0 ≤ eval x (coordinateQuadratic A B C)} in Subtype.val ⁻¹' s,
        ENNReal.ofNNReal (NNReal.mk (eval (x : Fin 2 → ℝ) (coordinateQuadratic A B C)) x.property)
        ∂((cartesianDiskMeasure
    (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β))).comap Subtype.val) :=
  MvPolynomial.polynomialDensityMeasure_apply _ _ hs

/-- Complete source degree basis in the fixed degree/inner ordering. -/
noncomputable def dunklXuDeformedDegreeBasis {β δ : ℝ} (hβ : -1 < β)
    (hδ : 0 < δ) (hδ1 : δ ≤ 1) (N : ℕ) :
    Basis (Fin (Fintype.card (DiskDegreeIndex N))) ℝ (restrictTotalDegree (Fin 2) ℝ N) := by
  have hβ' : -1 < β + 1 / 2 := by linarith
  let := isFiniteMeasure_jacobiMeasure hβ' hβ'
  let := isFiniteMeasure_jacobiMeasure hβ hβ
  exact (dunklXuDeformedCartesianSystem hβ hδ hδ1).orderedDegreeBasis
    (ae_mem_Icc_diskDeformedMeasure hδ.le hδ1 (ae_mem_Icc_jacobiMeasure _ _))
    (ae_mem_Icc_jacobiMeasure _ _) N

/-- Each requested degree admits an actual factor with inner gap 2r.
All bases, measures, positivity witnesses and sector embeddings are constructed. -/
theorem exists_integralGramFactor_quadratic_distant_entry {β : ℝ} (hβ : -1 < β)
    (A B C : ℝ) (hA : 0 ≤ A) (hAB : 0 < A + B) (hAC : 0 ≤ A + C) (hC : C ≠ 0)
    (r : ℕ) (hr : 2 ≤ r) :
    ∃ (δ : ℝ) (hδ : 0 < δ) (hδ1 : δ ≤ 1),
      ∃ (hM : MvPolynomial.HasFiniteMoments (dunklXuDeformedQuadraticMeasure β δ A B C))
        (hP : MvPolynomial.HasPositiveDefiniteIntegral (dunklXuDeformedQuadraticMeasure β δ A B C)),
      integralGramFactor (dunklXuDeformedQuadraticMeasure β δ A B C) hM hP
        (dunklXuDeformedDegreeBasis hβ hδ hδ1 (2 * r))
        (dunklXuEvenSectorEmbedding r (Sum.inr 0))
        (dunklXuEvenSectorEmbedding r (Sum.inr (Fin.last r))) ≠ 0 := by
  have he := eventually_integralGramFactor_admissible_quadratic_ne_zero
    hβ A B C hA hAB hAC hC r hr
  obtain ⟨ε, hε, he⟩ := Metric.mem_nhdsWithin_iff.mp he
  let δ := min ε 1 / 2
  have hδ : 0 < δ := div_pos (lt_min hε zero_lt_one) (by norm_num)
  have hδ1 : δ ≤ 1 := by dsimp [δ]; linarith [min_le_right ε 1]
  have hδε : δ < ε := by dsimp [δ]; linarith [min_le_left ε 1]
  have ht := he (show δ ∈ Metric.ball 0 ε ∩ ({0}ᶜ : Set ℝ) from
    ⟨by simpa [Metric.mem_ball, Real.dist_eq, abs_of_pos hδ] using hδε,
      by simpa using hδ.ne'⟩)
  obtain ⟨hM, hP⟩ := dunklXuDeformedQuadraticMeasure_admissible hβ hδ hδ1 A B C hA hAB hAC
  refine ⟨δ, hδ, hδ1, hM, hP, ?_⟩
  simp only [dunklXuDeformedQuadraticMeasure] at hM hP ⊢
  refine ht hδ hδ1 (dunklXuDeformedCartesianSystem hβ hδ hδ1) hM hP
    (dunklXuDeformedDegreeBasis hβ hδ hδ1 (2 * r))
    (fun i ↦ (diskDegreeOrder (2 * r) i).1 0 + (diskDegreeOrder (2 * r) i).1 1)
    (fun i ↦ (diskDegreeOrder (2 * r) i).1 1)
    (fun i ↦ dunklXuDeformedDegreeBasis_apply hβ hδ hδ1 (2 * r) i)
    (dunklXuEvenSectorEmbedding r) (dunklXuEvenSectorEmbedding_injective r) ?_ ?_
    (dunklXuEvenSectorEmbedding_before r) (dunklXuEvenSectorEmbedding_strictMono r)
  · intro i
    simp [diskDegreeOrder_evenSectorEmbedding]
  · intro i h0 h1
    apply dunklXuEvenSectorEmbedding_covers r i _ h1
    simpa using h0

/-- The complete finite connection in the fixed degree/inner ordering. -/
noncomputable def orderedConnectionMatrix (N : ℕ) :
    Matrix (Fin (Fintype.card (DiskDegreeIndex N))) (Fin (Fintype.card (DiskDegreeIndex N))) ℝ :=
  (T.orderedDegreeBasis hx' hz' N).toMatrix (S.orderedDegreeBasis hx hz N)

/-- The separated connection is upper triangular in total-degree/inner order. -/
theorem orderedConnectionMatrix_isUpperTriangular (N : ℕ) :
    (S.orderedConnectionMatrix U hx hz (ae_mem_diskOuterModification hx h l) hzU N).IsUpperTriangular := by
  intro i j hij
  rw [orderedConnectionMatrix_apply]
  apply S.connectionEntry_eq_zero_of_outside_band h g l U hμ hν hx hh hg hgpar hgdeg
    (by omega) (by omega)
  rintro ⟨s, r, hs, hr, hn, hk⟩
  have hord := (diskDegreeOrder N).strictMono hij
  rw [diskDegreeIndex_lt_iff] at hord
  dsimp only [id_eq] at hord
  rcases hord with hord | ⟨hord, hi⟩ <;> omega

/-- Its diagonal is strictly positive, with no independent sign choices. -/
theorem orderedConnectionMatrix_diag_pos (N : ℕ)
    (i : Fin (Fintype.card (DiskDegreeIndex N))) :
    0 < S.orderedConnectionMatrix U hx hz (ae_mem_diskOuterModification hx h l) hzU N i i := by
  rw [orderedConnectionMatrix_apply]
  exact S.connectionEntry_diag_pos h g l U hμ hν hx hh hg _ _

/-- Matching polynomial pairings identify the separated connection with the
canonical factor, including for an independently defined modified measure. -/
theorem integralGramFactor_eq_orderedConnectionMatrix
    (ω : Measure (Fin 2 → ℝ)) (hω : MvPolynomial.HasFiniteMoments ω)
    (hpω : MvPolynomial.HasPositiveDefiniteIntegral ω)
    (hpair : ∀ p q : MvPolynomial (Fin 2) ℝ, integralInner ω (IntegrablePair.of ω p q) =
      integralInner (cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskOuterModification h l (μ.comap Subtype.val)))
        (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskOuterModification h l (μ.comap Subtype.val))) (Polynomial.polynomialDensityMeasure (g) ((ν).comap Subtype.val)))) (IntegrablePair.of (cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskOuterModification h l (μ.comap Subtype.val)))
        (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskOuterModification h l (μ.comap Subtype.val))) (Polynomial.polynomialDensityMeasure (g) ((ν).comap Subtype.val)))) p q)) (N : ℕ) :
    integralGramFactor ω hω hpω (S.orderedDegreeBasis hx hz N) =
      S.orderedConnectionMatrix U hx hz (ae_mem_diskOuterModification hx h l) hzU N := by
  apply integralGramFactor_eq
  · exact S.orderedConnectionMatrix_isUpperTriangular hx hz h g l U hμ hν hh hg hgpar hgdeg hzU N
  · exact S.orderedConnectionMatrix_diag_pos hx hz h g l U hμ hν hh hg hzU N
  · rw [← S.integralGram_eq_orderedConnectionMatrix_transpose_mul U hx hz
      (ae_mem_diskOuterModification hx h l) hzU N]
    ext i j
    exact hpair _ _

/-- Separated polynomial modification bounds both directions of the canonical
factor's inner-index gap, uniformly over degree sections. -/
theorem integralGramFactor_eq_zero_of_separated_inner_gap
    (S : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν)
    (hμ : Polynomial.HasFiniteMoments μ) (hpμ : Polynomial.HasPositiveDefiniteIntegral μ)
    (hν : Polynomial.HasFiniteMoments ν) (hpν : Polynomial.HasPositiveDefiniteIntegral ν)
    (hx : ∀ᵐ x ∂μ, x ∈ Icc (-1 : ℝ) 1) (hz : ∀ᵐ x ∂ν, x ∈ Icc (-1 : ℝ) 1)
    (h g : ℝ[X]) (l : ℕ) (hh0 : h ≠ 0) (hg0 : g ≠ 0)
    (hh : ∀ᵐ x ∂μ, 0 ≤ h.eval x) (hg : ∀ᵐ x ∂ν, 0 ≤ g.eval x)
    (hgpar : ∀ x, g.eval (-x) = g.eval x)
    (hgcoeff : ∀ i, i % 2 ≠ 0 → g.coeff i = 0) (hgdeg : g.natDegree ≤ 2 * l)
    (u : MvPolynomial (Fin 2) ℝ) (hu : ∀ᵐ x ∂cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν), 0 ≤ eval x u)
    (hmod : u = cartesianDiskModification h g l)
    (hM : MvPolynomial.HasFiniteMoments
      (MvPolynomial.polynomialDensityMeasure (u) (((cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν))).comap Subtype.val)))
    (hP : MvPolynomial.HasPositiveDefiniteIntegral
      (MvPolynomial.polynomialDensityMeasure (u) (((cartesianDiskMeasure (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) (ProductMeasureInput.ofSFinite (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) μ) ν))).comap Subtype.val)))
    (N : ℕ) (i j : Fin (Fintype.card (DiskDegreeIndex N)))
    (hgap : (diskDegreeOrder N i).1 1 + g.natDegree < (diskDegreeOrder N j).1 1 ∨
      (diskDegreeOrder N j).1 1 + g.natDegree < (diskDegreeOrder N i).1 1) :
    integralGramFactor _ hM hP (S.orderedDegreeBasis hx hz N) i j = 0 := by
  let := (hμ.diskOuterModification hx h hh l).isFiniteMeasure
  let := (hν.polynomialDensityMeasure g hg).isFiniteMeasure
  have U := modifiedCartesianDiskSystem hμ hpμ hν hpν hx h g l hh0 hg0 hh hg hgpar
  have hzU : ∀ᵐ x ∂Polynomial.polynomialDensityMeasure (g) ((ν).comap Subtype.val),
      x ∈ Icc (-1 : ℝ) 1 := (polynomialDensityMeasure_comap_absolutelyContinuous g ν).ae_le hz
  rw [S.integralGramFactor_eq_orderedConnectionMatrix hx hz h g l U hμ hν hh hg
    hgpar hgdeg hzU _ hM hP (fun p q ↦ ?_) N, orderedConnectionMatrix_apply]
  · apply S.connectionEntry_eq_zero_of_inner_band h g l U hμ hν hx hh hg
    omega
  · rw [MvPolynomial.integralInner_withDensity_polynomial _ _ _ hu
      (hpq := (S.hasFiniteMoments hx).integrable_mul _ _)]
    simp only [hmod]
    exact (integralInner_cartesianDiskModification hμ hν hx hz h g l hh hg hgdeg hgcoeff p q).symm

/-- Evaluation of the numerator is the actual adjugate sandwich. -/
theorem eval_inversePencilNumerator (K : Matrix ι ι ℝ)
    (F : Matrix κ ι ℝ) (E : Matrix ι κ ℝ) (i j : κ) (t : ℝ) :
    (inversePencilNumerator K F E i j).eval t = (F * (1 + t • K).adjugate * E) i j := by
  change ((F.map Polynomial.C * (polynomialPencil K).adjugate * E.map Polynomial.C).map
    (Polynomial.evalRingHom t)) i j = _
  rw [Matrix.map_mul, Matrix.map_mul]
  change ((F.map Polynomial.C).map (Polynomial.evalRingHom t) *
    (Polynomial.evalRingHom t).mapMatrix (polynomialPencil K).adjugate *
    (E.map Polynomial.C).map (Polynomial.evalRingHom t)) i j = _
  rw [RingHom.map_adjugate, eval_polynomialPencil]
  simp [Matrix.map_map, Function.comp_def]

/-- Only finitely many real parameters are singular. -/
theorem eventually_cofinite_isUnit_pencil (K : Matrix ι ι ℝ) :
    ∀ᶠ t : ℝ in Filter.cofinite, IsUnit (1 + t • K) := by
  have hd : (polynomialPencil K).det ≠ 0 := by
    intro hz
    have h := congrArg (fun p : ℝ[X] ↦ p.eval 0) hz
    have he := ((Polynomial.evalRingHom (0 : ℝ)).map_det (polynomialPencil K)).trans
      (congrArg Matrix.det (eval_polynomialPencil K 0))
    simp only [zero_smul, add_zero, Matrix.det_one] at he
    change (polynomialPencil K).det.eval 0 = 1 at he
    simp [he] at h
  filter_upwards [Polynomial.eventually_cofinite_not_isRoot hd] with t ht
  exact (isUnit_pencil_iff_det_ne_zero K t).2 ht

/-- A nonzero classical inverse sandwich stays nonzero outside a finite exceptional set. -/
theorem eventually_cofinite_inversePencil_sandwich_ne_zero (K : Matrix ι ι ℝ)
    (F : Matrix κ ι ℝ) (E : Matrix ι κ ℝ) (i j : κ)
    (hn : ∃ t : InvertiblePencilParameter K, (F * inversePencil K t * E) i j ≠ 0) :
    ∀ᶠ t : ℝ in Filter.cofinite, ∃ ht : IsUnit (1 + t • K),
      (F * inversePencil K ⟨t, ht⟩ * E) i j ≠ 0 := by
  have hp : inversePencilNumerator K F E i j ≠ 0 := by
    obtain ⟨t, ht⟩ := hn
    intro hz
    apply ht
    rw [inversePencil_sandwich_eq, hz]
    simp
  filter_upwards [eventually_cofinite_isUnit_pencil K,
    Polynomial.eventually_cofinite_not_isRoot hp] with t ht hp
  refine ⟨ht, ?_⟩
  rw [inversePencil_sandwich_eq]
  exact mul_ne_zero (inv_ne_zero ((isUnit_pencil_iff_det_ne_zero K t).1 ht)) hp

/-- Only finitely many parameters can fail the full-factor implication. -/
theorem eventually_cofinite_dunklXu_full_factor_entry_ne_zero {β : ℝ} (hβ : -1 < β)
    (u v : ℝ) (hv : v ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ t in Filter.cofinite, ∀ {N : ℕ} {γ : Type}
      (R : Matrix (Fin N) (Fin N) ℝ) (sector : Fin N → γ)
      (e : DunklXuEvenInterior r ⊕ Fin (r + 1) → Fin N) (c : γ),
      R.IsUpperTriangular → (∀ i, R i i ≠ 0) →
      (∀ i j, sector i ≠ sector j → (R.transpose * R) i j = 0) →
      Function.Injective e → (∀ i, sector (e i) = c) →
      (∀ k, sector k = c → k ∈ Set.range e) →
      (∀ i j, e (Sum.inl i) < e (Sum.inr j)) → StrictMono (e ∘ Sum.inr) →
      (R.transpose * R).submatrix e e = dunklXuEvenGramPencil hβ u v t r →
      R (e (Sum.inr 0)) (e (Sum.inr (Fin.last r))) ≠ 0 := by
  have hn := eventually_dunklXuInterior_inverse_ne_zero_of_vertical_ne_zero hβ u v hv r hr
  have he : ∃ t : ℝ, (dunklXuQuadraticLeftBlock hβ u v r *
      (1 + t • dunklXuInteriorMatrix hβ u v r)⁻¹ * dunklXuQuadraticRightBlock hβ u v r)
      0 (Fin.last r) ≠ 0 := by
    simpa only [dunklXuInteriorFunctional_eq_sandwich] using hn.exists
  have he' : ∃ t : Matrix.InvertiblePencilParameter (dunklXuInteriorMatrix hβ u v r),
      (dunklXuQuadraticLeftBlock hβ u v r *
        Matrix.inversePencil (dunklXuInteriorMatrix hβ u v r) t *
        dunklXuQuadraticRightBlock hβ u v r) 0 (Fin.last r) ≠ 0 := by
    obtain ⟨t, ht⟩ := he
    have hu : IsUnit (1 + t • dunklXuInteriorMatrix hβ u v r) := by
      by_contra hu
      apply ht
      rw [Matrix.nonsing_inv_eq_ringInverse, Ring.inverse_non_unit _ hu]
      simp
    exact ⟨⟨t, hu⟩, by simpa only [Matrix.inversePencil_eq_inv] using ht⟩
  have hc := Matrix.eventually_cofinite_inversePencil_sandwich_ne_zero
    (dunklXuInteriorMatrix hβ u v r) (dunklXuQuadraticLeftBlock hβ u v r)
    (dunklXuQuadraticRightBlock hβ u v r) 0 (Fin.last r) he'
  filter_upwards [hc, Filter.eventually_cofinite_ne (0 : ℝ)] with t ht ht0
  obtain ⟨hu, ht⟩ := ht
  exact dunklXu_full_factor_entry_ne_zero_of_inverse hβ u v r hr ⟨t, hu⟩ ht0
    (by simpa only [dunklXuInteriorFunctional_eq_sandwich] using ht)

/-- Nonseparated admissible quadratics force distant canonical factor entries
outside a finite set of outer-factor deformation parameters. -/
theorem eventually_cofinite_integralGramFactor_admissible_quadratic_ne_zero {β : ℝ}
    (hβ : -1 < β) (A B C : ℝ) (hA : 0 ≤ A) (hAB : 0 < A + B) (hAC : 0 ≤ A + C) (hC : C ≠ 0) (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ δ in Filter.cofinite, ∀ (_hδ : 0 < δ) (_hδ1 : δ ≤ 1)
      (T : CartesianDiskSystem (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β)),
      let μ := cartesianDiskMeasure
        (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (ProductMeasureInput.ofSFinite (hν := by infer_instance) (Measure.comap (Subtype.val : Set.Icc (-1 : ℝ) 1 → ℝ) (diskDeformedMeasure δ ((jacobiMeasure (β + 1 / 2) (β + 1 / 2)).comap Subtype.val))) (jacobiMeasure β β))
      let ω := MvPolynomial.polynomialDensityMeasure ((coordinateQuadratic A B C)) ((μ).comap Subtype.val)
      ∀ (hM : MvPolynomial.HasFiniteMoments ω) (hP : MvPolynomial.HasPositiveDefiniteIntegral ω)
        {N : ℕ} {V : Submodule ℝ (MvPolynomial (Fin 2) ℝ)} (b : Basis (Fin N) ℝ V)
        (n k : Fin N → ℕ)
        (_hb : ∀ i, (b i : MvPolynomial (Fin 2) ℝ) = T.polynomial (n i) (k i))
        (e : DunklXuEvenInterior r ⊕ Fin (r + 1) → Fin N),
        Function.Injective e →
        (∀ i, n (e i) = 2 * (dunklXuEvenBlockCoordinates r i).1 +
          2 * (dunklXuEvenBlockCoordinates r i).2 ∧
          k (e i) = 2 * (dunklXuEvenBlockCoordinates r i).2) →
        (∀ i, (n i - k i) % 2 = 0 → k i % 2 = 0 → i ∈ Set.range e) →
        (∀ i j, e (Sum.inl i) < e (Sum.inr j)) → StrictMono (e ∘ Sum.inr) →
        integralGramFactor ω hM hP b (e (Sum.inr 0)) (e (Sum.inr (Fin.last r))) ≠ 0 := by
  have hnorm : A / (A + B) + B / (A + B) = 1 := by
    rw [← add_div, div_self hAB.ne']
  filter_upwards [eventually_cofinite_integralGramFactor_deformed_quadratic_ne_zero hβ
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hnorm
    (div_ne_zero hC hAB.ne') r hr] with δ ht
  intro hδ hδ1 T
  dsimp only
  intro hM hP N V b n k hb e he hlabels hcover hbefore hmono
  have hβ' : -1 < β + 1 / 2 := by linarith
  let := isFiniteMeasure_jacobiMeasure hβ' hβ'
  let := isFiniteMeasure_jacobiMeasure hβ hβ
  have hx := ae_mem_Icc_diskDeformedMeasure hδ.le hδ1
    (ae_mem_Icc_jacobiMeasure (β + 1 / 2) (β + 1 / 2))
  have hz := ae_mem_Icc_jacobiMeasure β β
  have hA' := div_nonneg hA hAB.le
  have hAC' : 0 ≤ A / (A + B) + C / (A + B) := by
    rw [← add_div]
    exact div_nonneg hAC hAB.le
  have hu := ae_coordinateQuadratic_nonneg_cartesianDiskMeasure hx hz
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hA' (by rw [hnorm]; norm_num) hAC'
  have hm := (hasFiniteMoments_cartesianDiskMeasure hx hz).withDensity_polynomial _ hu
  have hp := T.hasPositiveDefiniteIntegral_coordinateQuadratic hx hz
    (A / (A + B)) (B / (A + B)) (C / (A + B)) hA' (by rw [hnorm]; norm_num) hAC'
  have hn := ht hδ hδ1 (hasPositiveDefiniteIntegral_diskDeformed_jacobiMeasure
    hβ' hβ' hδ hδ1) T hu hm hp b n k hb e he hlabels hcover hbefore hmono
  have hscale : coordinateQuadratic A B C =
      (A + B) • coordinateQuadratic (A / (A + B)) (B / (A + B)) (C / (A + B)) := by
    simp only [coordinateQuadratic, smul_eq_C_mul, mul_add, ← mul_assoc, ← map_mul]
    field_simp
  intro hzR
  apply hn
  apply (integralGramFactor_eq_zero_iff_of_gram_scale _ _ hM hm hP hp b hAB ?_
    (e (Sum.inr 0)) (e (Sum.inr (Fin.last r)))).mp hzR
  · have hs := integralGram_withDensity_smul _ _
      (fun i => (b i : MvPolynomial (Fin 2) ℝ)) hu hAB.le
      (hi := fun _ _ => (T.hasFiniteMoments hx).integrable_mul _ _)
    ext i j
    have hij := congrFun (congrFun hs i) j
    dsimp only [integralGram, integralInner, IntegrablePair.of] at hij ⊢
    rw [hscale]
    exact hij

/-- Each complete degree section has only finitely many exceptional deformations. -/
theorem eventually_cofinite_integralGramFactor_quadratic_distant_entry {β : ℝ} (hβ : -1 < β)
    (A B C : ℝ) (hA : 0 ≤ A) (hAB : 0 < A + B) (hAC : 0 ≤ A + C) (hC : C ≠ 0)
    (r : ℕ) (hr : 2 ≤ r) :
    ∀ᶠ δ : ℝ in Filter.cofinite, ∀ (hδ : 0 < δ) (hδ1 : δ ≤ 1),
      ∃ (hM : MvPolynomial.HasFiniteMoments (dunklXuDeformedQuadraticMeasure β δ A B C))
        (hP : MvPolynomial.HasPositiveDefiniteIntegral (dunklXuDeformedQuadraticMeasure β δ A B C)),
      integralGramFactor (dunklXuDeformedQuadraticMeasure β δ A B C) hM hP
        (dunklXuDeformedDegreeBasis hβ hδ hδ1 (2 * r))
        (dunklXuEvenSectorEmbedding r (Sum.inr 0))
        (dunklXuEvenSectorEmbedding r (Sum.inr (Fin.last r))) ≠ 0 := by
  filter_upwards [eventually_cofinite_integralGramFactor_admissible_quadratic_ne_zero
    hβ A B C hA hAB hAC hC r hr] with δ ht
  intro hδ hδ1
  obtain ⟨hM, hP⟩ := dunklXuDeformedQuadraticMeasure_admissible hβ hδ hδ1 A B C hA hAB hAC
  refine ⟨hM, hP, ?_⟩
  simp only [dunklXuDeformedQuadraticMeasure] at hM hP ⊢
  refine ht hδ hδ1 (dunklXuDeformedCartesianSystem hβ hδ hδ1) hM hP
    (dunklXuDeformedDegreeBasis hβ hδ hδ1 (2 * r))
    (fun i ↦ (diskDegreeOrder (2 * r) i).1 0 + (diskDegreeOrder (2 * r) i).1 1)
    (fun i ↦ (diskDegreeOrder (2 * r) i).1 1)
    (fun i ↦ dunklXuDeformedDegreeBasis_apply hβ hδ hδ1 (2 * r) i)
    (dunklXuEvenSectorEmbedding r) (dunklXuEvenSectorEmbedding_injective r) ?_ ?_
    (dunklXuEvenSectorEmbedding_before r) (dunklXuEvenSectorEmbedding_strictMono r)
  · intro i
    simp [diskDegreeOrder_evenSectorEmbedding]
  · intro i h0 h1
    apply dunklXuEvenSectorEmbedding_covers r i _ h1
    simpa using h0