/-- Every admissible competitor bounds the best error from above. -/
theorem bestUniformError_le (N : ℕ) (f : C(Icc (-1 : ℝ) 1, ℂ)) {p : ℂ[X]}
    (hp : p.degree < ((N + 1 : ℕ) : WithBot ℕ)) :
    bestUniformError N f ≤ ‖f - polynomialOnInterval p‖ := by
  apply csInf_le
  · exact ⟨0, fun e ⟨q, _, he⟩ ↦ he ▸ norm_nonneg _⟩
  · exact ⟨p, hp, rfl⟩

/-- A polynomial changing sign in each of `n` consecutive intervals has degree
at least `n` (ATAP, proof of Theorem 10.1). -/
theorem le_natDegree_of_alternating_eval {n : ℕ} {p : ℝ[X]}
    {x : Fin (n + 1) → ℝ} (hx : StrictMono x)
    (hp : ∀ i : Fin n, p.eval (x i.castSucc) * p.eval (x i.succ) < 0) :
    n ≤ p.natDegree := by
  by_cases hn : n = 0
  · omega
  have hex (i : Fin n) := exists_eval_eq_zero_between_of_mul_neg
    (hx (Fin.castSucc_lt_succ (i := i))) (hp i)
  choose z hz he using hex
  have hzmono : StrictMono z := by
    intro i j hij
    have hmid : x i.succ ≤ x j.castSucc := hx.monotone (by simpa using hij)
    exact hz i |>.2.trans_le hmid |>.trans (hz j).1
  by_contra hdeg
  have hpzero := p.eq_zero_of_natDegree_lt_card_of_eval_eq_zero hzmono.injective he
    (by simpa using Nat.lt_of_not_ge hdeg)
  have ht := hp ⟨0, by omega⟩
  simp [hpzero] at ht

/-- In a finite-dimensional real normed space, the convex hull of a compact set
is compact. This is the compact form of Carathéodory's theorem. -/
theorem IsCompact.convexHull_finiteDimensional {E : Type*} [NormedAddCommGroup E]
    [NormedSpace ℝ E] [FiniteDimensional ℝ E] {s : Set E} (hs : IsCompact s) :
    IsCompact (convexHull ℝ s) := by
  have heq : convexHull ℝ s = ⋃ n ∈ Finset.range (Module.finrank ℝ E + 2), convexImage s n := by
    ext y
    constructor
    · intro hy
      obtain ⟨n, hn, w, z, hw, hsum, hz, heq⟩ := mem_convexHull_exists_fin_combination hy
      exact mem_iUnion.2 ⟨n, mem_iUnion.2 ⟨by simp only [Finset.mem_range]; omega,
        w, z, ⟨hw, hsum⟩, hz, heq⟩⟩
    · intro hy
      obtain ⟨n, _, w, z, hw, hz, rfl⟩ := mem_iUnion.1 hy |>.imp (fun _ h ↦ mem_iUnion.1 h)
      exact (convex_convexHull ℝ s).sum_mem (fun i _ ↦ hw.1 i) hw.2
        (fun i _ ↦ subset_convexHull ℝ s (hz i))
  rw [heq]
  apply Finset.isCompact_biUnion
  intro n _
  let S : Set (Fin n → ℝ) := {w | (∀ i, 0 ≤ w i) ∧ ∑ i, w i = 1}
  have hS : IsCompact S := by
    have h := isCompact_range (show Continuous
      (fun t : Convexity.StdSimplex ℝ (Fin n) ↦ (t.weights : Fin n → ℝ)) by fun_prop)
    rw [Convexity.StdSimplex.range_toFun_comp_weights] at h
    have heq : ((⋂ i : Fin n, {w : Fin n → ℝ | 0 ≤ w i}) ∩
        {w | ∑ i, w i = 1}) = S := by ext w; simp [S]
    exact heq ▸ h
  let K := S ×ˢ Set.pi univ (fun _ : Fin n ↦ s)
  have hK : IsCompact K := hS.prod (isCompact_univ_pi (fun _ ↦ hs))
  have him : convexImage s n = (fun wz : (Fin n → ℝ) × (Fin n → E) ↦
      ∑ i, wz.1 i • wz.2 i) '' K := by
    ext y
    simp [convexImage, K, S, Set.mem_pi, Set.mem_image, Prod.exists, and_assoc]
  rw [him]
  exact hK.image (by fun_prop)

/-- A nonzero best error has a positive norming certificate supported at at
most `n+2` distinct extreme points (ATAP 10.1). -/
theorem exists_extremeMoment_certificate {a b : ℝ} {n : ℕ} (e : C(Icc a b, ℝ))
    (he : 0 < ‖e‖)
    (hmin : ∀ q : ℝ[X], q.degree ≤ n → ‖e‖ ≤ ‖e - q.toContinuousMapOn (Icc a b)‖) :
    ∃ (m : ℕ) (x : Fin m → Icc a b) (w : Fin m → ℝ),
      m ≤ n + 2 ∧ Function.Injective x ∧ (∀ i, 0 < w i) ∧
      (∑ i, w i) = 1 ∧ (∀ i, |e (x i)| = ‖e‖) ∧
      ∑ i, w i • extremeMoment n e (x i) = 0 := by
  obtain ⟨ι, hι, z, w, hz, hi, hw, hsum, heq⟩ :=
    eq_pos_convex_span_of_mem_convexHull (zero_mem_convexHull_extremeMoment e he hmin)
  have hx (i : ι) : ∃ x : Icc a b, |e x| = ‖e‖ ∧ extremeMoment n e x = z i := hz ⟨i, rfl⟩
  choose x hxe hxz using hx
  have hinj : Function.Injective x := by
    intro i j hij
    apply hi.injective
    rw [← hxz i, ← hxz j, hij]
  let f := (Fintype.equivFin ι).symm
  refine ⟨Fintype.card ι, x ∘ f, w ∘ f, ?_, hinj.comp f.injective,
    fun i ↦ hw _, ?_, fun i ↦ hxe _, ?_⟩
  · have hcard := hi.card_le_finrank_succ.trans (Nat.add_le_add_right (Submodule.finrank_le _) 1)
    simpa using hcard
  · simpa only [Function.comp_apply] using (Equiv.sum_comp f w).trans hsum
  · simpa only [Function.comp_apply, hxz] using
      (Equiv.sum_comp f (fun i ↦ w i • z i)).trans heq

/-- An injective finite family in a linear order can be reindexed increasingly. -/
theorem exists_strictMono_reindex {α : Type*} [LinearOrder α] {m : ℕ}
    {x : Fin m → α} (hx : Function.Injective x) :
    ∃ e : Equiv.Perm (Fin m), StrictMono (x ∘ e) := by
  classical
  let s := Finset.univ.image x
  have hc : s.card = m := by simp [s, Finset.card_image_of_injective _ hx]
  let f : Fin m ≃ s := Equiv.ofBijective (fun i ↦ ⟨x i, by simp [s]⟩) (by
    constructor
    · intro i j h
      exact hx (congrArg Subtype.val h)
    · rintro ⟨y, hy⟩
      obtain ⟨i, _, rfl⟩ := Finset.mem_image.mp hy
      exact ⟨i, rfl⟩)
  let e := (s.orderIsoOfFin hc).toEquiv.trans f.symm
  refine ⟨e, ?_⟩
  intro i j hij
  have hi : x (e i) = (s.orderIsoOfFin hc i : α) := congrArg Subtype.val (f.apply_symm_apply _)
  have hj : x (e j) = (s.orderIsoOfFin hc j : α) := congrArg Subtype.val (f.apply_symm_apply _)
  change x (e i) < x (e j)
  rw [hi, hj]
  exact (s.orderIsoOfFin hc).strictMono hij

/-- A nontrivial annihilating functional on distinct nodes needs more than
`n+1` nodes. This is the Haar property for polynomial evaluation (ATAP 10.1). -/
theorem card_ge_of_nonzero_annihilating_weights {m n : ℕ} {x c : Fin m → ℝ}
    (hx : Function.Injective x) (hc : ∀ i, c i ≠ 0) (hm : 0 < m)
    (hz : ∀ p : ℝ[X], p.degree ≤ n → ∑ i, c i * p.eval (x i) = 0) :
    n + 2 ≤ m := by
  classical
  by_contra h
  let i : Fin m := ⟨0, hm⟩
  let p := Lagrange.basis Finset.univ x i
  have hp : p.degree ≤ n := by
    rw [Lagrange.degree_basis hx.injOn (Finset.mem_univ i), Finset.card_univ, Fintype.card_fin]
    exact_mod_cast (show m - 1 ≤ n by omega)
  have ht := hz p hp
  have he : ∑ j, c j * p.eval (x j) = c i := by
    rw [Finset.sum_eq_single i]
    · simp [p, Lagrange.eval_basis_self hx.injOn]
    · intro j _ hji
      simp [p, Lagrange.eval_basis_of_ne hji.symm]
    · simp
  exact hc i (he.symm.trans ht)

/-- At `n+2` increasing nodes, nonzero weights annihilating degree-`n`
polynomials have opposite signs at adjacent nodes (ATAP 10.1). -/
theorem mul_weights_neg_of_annihilating {n : ℕ} {x c : Fin (n + 2) → ℝ}
    (hx : StrictMono x) (hc : ∀ i, c i ≠ 0)
    (hz : ∀ p : ℝ[X], p.degree ≤ n → ∑ i, c i * p.eval (x i) = 0)
    (i : Fin (n + 1)) : c i.castSucc * c i.succ < 0 := by
  classical
  let s := (Finset.univ.erase i.castSucc).erase i.succ
  let p : ℝ[X] := ∏ j ∈ s, (X - C (x j))
  have hne : i.castSucc ≠ i.succ := ne_of_lt Fin.castSucc_lt_succ
  have hcard : s.card = n := by simp [s, hne.symm]
  have hp : p.degree ≤ n := by
    simp [p, degree_prod, hcard]
  have hzero (j : Fin (n + 2)) (hj : j ∈ s) : p.eval (x j) = 0 := by
    simp only [p, eval_prod, eval_sub, eval_X, eval_C]
    exact Finset.prod_eq_zero hj (sub_self _)
  have htwo : c i.castSucc * p.eval (x i.castSucc) + c i.succ * p.eval (x i.succ) = 0 := by
    have ht := hz p hp
    rw [← Finset.add_sum_erase _ _ (Finset.mem_univ i.castSucc),
      ← Finset.add_sum_erase _ _ (show i.succ ∈ Finset.univ.erase i.castSucc by simp [hne.symm])] at ht
    have hr : ∑ j ∈ s, c j * p.eval (x j) = 0 := Finset.sum_eq_zero (fun j hj ↦ by rw [hzero j hj, mul_zero])
    change _ + (_ + ∑ j ∈ s, c j * p.eval (x j)) = 0 at ht
    rwa [hr, add_zero] at ht
  have hprod : 0 < p.eval (x i.castSucc) * p.eval (x i.succ) := by
    simp only [p, eval_prod, eval_sub, eval_X, eval_C, ← Finset.prod_mul_distrib]
    apply Finset.prod_pos
    intro j hj
    have hj1 : j ≠ i.castSucc := (Finset.mem_erase.mp (Finset.mem_erase.mp hj).2).1
    have hj2 : j ≠ i.succ := (Finset.mem_erase.mp hj).1
    have hji : j < i.castSucc ∨ i.succ < j := by
      simp only [ne_eq, Fin.ext_iff, Fin.val_castSucc, Fin.val_succ] at hj1 hj2
      change (j : ℕ) < (i : ℕ) ∨ (i : ℕ) + 1 < (j : ℕ)
      omega
    rcases hji with h | h
    · exact mul_pos (sub_pos.mpr (hx h))
        (sub_pos.mpr (hx (h.trans Fin.castSucc_lt_succ)))
    · exact mul_pos_of_neg_of_neg
        (sub_neg.mpr (hx (Fin.castSucc_lt_succ.trans h))) (sub_neg.mpr (hx h))
  have hn : c i.castSucc * p.eval (x i.castSucc) ≠ 0 :=
    mul_ne_zero (hc _) (by intro h; rw [h, zero_mul] at hprod; exact (lt_irrefl 0) hprod)
  have hsq : 0 < (c i.castSucc * p.eval (x i.castSucc)) ^ 2 := sq_pos_of_ne_zero hn
  have hh : (c i.castSucc * c i.succ) * (p.eval (x i.castSucc) * p.eval (x i.succ)) < 0 := by
    have hj : c i.succ * p.eval (x i.succ) = -(c i.castSucc * p.eval (x i.castSucc)) := by linarith
    calc
      _ = (c i.castSucc * p.eval (x i.castSucc)) * (c i.succ * p.eval (x i.succ)) := by ring
      _ = -((c i.castSucc * p.eval (x i.castSucc)) ^ 2) := by rw [hj]; ring
      _ < 0 := neg_neg_of_pos hsq
  by_contra h
  exact (not_le_of_gt hh) (mul_nonneg (le_of_not_gt h) hprod.le)