/-- Exactness of Clenshaw evaluation for the first `n` coefficients. -/
theorem clenshaw_eq_sum (a b c p : ℕ → R)
    (hp : ∀ k, p (k + 2) = a k * p (k + 1) + b k * p k) (n : ℕ) :
    clenshaw a b c (p 0) (p 1) n = ∑ i ∈ Finset.range n, c i * p i := by
  simpa [clenshaw] using clenshawState_spec a b c p hp n 0

/-- An empty expansion evaluates to zero. -/
@[simp] theorem clenshaw_zero (a b c : ℕ → R) (p₀ p₁ : R) :
    clenshaw a b c p₀ p₁ 0 = 0 := by
  simp [clenshaw, clenshawState]

/-- A constant expansion uses the full initial coefficient. -/
@[simp] theorem clenshaw_one (a b c : ℕ → R) (p₀ p₁ : R) :
    clenshaw a b c p₀ p₁ 1 = c 0 * p₀ := by
  simp [clenshaw, clenshawState]

/-- Executable Clenshaw evaluation of the first `n` Chebyshev coefficients. -/
def clenshawT (c : ℕ → R) (x : R) (n : ℕ) : R :=
  clenshaw (fun _ ↦ 2 * x) (fun _ ↦ -1) c 1 x n

/-- Finite coefficient integration is exact (Chebfun cumsum's polynomial core). -/
theorem derivative_sum_primitiveT (s : Finset ℕ) (a : ℕ → K) :
    derivative (∑ n ∈ s, Polynomial.C (a n) * primitiveT K n) = ∑ n ∈ s, Polynomial.C (a n) * T K n := by
  simp only [derivative_sum, derivative_C_mul, derivative_primitiveT]

/-- Anchoring does not alter the derivative. -/
@[simp] theorem derivative_primitiveTAt (a : K) (n : ℕ) :
    derivative (primitiveTAt a n) = T K n := by
  simp [primitiveTAt, derivative_primitiveT]

/-- The anchored primitive vanishes at the chosen starting point. -/
@[simp] theorem eval_primitiveTAt_self (a : K) (n : ℕ) (d : PrimitiveInput K n) :
    (primitiveTAt a n d).eval a = 0 := by simp [primitiveTAt]

/-- The negative-endpoint Neumann row in the full Chebyshev basis
(Olver–Townsend, boundary bordering). -/
theorem derivative_T_succ_eval_neg_one (n : ℕ) :
    (derivative (T K (n + 1 : ℕ))).eval (-1) = (-1 : K) ^ n * (n + 1 : K) ^ 2 := by
  rw [T_derivative_eq_U]
  simp only [eval_mul, eval_intCast, Nat.cast_add, Nat.cast_one,
    Int.add_sub_cancel, U_eval_neg_one, Int.cast_natCast]
  simp only [Int.cast_negOnePow_natCast, Int.cast_add, Int.cast_natCast, Int.cast_one]
  ring

/-- The first derivative enters the C^(1) basis, with the constant mode removed. -/
theorem hasSum_gegenbauer_derivative_complexChebyshevSeries {a : ℕ → ℂ} {C ρ : ℝ}
    (hρ : 1 < ρ) (ha : ∀ n, ‖a n‖ ≤ C / ρ ^ n) (z : bernsteinEllipse ⟨ρ, hρ⟩) :
    HasSum (fun n ↦ a (n + 1) * (n + 1 : ℂ) * (gegenbauer (1 : ℂ) n).eval z.val)
      (iteratedDerivOnOpen (complexChebyshevSeriesDerivInput hρ ha 1) z) := by
  have h := (hasSum_nat_add_iff' 1).2 (hasSum_derivative_complexChebyshevSeries hρ ha z)
  simpa only [Finset.sum_range_one, Nat.cast_zero, T_zero, derivative_one, eval_zero,
    mul_zero, sub_zero, Nat.cast_add, Nat.cast_one, derivative_chebyshev_T_succ,
    eval_mul, eval_C, mul_assoc] using h

/-- Uniform convergence on compact subsets, with the limit evaluated only on its domain. -/
theorem tendstoUniformlyOn_derivative_complexChebyshevSeries {a : ℕ → ℂ} {C ρ r : ℝ}
    (hρ : 1 < ρ) (hr : 1 < r) (hrρ : r < ρ) (ha : ∀ n, ‖a n‖ ≤ C / ρ ^ n)
    {K : Set ℂ} (hK : IsCompact K) (hKr : K ⊆ bernsteinEllipse ⟨r, hr⟩) :
    TendstoUniformlyOn
      (fun N (z : bernsteinEllipse ⟨ρ, hρ⟩) ↦
        ∑ n ∈ Finset.range N, a n * (T ℂ n).derivative.eval z.val)
      (iteratedDerivOnOpen (complexChebyshevSeriesDerivInput hρ ha 1))
      atTop {z | z.val ∈ K} := by
  have h := (tendstoUniformlyOn_derivative_series (isOpen_bernsteinEllipse ⟨r, hr⟩)
    ((summable_geometric_of_lt_one (by positivity : 0 ≤ r / ρ)
      ((div_lt_one (by linarith : 0 < ρ)).2 hrρ)).mul_left C)
    (norm_chebyshev_term_le_geometric hρ hr ha) hK hKr).comp
    (fun z : bernsteinEllipse ⟨ρ, hρ⟩ ↦ z.val)
  have he (N : ℕ) (z : ℂ) :
      deriv (fun w ↦ ∑ n ∈ Finset.range N, a n * (T ℂ n).eval w) z =
        ∑ n ∈ Finset.range N, a n * (T ℂ n).derivative.eval z := by
    rw [deriv_fun_sum (fun (n : ℕ) _ ↦ ((T ℂ n).differentiable.const_mul (a n)).differentiableAt)]
    simp only [deriv_const_mul_field, Polynomial.deriv]
  simp only [Function.comp_def, he] at h
  exact h.congr_right (fun z _ ↦ by
    rw [iteratedDeriv_complexChebyshevSeries_eq, iteratedDeriv_one])

/-- The corresponding uniform first-order bound. -/
theorem exists_uniform_derivative_series_bound {a : ℕ → ℂ} {C ρ r : ℝ}
    (hρ : 1 < ρ) (hr : 1 < r) (hrρ : r < ρ) (ha : ∀ n, ‖a n‖ ≤ C / ρ ^ n)
    {K : Set ℂ} (hK : IsCompact K) (hKr : K ⊆ bernsteinEllipse ⟨r, hr⟩) :
    ∃ δ : ℝ, 0 < δ ∧ ∀ z : bernsteinEllipse ⟨ρ, hρ⟩, z.val ∈ K → ∀ N : ℕ,
      ‖iteratedDerivOnOpen (complexChebyshevSeriesDerivInput hρ ha 1) z -
        ∑ n ∈ Finset.range N, a n * (T ℂ n).derivative.eval z.val‖ ≤
        C * (r / ρ) ^ N / (1 - r / ρ) / δ := by
  obtain ⟨δ, hδ, hsub⟩ := hK.exists_cthickening_subset_open (isOpen_bernsteinEllipse ⟨r, hr⟩) hKr
  refine ⟨δ, hδ, fun z hz N ↦ ?_⟩
  apply norm_derivative_complexChebyshevSeries_sub_sum_le hρ hr hrρ ha z (hKr hz) hδ
  exact Metric.sphere_subset_closedBall.trans ((Metric.closedBall_subset_cthickening hz δ).trans hsub)

/-- Anchored modes are uniformly bounded by the interval length. -/
theorem norm_eval_primitiveTAt_le_two (n : ℕ) {x : ℝ} (hx : x ∈ Icc (-1 : ℝ) 1) :
    ‖(primitiveTAt (-1 : ℂ) n).eval (x : ℂ)‖ ≤ 2 := by
  rw [eval_primitiveTAt_eq_integral]
  have h := intervalIntegral.norm_integral_le_of_norm_le_const
    (a := (-1 : ℝ)) (b := x) (C := 1) (f := fun t : ℝ ↦ (T ℂ n).eval (t : ℂ))
    (fun t ht ↦ norm_eval_T_complex_le_one n (by
      rw [uIoc_of_le hx.1] at ht
      exact ⟨ht.1.le, ht.2.trans hx.2⟩))
  have he : |x - (-1)| ≤ 2 := by rw [abs_of_nonneg (by linarith [hx.1])]; linarith [hx.2]
  exact h.trans (by simpa using he)

/-- Termwise anchored integration, with actual integrability included. -/
theorem hasSum_primitiveTAt {a : ℕ → ℂ} (ha : Summable (fun n ↦ ‖a n‖))
    {f : ℝ → ℂ}
    (hf : ∀ x : Icc (-1 : ℝ) 1, HasSum (fun n ↦ a n * (T ℂ n).eval (x.val : ℂ)) (f x.val))
    (x : Icc (-1 : ℝ) 1) :
    IntervalIntegrable f volume (-1) x.val ∧
    HasSum (fun n ↦ a n * (primitiveTAt (-1 : ℂ) n).eval (x.val : ℂ))
      (∫ t in (-1 : ℝ)..x.val, f t) := by
  have hsub : uIcc (-1 : ℝ) x.val ⊆ Icc (-1 : ℝ) 1 := by
    rw [uIcc_of_le x.property.1]
    exact Icc_subset_Icc le_rfl x.property.2
  refine ⟨((continuousOn_of_hasSum_chebyshev ha hf).mono hsub).intervalIntegrable, ?_⟩
  have hb (n : ℕ) (t : ℝ) (ht : t ∈ uIcc (-1) x.val) :
      ‖a n * (T ℂ n).eval (t : ℂ)‖ ≤ ‖a n‖ := by
    rw [norm_mul]
    exact mul_le_of_le_one_right (norm_nonneg _) (norm_eval_T_complex_le_one n (hsub ht))
  have h := hasSum_integral_series (F := fun n t ↦ a n * (T ℂ n).eval (t : ℂ))
    (fun n ↦ by fun_prop) ha hb (fun t ht ↦ hf ⟨t, hsub ht⟩)
  simpa only [intervalIntegral.integral_const_mul, ← eval_primitiveTAt_eq_integral] using h

/-- Uniform convergence includes both endpoints of the typed interval. -/
theorem tendstoUniformlyOn_primitiveTAt {a : ℕ → ℂ} (ha : Summable (fun n ↦ ‖a n‖))
    {f : ℝ → ℂ}
    (hf : ∀ x : Icc (-1 : ℝ) 1, HasSum (fun n ↦ a n * (T ℂ n).eval (x.val : ℂ)) (f x.val)) :
    TendstoUniformly
      (fun N (x : Icc (-1 : ℝ) 1) ↦ ∑ n ∈ Finset.range N, a n * (primitiveTAt (-1 : ℂ) n).eval (x.val : ℂ))
      (fun x ↦ ∫ t in (-1 : ℝ)..x.val, f t) atTop := by
  have h := tendstoUniformlyOn_tsum_nat (s := Set.univ) (ha.mul_left 2)
    (fun n (x : Icc (-1 : ℝ) 1) _ ↦
      (show ‖a n * (primitiveTAt (-1 : ℂ) n).eval (x.val : ℂ)‖ ≤ 2 * ‖a n‖ by
        rw [norm_mul]; simpa only [mul_comm] using
          mul_le_mul_of_nonneg_left (norm_eval_primitiveTAt_le_two n x.property) (norm_nonneg (a n))))
  exact tendstoUniformlyOn_univ.mp (h.congr_right (fun x _ ↦ (hasSum_primitiveTAt ha hf x).2.tsum_eq))

/-- Geometric input coefficients give a geometric anchored-integration error,
uniformly including both endpoints. -/
theorem norm_integral_sub_sum_primitiveTAt_le_geometric {a : ℕ → ℂ} {C ρ : ℝ}
    (hρ : 1 < ρ) (ha : ∀ n, ‖a n‖ ≤ C / ρ ^ n)
    {f : ℝ → ℂ}
    (hf : ∀ x : Icc (-1 : ℝ) 1, HasSum (fun n ↦ a n * (T ℂ n).eval (x.val : ℂ)) (f x.val))
    (x : Icc (-1 : ℝ) 1) (N : ℕ) :
    ‖(∫ t in (-1 : ℝ)..x.val, f t) -
      ∑ n ∈ Finset.range N, a n * (primitiveTAt (-1 : ℂ) n).eval (x.val : ℂ)‖ ≤
      2 * (C * (ρ⁻¹) ^ N / (1 - ρ⁻¹)) := by
  have hq0 : 0 ≤ ρ⁻¹ := by positivity
  have hq1 : ρ⁻¹ < 1 := (inv_lt_one₀ (by linarith : 0 < ρ)).2 hρ
  have hb (n : ℕ) : ‖a n‖ ≤ C * (ρ⁻¹) ^ n := by
    simpa only [inv_pow, div_eq_mul_inv] using ha n
  have hs := Summable.of_nonneg_of_le (fun _ ↦ norm_nonneg _) hb
    ((summable_geometric_of_lt_one hq0 hq1).mul_left C)
  exact (norm_integral_sub_sum_primitiveTAt_le hs hf x N).trans
    (mul_le_mul_of_nonneg_left
      (tsum_tail_le_geometric hq0 hq1 (fun _ ↦ norm_nonneg _) hb N) (by norm_num))

/-- The absolute coefficient norm is submultiplicative. -/
theorem tsum_norm_chebyshevProductCoefficients_le {a b : ℕ → ℂ}
    (ha : Summable (fun n ↦ ‖a n‖)) (hb : Summable (fun n ↦ ‖b n‖)) :
    ∑' k, ‖chebyshevProductCoefficients (ChebyshevProductInput.ofSummableNorm ha hb k)‖ ≤
      (∑' n, ‖a n‖) * ∑' n, ‖b n‖ := by
  let c (p : ℕ × ℕ) := a p.1 * b p.2
  let e (p : ℕ × ℕ) := p.1 + p.2
  let d (p : ℕ × ℕ) := ((p.1 : ℤ) - p.2).natAbs
  have hc : Summable (fun p ↦ ‖c p‖) := ha.mul_norm hb
  have he := Polynomial.summable_norm_regroupCoefficients hc e
  have hd := Polynomial.summable_norm_regroupCoefficients hc d
  have ht := (summable_norm_chebyshevProductCoefficients ha hb).tsum_le_tsum
    (g := fun k ↦ (‖Polynomial.regroupCoefficients (RegroupCoefficientInput.ofSummableNorm hc e k)‖ + ‖Polynomial.regroupCoefficients (RegroupCoefficientInput.ofSummableNorm hc d k)‖) / 2)
    (fun k ↦ by
      rw [chebyshevProductCoefficients, norm_div, norm_ofNat]
      exact div_le_div_of_nonneg_right (norm_add_le _ _) (by norm_num))
    ((he.add hd).div_const 2)
  rw [tsum_div_const, he.tsum_add hd] at ht
  have h1 := Polynomial.tsum_norm_regroupCoefficients_le hc e
  have h2 := Polynomial.tsum_norm_regroupCoefficients_le hc d
  have hm : (∑' p, ‖c p‖) = (∑' n, ‖a n‖) * ∑' n, ‖b n‖ := by
    simpa only [c, norm_mul] using (tsum_mul_tsum_of_summable_norm ha.norm hb.norm).symm
  rw [hm] at h1 h2
  linarith

/-- The collected product coefficients represent the pointwise product on the
whole closed interval. All rearrangements are absolutely convergent. -/
theorem hasSum_chebyshevProductCoefficients {a b : ℕ → ℂ}
    (ha : Summable (fun n ↦ ‖a n‖)) (hb : Summable (fun n ↦ ‖b n‖))
    (x : Icc (-1 : ℝ) 1) :
    HasSum (fun k ↦ chebyshevProductCoefficients (ChebyshevProductInput.ofSummableNorm ha hb k) * (T ℂ k).eval (x : ℂ))
      (complexChebyshevSeries a (.ofAbsolute ha x) * complexChebyshevSeries b (.ofAbsolute hb x)) := by
  let c (p : ℕ × ℕ) := a p.1 * b p.2
  let e (p : ℕ × ℕ) := p.1 + p.2
  let d (p : ℕ × ℕ) := ((p.1 : ℤ) - p.2).natAbs
  have hc : Summable (fun p ↦ ‖c p‖) := ha.mul_norm hb
  have hs := Polynomial.hasSum_regrouped_series hc e (fun k ↦ norm_eval_T_complex_le_one k x.property)
  have hd := Polynomial.hasSum_regrouped_series hc d (fun k ↦ norm_eval_T_complex_le_one k x.property)
  have heval (v : ℕ → ℂ) (hv : Summable (fun n ↦ ‖v n‖)) :
      Summable (fun n ↦ ‖v n * (T ℂ n).eval (x : ℂ)‖) := by
    apply Summable.of_nonneg_of_le (fun _ ↦ norm_nonneg _) _ hv
    intro n
    rw [norm_mul]
    exact mul_le_of_le_one_right (norm_nonneg _) (norm_eval_T_complex_le_one n x.property)
  have htotal : (∑' p, c p * (T ℂ (e p)).eval (x : ℂ)) +
      (∑' p, c p * (T ℂ (d p)).eval (x : ℂ)) =
      2 * (complexChebyshevSeries a (.ofAbsolute ha x) * complexChebyshevSeries b (.ofAbsolute hb x)) := by
    have hsum (j : ℕ × ℕ → ℕ) : Summable (fun p ↦ c p * (T ℂ (j p)).eval (x : ℂ)) := by
      apply Summable.of_norm
      apply Summable.of_nonneg_of_le (fun _ ↦ norm_nonneg _) _ hc
      intro p
      rw [norm_mul]
      exact mul_le_of_le_one_right (norm_nonneg _) (norm_eval_T_complex_le_one _ x.property)
    rw [← (hsum e).tsum_add (hsum d)]
    rw [complexChebyshevSeries_eq_tsum (.ofAbsolute ha x) (summable_chebyshevSeries_of_absolute ha x),
      complexChebyshevSeries_eq_tsum (.ofAbsolute hb x) (summable_chebyshevSeries_of_absolute hb x)]
    simp only [ComplexChebyshevSeriesPoint.ofAbsolute_val]
    rw [tsum_mul_tsum_of_summable_norm (heval a ha) (heval b hb), ← tsum_mul_left]
    apply tsum_congr
    intro p
    dsimp [c, e, d]
    have hm := product_mode p.1 p.2 (x : ℂ)
    simp only [Nat.cast_add] at hm ⊢
    linear_combination -2 * a p.1 * b p.2 * hm
  have h := (hs.add hd).div_const 2
  rw [htotal] at h
  convert h using 1
  · ext k
    dsimp [chebyshevProductCoefficients, ChebyshevProductInput.ofSummableNorm,
      ChebyshevProductInput.addInput, Polynomial.regroupCoefficients,
      RegroupCoefficientInput.ofSummableNorm, RegroupCoefficientInput.ofSummable, c, e, d]
    ring
  · ring

/-- The product expansion is an equality of the represented functions. -/
theorem complexChebyshevSeries_product {a b : ℕ → ℂ}
    (ha : Summable (fun n ↦ ‖a n‖)) (hb : Summable (fun n ↦ ‖b n‖))
    (x : Icc (-1 : ℝ) 1) :
    complexChebyshevSeries ((fun k ↦ chebyshevProductCoefficients (ChebyshevProductInput.ofSummableNorm ha hb k)))
      (.ofAbsolute (summable_norm_chebyshevProductCoefficients ha hb) x) =
      complexChebyshevSeries a (.ofAbsolute ha x) * complexChebyshevSeries b (.ofAbsolute hb x) := by
  rw [complexChebyshevSeries_eq_tsum
    (.ofAbsolute (summable_norm_chebyshevProductCoefficients ha hb) x)
    (summable_chebyshevSeries_of_absolute (summable_norm_chebyshevProductCoefficients ha hb) x)]
  exact (hasSum_chebyshevProductCoefficients ha hb x).tsum_eq

/-- Integral coefficients of an interval representative of the product. -/
theorem complexCoefficient_product_series {a b : ℕ → ℂ}
    (ha : Summable (fun n ↦ ‖a n‖)) (hb : Summable (fun n ↦ ‖b n‖))
    (k : ℕ) (f : ComplexCoefficientInput k)
    (he : ∀ x : Icc (-1 : ℝ) 1, f.val x.val =
      complexChebyshevSeries a (.ofAbsolute ha x) * complexChebyshevSeries b (.ofAbsolute hb x)) :
    complexCoefficient f = chebyshevProductCoefficients (ChebyshevProductInput.ofSummableNorm ha hb k) := by
  apply complexCoefficient_complexChebyshevSeries (summable_norm_chebyshevProductCoefficients ha hb) k f
  intro x
  exact (he x).trans (complexChebyshevSeries_product ha hb x).symm

/-- Absolutely summable Chebyshev factors have a vanishing product-coefficient tail. -/
theorem tendsto_chebyshevProductCoefficients_tail {a b : ℕ → ℂ}
    (ha : Summable (fun n ↦ ‖a n‖)) (hb : Summable (fun n ↦ ‖b n‖)) :
    Tendsto (fun N ↦ ∑' k, ‖chebyshevProductCoefficients
      (ChebyshevProductInput.ofSummableNorm ha hb (k + N))‖) atTop (𝓝 0) :=
  tendsto_sum_nat_add (fun k ↦ ‖chebyshevProductCoefficients
    (ChebyshevProductInput.ofSummableNorm ha hb k)‖)

@[simp] theorem intervalMap_one (a b : ℝ) : intervalMap a b 1 = b := by
  unfold intervalMap
  ring

/-- Physical and canonical coordinates cancel for a nondegenerate interval. -/
@[simp] theorem intervalMap_intervalCoord {a b : ℝ} (hab : a < b) (x : ℝ) :
    intervalMap a b (intervalCoord a ⟨b, ne_of_lt hab⟩ x) = x := by
  unfold intervalMap intervalCoord
  field_simp [ne_of_gt (sub_pos.mpr hab)]
  ring

/-- Canonical and physical coordinates cancel for a nondegenerate interval. -/
@[simp] theorem intervalCoord_intervalMap {a b : ℝ} (hab : a < b) (t : ℝ) :
    intervalCoord a ⟨b, ne_of_lt hab⟩ (intervalMap a b t) = t := by
  unfold intervalMap intervalCoord
  field_simp [ne_of_gt (sub_pos.mpr hab)]
  ring

/-- The forward map preserves closed-interval membership. -/
theorem intervalMap_mem_Icc {a b t : ℝ} (hab : a < b) (ht : t ∈ Icc (-1 : ℝ) 1) :
    intervalMap a b t ∈ Icc a b := by
  unfold intervalMap
  constructor <;> nlinarith [mul_nonneg (sub_nonneg.mpr hab.le) (sub_nonneg.mpr ht.2),
    mul_nonneg (sub_nonneg.mpr hab.le) (show 0 ≤ t + 1 by linarith [ht.1])]

/-- The inverse map preserves closed-interval membership. -/
theorem intervalCoord_mem_Icc {a b x : ℝ} (hab : a < b) (hx : x ∈ Icc a b) :
    intervalCoord a ⟨b, ne_of_lt hab⟩ x ∈ Icc (-1 : ℝ) 1 := by
  unfold intervalCoord
  constructor
  · apply (le_div_iff₀ (sub_pos.mpr hab)).2
    linarith [hx.1]
  · apply (div_le_iff₀ (sub_pos.mpr hab)).2
    linarith [hx.2]

/-- Derivative of the physical coordinate map. -/
theorem hasDerivAt_intervalMap (a b t : ℝ) :
    HasDerivAt (intervalMap a b) ((b - a) / 2) t := by
  unfold intervalMap
  convert ((hasDerivAt_id t).const_mul ((b - a) / 2)).add_const ((a + b) / 2) using 1 <;> simp

/-- Derivative of the canonical coordinate, giving the factor `2/(b-a)`. -/
theorem hasDerivAt_intervalCoord (a : ℝ) (b : {b : ℝ // a ≠ b}) (x : ℝ) :
    HasDerivAt (intervalCoord a b) (2 / (b.val - a)) x := by
  unfold intervalCoord
  convert (((hasDerivAt_id x).const_mul 2).sub_const (a + b.val)).div_const (b.val - a) using 1 <;> simp

/-- Differentiation on a physical interval is scaled by `2/(b-a)`. -/
theorem hasDerivAt_comp_intervalCoord {a x d : ℝ} {b : {b : ℝ // a ≠ b}} {f : ℝ → ℝ}
    (hf : HasDerivAt f d (intervalCoord a b x)) :
    HasDerivAt (fun y ↦ f (intervalCoord a b y)) (d * (2 / (b.val - a))) x :=
  hf.comp x (hasDerivAt_intervalCoord a b x)

/-- Integrability transports through the affine map, including its constant case. -/
theorem intervalIntegrable_comp_intervalMap (a b : ℝ) (f : IntervalIntegrableInput a b) :
    IntervalIntegrable (fun t ↦ f.val (intervalMap a b t)) volume (-1) 1 := by
  by_cases hab : a = b
  · subst b
    simp [intervalMap]
  have h := (f.property.comp_add_right ((a + b) / 2)).comp_mul_left (c := (b - a) / 2)
  have ha : (a - (a + b) / 2) / ((b - a) / 2) = -1 := by
    field_simp [sub_ne_zero.mpr (Ne.symm hab)]; ring
  have hb : (b - (a + b) / 2) / ((b - a) / 2) = 1 := by
    field_simp [sub_ne_zero.mpr (Ne.symm hab)]; ring
  simpa only [ha, hb, intervalMap] using h

/-- The classical oriented change of variables for integrable input. -/
theorem integral_comp_intervalMap (a b : ℝ) (f : IntervalIntegrableInput a b) :
    (∫ x in a..b, f.val x) =
      (b - a) / 2 * ∫ t in (-1 : ℝ)..1, f.val (intervalMap a b t) := by
  have h := intervalIntegral.smul_integral_comp_mul_add
    (a := (-1 : ℝ)) (b := 1) f.val ((b - a) / 2) ((a + b) / 2)
  have h₀ : (b - a) / 2 * (-1) + (a + b) / 2 = a := by ring
  have h₁ : (b - a) / 2 * 1 + (a + b) / 2 = b := by ring
  simpa only [h₀, h₁, smul_eq_mul, intervalMap] using h.symm

/-- A uniform approximation bound transports without changing its constant. -/
theorem abs_sub_comp_intervalCoord_le {a b ε : ℝ} (hab : a < b)
    {f g : ℝ → ℝ}
    (he : ∀ t ∈ Icc (-1 : ℝ) 1, |f (intervalMap a b t) - g t| ≤ ε)
    {x : ℝ} (hx : x ∈ Icc a b) : |f x - g (intervalCoord a ⟨b, ne_of_lt hab⟩ x)| ≤ ε := by
  simpa only [intervalMap_intervalCoord hab] using he _ (intervalCoord_mem_Icc hab hx)

/-- Transported polynomial values are ordinary polynomial composition; no new
interval-parametrised polynomial family is introduced. -/
theorem eval_comp_intervalCoord (p : ℝ[X]) (a : ℝ) (b : {b : ℝ // a ≠ b}) (x : ℝ) :
    (p.comp (Polynomial.C (2 / (b.val - a)) * X - Polynomial.C ((a + b.val) / (b.val - a)))).eval x =
      p.eval (intervalCoord a b x) := by
  rw [eval_comp]
  congr 1
  simp only [eval_sub, eval_mul, eval_C, eval_X, intervalCoord]
  ring

/-- Nondegenerate affine transport preserves polynomial degree. -/
theorem degree_comp_intervalCoord (p : ℝ[X]) {a b : ℝ} (hab : a < b) :
    (p.comp (Polynomial.C (2 / (b - a)) * X - Polynomial.C ((a + b) / (b - a)))).degree =
      p.degree := by
  have hc : 2 / (b - a) ≠ 0 := div_ne_zero (by norm_num) (ne_of_gt (sub_pos.mpr hab))
  have hd : (Polynomial.C (2 / (b - a)) * X -
      Polynomial.C ((a + b) / (b - a))).degree = 1 := by
    rw [degree_sub_eq_left_of_degree_lt (degree_C_lt_degree_C_mul_X hc), degree_C_mul_X hc]
  rw [degree_comp (by rw [hd]; norm_num), hd, mul_one]