/-- Mathlib's first-kind Chebyshev sequence satisfies the general orthogonality
predicate for its standard weight (DLMF Table 18.3.1). -/
theorem isOrthogonal_chebyshevTsequence :
    (chebyshevTsequence ℝ).IsOrthogonal classicalMeasureT where
  integrable m n := integrable_classicalMeasureT (by fun_prop)
  orthogonal m n hmn := by simpa only [chebyshevTsequence, hmn, ite_false] using (classicalMeasureT_orthogonality m n).2

/-- Finite polynomial moments make all polynomial products integrable
(DLMF 18.2(i), Szegő, Chapter 3). -/
@[aesop (rule_sets := [PolynomialIntegrability]) unsafe 80% apply]
theorem HasFiniteMoments.integrable_mul {μ : Measure ℝ} (hμ : HasFiniteMoments μ)
    (p q : ℝ[X]) : Integrable (fun x => p.eval x * q.eval x) μ := by
  have hp : MemLp (fun x => p.eval x) 2 μ :=
    (memLp_two_iff_integrable_sq (by fun_prop)).2 (hμ p)
  have hq : MemLp (fun x => q.eval x) 2 μ :=
    (memLp_two_iff_integrable_sq (by fun_prop)).2 (hμ q)
  exact hp.integrable_mul hq

/-- Finite interval support supplies every polynomial moment. -/
@[aesop (rule_sets := [PolynomialIntegrability]) unsafe 70% apply]
theorem hasFiniteMoments_of_ae_mem_Icc {μ : Measure ℝ} [IsFiniteMeasure μ]
    {a b : ℝ} (hx : ∀ᵐ x ∂μ, x ∈ Icc a b) : HasFiniteMoments μ := by
  intro p
  have hi := (p.continuous.pow 2).integrableOn_Icc (μ := μ) (a := a) (b := b)
  rwa [IntegrableOn, Measure.restrict_eq_self_of_ae_mem hx] at hi

/-- Every measure with finite polynomial moments and positive squared integrals
admits an orthogonal polynomial sequence of exact degrees (DLMF 18.2(i), Szegő, Chapter 3). -/
theorem exists_isOrthogonal (μ : Measure ℝ) (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ) : ∃ P : Sequence ℝ, P.IsOrthogonal μ := by
  let core : InnerProductSpace.Core ℝ ℝ[X] := integralInnerProductCore μ hμ hpos
  let : NormedAddCommGroup ℝ[X] := InnerProductSpace.Core.toNormedAddCommGroup (𝕜 := ℝ)
  let : InnerProductSpace ℝ ℝ[X] :=
    InnerProductSpace.ofCore (inferInstance : PreInnerProductSpace.Core ℝ ℝ[X])
  let P : Sequence ℝ := ⟨fun n => X ^ n, fun n => by simp⟩
  have hdegree (n : ℕ) : (gramSchmidt ℝ P n).degree = n := by
    have hunit (i : ℕ) : IsUnit (P i).leadingCoeff :=
      P.isUnit_leadingCoeff i
    have hspan (k : ℕ) : span ℝ (P '' Set.Iio k) = degreeLT ℝ k :=
      P.span_degreeLT (fun i _ => hunit i)
    have hmem : gramSchmidt ℝ P n ∈ span ℝ (P '' Set.Iio (n + 1)) := by
      rw [← span_gramSchmidt_Iio ℝ (P : ℕ → ℝ[X]) (n + 1)]
      exact subset_span ⟨n, Nat.lt_succ_self n, rfl⟩
    rw [hspan, mem_degreeLT] at hmem
    have hnot : ¬ (gramSchmidt ℝ P n).degree < n := by
      intro h
      have hm : gramSchmidt ℝ P n ∈ span ℝ (P '' Set.Iio n) := by
        rw [hspan, mem_degreeLT]; exact h
      rw [← span_gramSchmidt_Iio ℝ (P : ℕ → ℝ[X]) n] at hm
      exact (gramSchmidt_linearIndependent P.linearIndependent).notMem_span_image
        (by simp : n ∉ Set.Iio n) hm
    have hn := gramSchmidt_ne_zero n P.linearIndependent
    rw [degree_eq_natDegree hn] at hmem hnot ⊢
    norm_cast at hmem hnot ⊢
    omega
  let Q : Sequence ℝ := ⟨gramSchmidt ℝ P, hdegree⟩
  refine ⟨Q, ⟨fun m n => hμ.integrable_mul (Q m) (Q n), ?_⟩⟩
  intro m n hmn
  exact gramSchmidt_orthogonal ℝ P hmn

/-- At each index, two orthogonal sequences differ by the ratio of their leading
coefficients (DLMF 18.2(i), Szegő, Chapter 3). -/
theorem IsOrthogonal.eq_leadingCoeff_div_smul {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
    {P Q : Sequence ℝ} (hP : P.IsOrthogonal μ) (hQ : Q.IsOrthogonal μ) (n : ℕ) :
    Q n = ((Q n).leadingCoeff / (P n).leadingCoeff) • P n := by
  let c : ℝ := (Q n).leadingCoeff / (P n).leadingCoeff
  have hP0 := leadingCoeff_ne_zero.mpr (P.ne_zero n)
  have hQ0 := leadingCoeff_ne_zero.mpr (Q.ne_zero n)
  have hc : c ≠ 0 := div_ne_zero hQ0 hP0
  have hdeg : (Q n - c • P n).degree < n := by
    have h := degree_sub_lt_left
      (show (Q n).degree = (c • P n).degree by
        simp [degree_smul_of_smul_regular _ (IsSMulRegular.of_ne_zero hc)])
      (Q.ne_zero n) (show (Q n).leadingCoeff = (c • P n).leadingCoeff by
        rw [leadingCoeff_smul_of_smul_regular _ (IsSMulRegular.of_ne_zero hc)]
        simp [c, hP0])
    simpa using h
  have hq := hQ.integralInner_eq_zero_of_degree_lt hμ hdeg
  have hp := hP.integralInner_eq_zero_of_degree_lt hμ hdeg
  have hzero : integralInner μ (IntegrablePair.of μ (Q n - c • P n) (Q n - c • P n)) = 0 := by
    rw [integralInner_sub_left hμ, integralInner_smul_left,
      integralInner_comm μ (Q n), integralInner_comm μ (P n), hq, hp]
    all_goals first | exact hμ.integrable_mul _ _ | ring
  exact sub_eq_zero.mp (hpos.eq_zero_of_integralInner_self_eq_zero hzero)

/-- A sequence element is orthogonal to every polynomial of lower degree
(DLMF 18.2(i), Szegő, Chapter 3). -/
theorem IsOrthogonal.integralInner_eq_zero_of_degree_lt {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) {P : Sequence ℝ} (hP : P.IsOrthogonal μ)
    {p : ℝ[X]} {n : ℕ} (hp : p.degree < n) : integralInner μ (IntegrablePair.of μ p (P n)) = 0 := by
  have hm : p ∈ span ℝ (P '' Set.Iio n) := by
    rw [P.span_degreeLT (fun i _ =>
      P.isUnit_leadingCoeff i), mem_degreeLT]
    exact hp
  clear hp
  induction hm using Submodule.span_induction with
  | mem q hq =>
    obtain ⟨m, hm, rfl⟩ := hq
    exact hP.orthogonal m n (ne_of_lt hm)
  | zero => exact integralInner_zero_left μ (P n)
  | add p q _ _ hp hq => rw [integralInner_add_left hμ, hp, hq, add_zero]
  | smul a p _ hp =>
    rw [integralInner_smul_left (h := hμ.integrable_mul p (P n)), hp, mul_zero]

/-- Coordinates in an orthogonal polynomial basis are integral projections
(DLMF 18.2(i), Szegő, Chapter 3). -/
theorem IsOrthogonal.integralInner_eq_basis_repr_mul {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) {P : Sequence ℝ} (hP : P.IsOrthogonal μ)
    (hu : ∀ i, IsUnit (P i).leadingCoeff) (p : ℝ[X]) (n : ℕ) :
    integralInner μ (IntegrablePair.of μ p (P n)) = (P.basis hu).repr p n * integralInner μ (IntegrablePair.of μ (P n) (P n)) := by
  classical
  let b := P.basis hu
  have hb (i : ℕ) : b i = P i := P.basis_eq_self hu i
  have hsum : (∑ i ∈ (b.repr p).support, b.repr p i • b i) = p := by
    simpa only [Finsupp.linearCombination_apply, Finsupp.sum] using b.linearCombination_repr p
  have h := congrArg (fun q => integralInner μ (IntegrablePair.of μ q (P n))) hsum
  rw [integralInner_sum_left hμ] at h
  have hsmul (p q : ℝ[X]) (a : ℝ) :=
    integralInner_smul_left μ p q a (h := hμ.integrable_mul p q)
  simp only [hsmul, hb] at h
  have hs : (∑ i ∈ (b.repr p).support,
      b.repr p i * integralInner μ (IntegrablePair.of μ (P i) (P n))) =
      b.repr p n * integralInner μ (IntegrablePair.of μ (P n) (P n)) := by
    apply Finset.sum_eq_single n
    · intro i _ hi
      rw [show integralInner μ (IntegrablePair.of μ (P i) (P n)) = 0 from hP.orthogonal i n hi, mul_zero]
    · intro hn
      rw [Finsupp.notMem_support_iff.mp hn, zero_mul]
  exact h.symm.trans hs

/-- Forward coefficients never vanish (DLMF 18.2(iv)). -/
theorem recurrenceA_ne_zero (P : Sequence ℝ) (n : ℕ) : P.recurrenceA n ≠ 0 :=
  div_ne_zero (leadingCoeff_ne_zero.mpr (P.ne_zero n))
    (leadingCoeff_ne_zero.mpr (P.ne_zero (n + 1)))

/-- Christoffel–Darboux with the factor `x-y` cleared (DLMF 18.2.12). -/
theorem IsOrthogonal.christoffelDarboux {μ : Measure ℝ} (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ) {P : Sequence ℝ} (hP : P.IsOrthogonal μ)
    (n : ℕ) (x y : ℝ) :
    (x - y) * P.kernel μ n x y (Polynomial.Sequence.KernelInput.of _ _ _) = P.recurrenceA n / P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _) *
      ((P (n + 1)).eval x * (P n).eval y - (P n).eval x * (P (n + 1)).eval y) := by
  induction n with
  | zero =>
    have hx := congrArg (fun p : ℝ[X] => p.eval x) (hP.two_term_zero hμ hpos)
    have hy := congrArg (fun p : ℝ[X] => p.eval y) (hP.two_term_zero hμ hpos)
    simp only [eval_mul, eval_X, eval_add, eval_smul, smul_eq_mul] at hx hy
    have h : (x - y) * (P 0).eval x * (P 0).eval y = P.recurrenceA 0 *
        ((P 1).eval x * (P 0).eval y - (P 0).eval x * (P 1).eval y) := by
      linear_combination (P 0).eval y * hx - (P 0).eval x * hy
    simp only [kernel_of_positive hpos P, Nat.zero_add, Finset.sum_range_one]
    calc
      _ = ((x - y) * (P 0).eval x * (P 0).eval y) / P.squaredNorm μ 0 (Polynomial.Sequence.IntegrableSquare.of _ _) := by ring
      _ = _ := by rw [h]; ring
  | succ n ih =>
    have hx := congrArg (fun p : ℝ[X] => p.eval x) (hP.three_term hμ hpos n)
    have hy := congrArg (fun p : ℝ[X] => p.eval y) (hP.three_term hμ hpos n)
    simp only [eval_mul, eval_X, eval_add, eval_smul, smul_eq_mul] at hx hy
    have h : (x - y) * (P (n + 1)).eval x * (P (n + 1)).eval y =
        P.recurrenceA (n + 1) * ((P (n + 2)).eval x * (P (n + 1)).eval y -
          (P (n + 1)).eval x * (P (n + 2)).eval y) -
        P.recurrenceC μ n (Sequence.RecurrenceCInput.of _ _ _) * ((P (n + 1)).eval x * (P n).eval y -
          (P n).eval x * (P (n + 1)).eval y) := by
      linear_combination (P (n + 1)).eval y * hx - (P (n + 1)).eval x * hy
    have hk : P.kernel μ (n + 1) x y (Polynomial.Sequence.KernelInput.of _ _ _) = P.kernel μ n x y (Polynomial.Sequence.KernelInput.of _ _ _) +
        (P (n + 1)).eval x * (P (n + 1)).eval y / P.squaredNorm μ (n + 1) (Polynomial.Sequence.IntegrableSquare.of _ _) := by
      simp only [kernel_of_positive hpos P]
      exact Finset.sum_range_succ _ _
    rw [hk, mul_add, ih]
    rw [show (x - y) * ((P (n + 1)).eval x * (P (n + 1)).eval y /
      P.squaredNorm μ (n + 1) (Polynomial.Sequence.IntegrableSquare.of _ _)) =
      ((x - y) * (P (n + 1)).eval x * (P (n + 1)).eval y) /
      P.squaredNorm μ (n + 1) (Polynomial.Sequence.IntegrableSquare.of _ _) by ring, h]
    simp only [recurrenceC]
    field_simp [ne_of_gt (squaredNorm_pos hpos P n),
      ne_of_gt (squaredNorm_pos hpos P (n + 1))]
    ring

/-- The diagonal Christoffel–Darboux identity (DLMF 18.2.13). -/
theorem IsOrthogonal.christoffelDarboux_diag {μ : Measure ℝ} (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ) {P : Sequence ℝ} (hP : P.IsOrthogonal μ)
    (n : ℕ) (x : ℝ) :
    P.kernel μ n x x (Polynomial.Sequence.KernelInput.of _ _ _) = P.recurrenceA n / P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _) *
      ((P (n + 1)).derivative.eval x * (P n).eval x -
        (P n).derivative.eval x * (P (n + 1)).eval x) := by
  have hk : DifferentiableAt ℝ (fun y => P.kernel μ n y x (Polynomial.Sequence.KernelInput.of _ _ _)) x := by
    unfold kernel
    exact DifferentiableAt.fun_sum fun k _ =>
      (((P k).differentiableAt).mul_const ((P k).eval x)).div_const _
  have hid : HasDerivAt (fun y : ℝ => y - x) 1 x := (hasDerivAt_id x).sub_const x
  have hl := hid.mul hk.hasDerivAt
  have hr := ((((P (n + 1)).hasDerivAt x).mul_const ((P n).eval x)).sub
    (((P n).hasDerivAt x).mul_const ((P (n + 1)).eval x))).const_mul
      (P.recurrenceA n / P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _))
  have heq : (fun y => (y - x) * P.kernel μ n y x (Polynomial.Sequence.KernelInput.of _ _ _)) =
      (fun y => P.recurrenceA n / P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _) *
        ((P (n + 1)).eval y * (P n).eval x -
          (P n).eval y * (P (n + 1)).eval x)) :=
    _root_.funext fun y => hP.christoffelDarboux hμ hpos n y x
  change HasDerivAt (fun y => (y - x) * P.kernel μ n y x (Polynomial.Sequence.KernelInput.of _ _ _)) _ x at hl
  rw [heq] at hl
  simpa using hl.unique hr

/-- Every irreducible monic factor of an orthogonal polynomial is linear
(Szegő, Chapter 3). -/
theorem IsOrthogonal.natDegree_eq_one_of_mem_normalizedFactors {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
    {P : Sequence ℝ} (hP : P.IsOrthogonal μ) (n : ℕ) {f : ℝ[X]}
    (hf : f ∈ normalizedFactors (P n)) : f.natDegree = 1 := by
  classical
  obtain ⟨hi, hm, q, hpq⟩ := (Polynomial.mem_normalizedFactors_iff (P.ne_zero n)).mp hf
  by_contra hdeg
  have hfpos (x : ℝ) : 0 < f.eval x :=
    zero_lt_eval_of_roots_lt_of_leadingCoeff_nonneg
      (fun _ hr => (hi.not_isRoot_of_natDegree_ne_one hdeg hr).elim)
      (by rw [hm.leadingCoeff]; norm_num)
  have hq : q ≠ 0 := by
    intro h
    exact P.ne_zero n (by rw [hpq, h, mul_zero])
  have hfd : 0 < f.natDegree := by
    exact natDegree_pos_iff_degree_pos.mpr (degree_pos_of_irreducible hi)
  have hqd : q.degree < (n : WithBot ℕ) := by
    rw [degree_eq_natDegree hq, ← P.natDegree_eq n, hpq,
      natDegree_mul hi.ne_zero hq]
    exact_mod_cast (show q.natDegree < f.natDegree + q.natDegree by omega)
  have hzero := hP.integralInner_eq_zero_of_degree_lt hμ hqd
  have hsupp : Function.support (fun x => q.eval x * (P n).eval x) =
      Function.support (fun x => q.eval x ^ 2) := by
    ext x
    simp only [Function.mem_support, hpq, eval_mul, ne_eq, mul_eq_zero]
    simp [ne_of_gt (hfpos x)]
  have hpositive : 0 < integralInner μ (IntegrablePair.of μ q (P n)) := by
    apply (integral_pos_iff_support_of_nonneg ?_ (hμ.integrable_mul q (P n))).mpr
    · rw [hsupp]
      exact (integral_pos_iff_support_of_nonneg (fun x => sq_nonneg (q.eval x))
        (hμ q)).mp (hpos q hq)
    · intro x
      change 0 ≤ q.eval x * (P n).eval x
      rw [hpq, eval_mul]
      nlinarith [mul_nonneg (le_of_lt (hfpos x)) (sq_nonneg (q.eval x))]
  exact (ne_of_gt hpositive) hzero

/-- A nowhere-zero Wronskian forces a zero of the denominator between two zeros
of the numerator. This is the Rolle argument behind the separation theorem
(Szegő, Chapter 3). -/
theorem exists_eval_eq_zero_between_of_wronskian_ne_zero {p q : ℝ[X]}
    (hw : ∀ x : ℝ, p.derivative.eval x * q.eval x - p.eval x * q.derivative.eval x ≠ 0)
    {a b : ℝ} (hab : a < b) (ha : p.eval a = 0) (hb : p.eval b = 0) :
    ∃ c ∈ Ioo a b, q.eval c = 0 := by
  by_contra hno
  push Not at hno
  have hq : ∀ x ∈ Icc a b, q.eval x ≠ 0 := by
    intro x hx
    by_cases hxa : x = a
    · subst x
      intro h
      exact hw a (by rw [ha, h]; ring)
    by_cases hxb : x = b
    · subst x
      intro h
      exact hw b (by rw [hb, h]; ring)
    exact hno x ⟨lt_of_le_of_ne hx.1 (Ne.symm hxa), lt_of_le_of_ne hx.2 hxb⟩
  obtain ⟨c, hc, hd⟩ := exists_hasDerivAt_eq_zero hab
    (p.differentiable.continuous.continuousOn.div
      q.differentiable.continuous.continuousOn hq)
    (show p.eval a / q.eval a = p.eval b / q.eval b by rw [ha, hb]; simp)
    (fun x hx => (p.hasDerivAt x).div (q.hasDerivAt x) (hq x (Ioo_subset_Icc_self hx)))
  exact hw c ((div_eq_zero_iff).mp hd |>.resolve_right
    (pow_ne_zero 2 (hq c (Ioo_subset_Icc_self hc))))

/-- The degree-`n+1` polynomial has a root strictly below every root of the degree-`n`
polynomial, and a distinct root strictly above every such root (DLMF 18.2(vi)). -/
theorem IsOrthogonal.exists_outer_roots {μ : Measure ℝ}
    (hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
    {P : Sequence ℝ} (hP : P.IsOrthogonal μ) {n : ℕ} (hn : 1 ≤ n) :
    ∃ a b : ℝ, a < b ∧ (P (n + 1)).eval a = 0 ∧ (P (n + 1)).eval b = 0 ∧
      ∀ x : ℝ, (P n).eval x = 0 → a < x ∧ x < b := by
  classical
  let s := (P n).roots.toFinset
  let t := (P (n + 1)).roots.toFinset
  have hs (x : ℝ) : x ∈ s ↔ (P n).eval x = 0 := by
    simp only [s, Multiset.mem_toFinset, mem_roots (P.ne_zero n), IsRoot]
  have ht (x : ℝ) : x ∈ t ↔ (P (n + 1)).eval x = 0 := by
    simp only [t, Multiset.mem_toFinset, mem_roots (P.ne_zero (n + 1)), IsRoot]
  have hc : s.card < t.card := by
    dsimp [s, t]
    rw [Multiset.toFinset_card_of_nodup (hP.roots_nodup hμ hpos n),
      Multiset.toFinset_card_of_nodup (hP.roots_nodup hμ hpos (n + 1)),
      hP.card_roots hμ hpos n, hP.card_roots hμ hpos (n + 1)]
    omega
  have hsep : ∀ a ∈ t, ∀ b ∈ t, a < b → ∃ c ∈ s, a < c ∧ c < b := by
    intro a ha b hb hab
    obtain ⟨c, hc, hz⟩ := hP.exists_root_between_succ hμ hpos n hab
      ((ht a).mp ha) ((ht b).mp hb)
    exact ⟨c, (hs c).mpr hz, hc⟩
  obtain ⟨a, ha, hal⟩ := exists_mem_lt_all_of_separates s t hc hsep
  obtain ⟨b, hb, hbr⟩ := exists_mem_lt_all_of_separates
    (α := OrderDual ℝ) s t hc (by
      intro u hu v hv hvu
      obtain ⟨c, hc, hvc, hcu⟩ := hsep v hv u hu hvu
      exact ⟨c, hc, hcu, hvc⟩)
  obtain ⟨x, hx⟩ := Multiset.card_pos_iff_exists_mem.mp (show 0 < (P n).roots.card by
    rw [hP.card_roots hμ hpos n]
    omega)
  have hxs : x ∈ s := Multiset.mem_toFinset.mpr hx
  exact ⟨a, b, (hal x hxs).trans (hbr x hxs), (ht a).mp ha, (ht b).mp hb,
    fun x hx ↦ ⟨hal x ((hs x).mpr hx), hbr x ((hs x).mpr hx)⟩⟩

/-- Positive-leading orthonormal sequences exist for every measure with finite
moments and a strictly positive polynomial pairing. -/
theorem exists_positive_orthonormal_sequence (μ : Measure ℝ) (hμ : HasFiniteMoments μ)
    (hpos : HasPositiveDefiniteIntegral μ) :
    ∃ P : Sequence ℝ, P.IsOrthogonal μ ∧
      (∀ n, 0 < (P n).leadingCoeff) ∧ (∀ n, ∫ x, (P n).eval x ^ 2 ∂μ = 1) := by
  obtain ⟨R, hR⟩ := exists_isOrthogonal μ hμ hpos
  have hr (n : ℕ) : (R n).leadingCoeff ≠ 0 := leadingCoeff_ne_zero.mpr (R.ne_zero n)
  let Q := R.rescale (fun n => (R n).leadingCoeff⁻¹) (fun n => inv_ne_zero (hr n))
  have hQ : Q.IsOrthogonal μ := hR.rescale _ _
  have hlead (n : ℕ) : (Q n).leadingCoeff = 1 := by
    change ((R n).leadingCoeff⁻¹ • R n).leadingCoeff = 1
    rw [leadingCoeff_smul_of_smul_regular _
      (IsSMulRegular.of_ne_zero (inv_ne_zero (hr n))), smul_eq_mul, inv_mul_cancel₀ (hr n)]
  let a (n : ℕ) := Real.sqrt (∫ x, (Q n).eval x ^ 2 ∂μ)
  have ha (n : ℕ) : 0 < a n := Real.sqrt_pos.mpr (hpos (Q n) (Q.ne_zero n))
  let P := Q.rescale (fun n => (a n)⁻¹) (fun n => inv_ne_zero (ha n).ne')
  refine ⟨P, hQ.rescale _ _, ?_, ?_⟩
  · intro n
    change 0 < ((a n)⁻¹ • Q n).leadingCoeff
    rw [leadingCoeff_smul_of_smul_regular _
      (IsSMulRegular.of_ne_zero (inv_ne_zero (ha n).ne')), smul_eq_mul, hlead, mul_one]
    exact inv_pos.mpr (ha n)
  · intro n
    change ∫ x, ((a n)⁻¹ • Q n).eval x ^ 2 ∂μ = 1
    simp only [eval_smul, smul_eq_mul, mul_pow, integral_const_mul]
    have hnorm : a n ^ 2 = ∫ x, (Q n).eval x ^ 2 ∂μ :=
      Real.sq_sqrt (hpos (Q n) (Q.ne_zero n)).le
    rw [← hnorm, ← mul_pow, inv_mul_cancel₀ (ha n).ne', one_pow]

/-- Weierstrass approximation gives dense polynomial classes in weighted L² (ATAP, Chapter 1). -/
theorem polynomialToL2_denseRange (a b : ℝ)
    (hsupp : ∀ᵐ x ∂μ, x ∈ Set.Icc a b) : DenseRange (polynomialToL2 hμ) := by
  let _ := hμ.isFiniteMeasure
  have hd := BoundedContinuousFunction.toLp_denseRange ℝ μ ℝ
    (by norm_num : (2 : ℝ≥0∞) ≠ ∞)
  apply Dense.of_closure (hd.mono ?_)
  rintro _ ⟨f, rfl⟩
  rw [Metric.mem_closure_iff]
  intro ε hε
  let K : ℝ := (measureUnivNNReal μ : ℝ) ^ (2 : ℝ≥0∞).toReal⁻¹
  have hK : 0 ≤ K := by positivity
  let δ := ε / (K + 1)
  have hδ : 0 < δ := div_pos hε (by positivity)
  obtain ⟨p, hp⟩ := exists_polynomial_near_of_continuousOn a b f
    f.continuous.continuousOn δ hδ
  refine ⟨polynomialToL2 hμ p, ⟨p, rfl⟩, ?_⟩
  rw [dist_eq_norm, norm_sub_rev]
  have hb : ‖polynomialToL2 hμ p - BoundedContinuousFunction.toLp 2 μ ℝ f‖ ≤ K * δ := by
    apply Lp.norm_le_of_ae_bound hδ.le
    filter_upwards [hsupp, coe_polynomialToL2 hμ p,
      f.coeFn_toLp 2 μ ℝ,
      Lp.coeFn_sub (polynomialToL2 hμ p) (BoundedContinuousFunction.toLp 2 μ ℝ f)]
      with x hx hpx hfx hsub
    rw [hsub]
    simp only [Pi.sub_apply, hpx, hfx, Real.norm_eq_abs]
    exact (hp x hx).le
  apply hb.trans_lt
  dsimp [δ]
  rw [← mul_div_assoc]
  apply (div_lt_iff₀ (by positivity : 0 < K + 1)).mpr
  nlinarith

/-- Normalising an orthogonal sequence gives an orthonormal family (DLMF 18.2(viii)). -/
theorem IsOrthogonal.orthonormal_normalizedToL2 (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) : Orthonormal ℝ (fun n ↦ P.normalizedToL2 hμ ⟨n, squaredNorm_pos hpos P n⟩) := by
  rw [orthonormal_iff_ite]
  intro n m
  simp only [normalizedToL2, inner_smul_left, inner_smul_right, conj_trivial,
    inner_polynomialToL2]
  by_cases h : n = m
  · subst m
    simp only [↓reduceIte]
    have hp := squaredNorm_pos hpos P n
    have hn : (∫ x, (P n).eval x * (P n).eval x ∂μ) = P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _) := rfl
    rw [hn]
    calc
      _ = ((√(P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _)))⁻¹ * (√(P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _)))⁻¹) * P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _) := by ring
      _ = 1 := by rw [← mul_inv, Real.mul_self_sqrt hp.le, inv_mul_cancel₀ hp.ne']
  · rw [hP.orthogonal n m h]
    simp [h]

/-- Hilbert-coordinate terms equal the integral-coefficient terms (DLMF 18.2(viii)). -/
theorem IsOrthogonal.hilbertBasis_repr_smul (hP : P.IsOrthogonal μ)
    (hpos : HasPositiveDefiniteIntegral μ) (hdense : DenseRange (polynomialToL2 hμ))
    (f : Lp ℝ 2 μ) (n : ℕ) :
    (hP.hilbertBasis hμ hpos hdense).repr f n • hP.hilbertBasis hμ hpos hdense n =
      P.integralCoefficient μ (IntegralCoefficientInput.ofL2 hμ hpos (P n) (P.ne_zero n) f) • polynomialToL2 hμ (P n) := by
  rw [HilbertBasis.repr_apply_apply, IsOrthogonal.hilbertBasis_apply, normalizedToL2,
    inner_smul_left, inner_polynomialToL2_left]
  simp only [conj_trivial, smul_smul]
  congr 1
  calc
    _ = ((√(P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _)))⁻¹ * (√(P.squaredNorm μ n (Polynomial.Sequence.IntegrableSquare.of _ _)))⁻¹) *
        (∫ x, f x * (P n).eval x ∂μ) := by ring
    _ = _ := by rw [← mul_inv, Real.mul_self_sqrt (squaredNorm_pos hpos P n).le]; rfl