feat(linear_algebra/lagrange): Lagrange interpolation (#2764)

Fixes #2600.

src/data/polynomial.lean

@@ -1374,6 +1374,9 @@ is_ring_hom.map_neg _
 @[simp] lemma degree_neg (p : polynomial R) : degree (-p) = degree p :=
 by unfold degree; rw support_neg
 
+lemma degree_sub_le (p q : polynomial R) : degree (p - q) ≤ max (degree p) (degree q) :=
+degree_neg q ▸ degree_add_le p (-q)
+
 @[simp] lemma nat_degree_neg (p : polynomial R) : nat_degree (-p) = nat_degree p :=
 by simp [nat_degree]
 

src/linear_algebra/lagrange.lean

@@ -0,0 +1,157 @@
+/-
+Copyright (c) 2020 Kenny Lau. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Author: Kenny Lau.
+-/
+
+import ring_theory.polynomial algebra.big_operators
+
+/-!
+# Lagrange interpolation
+
+## Main definitions
+
+* `lagrange.basis s x` where `s : finset F` and `x : F`: the Lagrange basis polynomial
+  that evaluates to `1` at `x` and `0` at other elements of `s`.
+* `lagrange.interpolate s f` where `s : finset F` and `f : s → F`: the Lagrange interpolant
+  that evaluates to `f x` at `x` for `x ∈ s`.
+
+-/
+
+noncomputable theory
+open_locale big_operators classical
+
+universe u
+
+namespace lagrange
+
+variables {F : Type u} [decidable_eq F] [field F] (s : finset F)
+variables {F' : Type u} [field F'] (s' : finset F')
+
+open polynomial
+
+/-- Lagrange basis polynomials that evaluate to 1 at `x` and 0 at other elements of `s`. -/
+def basis (x : F) : polynomial F :=
+∏ y in s.erase x, C (x - y)⁻¹ * (X - C y)
+
+@[simp] theorem basis_empty (x : F) : basis ∅ x = 1 :=
+rfl
+
+@[simp] theorem eval_basis_self (x : F) : (basis s x).eval x = 1 :=
+begin
+  rw [basis, ← (s.erase x).prod_hom (eval x), finset.prod_eq_one],
+  intros y hy, simp_rw [eval_mul, eval_sub, eval_C, eval_X],
+  exact inv_mul_cancel (sub_ne_zero_of_ne (finset.ne_of_mem_erase hy).symm)
+end
+
+@[simp] theorem eval_basis_ne (x y : F) (h1 : y ∈ s) (h2 : y ≠ x) : (basis s x).eval y = 0 :=
+begin
+  rw [basis, ← (s.erase x).prod_hom (eval y), finset.prod_eq_zero (finset.mem_erase.2 ⟨h2, h1⟩)],
+  simp_rw [eval_mul, eval_sub, eval_C, eval_X, sub_self, mul_zero]
+end
+
+theorem eval_basis (x y : F) (h : y ∈ s) : (basis s x).eval y = if y = x then 1 else 0 :=
+by { split_ifs with H, { subst H, apply eval_basis_self }, { exact eval_basis_ne s x y h H } }
+
+@[simp] theorem nat_degree_basis (x : F) (hx : x ∈ s) : (basis s x).nat_degree = s.card - 1 :=
+begin
+  unfold basis, generalize hsx : s.erase x = sx,
+  have : x ∉ sx := hsx ▸ finset.not_mem_erase x s,
+  rw [← finset.insert_erase hx, hsx, finset.card_insert_of_not_mem this, nat.add_sub_cancel],
+  clear hx hsx s, revert this, apply sx.induction_on,
+  { intros hx, rw [finset.prod_empty, nat_degree_one], refl },
+  { intros y s hys ih hx, rw [finset.mem_insert, not_or_distrib] at hx,
+    have h1 : C (x - y)⁻¹ ≠ C 0 := λ h, hx.1 (eq_of_sub_eq_zero $ inv_eq_zero.1 $ C_inj.1 h),
+    have h2 : X ^ 1 - C y ≠ 0 := by convert X_pow_sub_C_ne_zero zero_lt_one _,
+    rw C_0 at h1, rw pow_one at h2,
+    rw [finset.prod_insert hys, nat_degree_mul_eq (mul_ne_zero h1 h2), ih hx.2,
+        finset.card_insert_of_not_mem hys, nat_degree_mul_eq h1 h2,
+        nat_degree_C, zero_add, nat_degree, degree_X_sub_C, add_comm], refl,
+    rw [ne, finset.prod_eq_zero_iff], rintro ⟨z, hzs, hz⟩,
+    rw mul_eq_zero at hz, cases hz with hz hz,
+    { rw [← C_0, C_inj, inv_eq_zero, sub_eq_zero] at hz, exact hx.2 (hz.symm ▸ hzs) },
+    { rw ← pow_one (X : polynomial F) at hz, exact X_pow_sub_C_ne_zero zero_lt_one _ hz } }
+end
+
+variables (f : (↑s : set F) → F)
+
+/-- Lagrange interpolation: given a finset `s` and a function `f : s → F`,
+`interpolate s f` is the unique polynomial of degree `< s.card`
+that takes value `f x` on all `x` in `s`. -/
+def interpolate : polynomial F :=
+∑ x in s.attach, C (f x) * basis s x
+
+@[simp] theorem interpolate_empty (f) : interpolate (∅ : finset F) f = 0 :=
+rfl
+
+@[simp] theorem eval_interpolate (x) (H : x ∈ s) : eval x (interpolate s f) = f ⟨x, H⟩ :=
+begin
+  rw [interpolate, ← finset.sum_hom _ (eval x), finset.sum_eq_single (⟨x, H⟩ : { x // x ∈ s })],
+  { rw [eval_mul, eval_C, subtype.coe_mk, eval_basis_self, mul_one] },
+  { rintros ⟨y, hy⟩ _ hyx, rw [eval_mul, subtype.coe_mk, eval_basis_ne s y x H, mul_zero],
+    { rintros rfl, exact hyx rfl } },
+  { intro h, exact absurd (finset.mem_attach _ _) h }
+end
+
+theorem degree_interpolate_lt : (interpolate s f).degree < s.card :=
+if H : s = ∅ then by { subst H, rw [interpolate_empty, degree_zero], exact with_bot.bot_lt_coe _ }
+else lt_of_le_of_lt (degree_sum_le _ _) $ (finset.sup_lt_iff $ with_bot.bot_lt_coe s.card).2 $ λ b _,
+calc  (C (f b) * basis s b).degree
+    ≤ (C (f b)).degree + (basis s b).degree : degree_mul_le _ _
+... ≤ 0 + (basis s b).degree : add_le_add_right' degree_C_le
+... = (basis s b).degree : zero_add _
+... ≤ (basis s b).nat_degree : degree_le_nat_degree
+... = (s.card - 1 : ℕ) : by { rw nat_degree_basis s b b.2 }
+... < s.card : with_bot.coe_lt_coe.2 (nat.pred_lt $ mt finset.card_eq_zero.1 H)
+
+/-- Linear version of `interpolate`. -/
+def linterpolate : ((↑s : set F) → F) →ₗ[F] polynomial F :=
+{ to_fun := interpolate s,
+  add := λ f g, by { simp_rw [interpolate, ← finset.sum_add_distrib, ← add_mul, ← C_add], refl },
+  smul := λ c f, by { simp_rw [interpolate, finset.smul_sum, C_mul', smul_smul], refl } }
+
+@[simp] lemma interpolate_add (f g) : interpolate s (f + g) = interpolate s f + interpolate s g :=
+(linterpolate s).map_add f g
+
+@[simp] lemma interpolate_zero : interpolate s 0 = 0 :=
+(linterpolate s).map_zero
+
+@[simp] lemma interpolate_neg (f) : interpolate s (-f) = -interpolate s f :=
+(linterpolate s).map_neg f
+
+@[simp] lemma interpolate_sub (f g) : interpolate s (f - g) = interpolate s f - interpolate s g :=
+(linterpolate s).map_sub f g
+
+@[simp] lemma interpolate_smul (c : F) (f) : interpolate s (c • f) = c • interpolate s f :=
+(linterpolate s).map_smul c f
+
+theorem eq_zero_of_eval_eq_zero {f : polynomial F'} (hf1 : f.degree < s'.card)
+  (hf2 : ∀ x ∈ s', f.eval x = 0) : f = 0 :=
+by_contradiction $ λ hf3, not_le_of_lt hf1 $
+calc  (s'.card : with_bot ℕ)
+    ≤ f.roots.card : with_bot.coe_le_coe.2 $ finset.card_le_of_subset $ λ x hx,
+        (mem_roots hf3).2 $ hf2 x hx
+... ≤ f.degree : card_roots hf3
+
+theorem eq_of_eval_eq {f g : polynomial F'} (hf : f.degree < s'.card) (hg : g.degree < s'.card)
+  (hfg : ∀ x ∈ s', f.eval x = g.eval x) : f = g :=
+eq_of_sub_eq_zero $ eq_zero_of_eval_eq_zero s'
+  (lt_of_le_of_lt (degree_sub_le f g) $ max_lt hf hg)
+  (λ x hx, by rw [eval_sub, hfg x hx, sub_self])
+
+theorem eq_interpolate (f : polynomial F) (hf : f.degree < s.card) :
+  interpolate s (λ x, f.eval x) = f :=
+eq_of_eval_eq s (degree_interpolate_lt s _) hf $ λ x hx, eval_interpolate s _ x hx
+
+/-- Lagrange interpolation induces isomorphism between functions from `s` and polynomials
+of degree less than `s.card`. -/
+def fun_equiv_degree_lt : degree_lt F s.card ≃ₗ[F] ((↑s : set F) → F) :=
+{ to_fun := λ f x, f.1.eval x,
+  add := λ f g, funext $ λ x, eval_add,
+  smul := λ c f, funext $ λ x, by { rw [pi.smul_apply, smul_eq_mul, ← @eval_C F c _ x,
+      ← eval_mul, eval_C, C_mul'], refl },
+  inv_fun := λ f, ⟨interpolate s f, mem_degree_lt.2 $ degree_interpolate_lt s f⟩,
+  left_inv := λ f, subtype.eq $ eq_interpolate s f $ mem_degree_lt.1 f.2,
+  right_inv := λ f, funext $ λ ⟨x, hx⟩, eval_interpolate s f x hx }
+
+end lagrange

src/ring_theory/polynomial.lean

@@ -23,13 +23,17 @@ variables (R : Type u) [comm_ring R]
 def degree_le (n : with_bot ℕ) : submodule R (polynomial R) :=
 ⨅ k : ℕ, ⨅ h : ↑k > n, (lcoeff R k).ker
 
+/-- The `R`-submodule of `R[X]` consisting of polynomials of degree < `n`. -/
+def degree_lt (n : ℕ) : submodule R (polynomial R) :=
+⨅ k : ℕ, ⨅ h : k ≥ n, (lcoeff R k).ker
+
 variable {R}
 
 theorem mem_degree_le {n : with_bot ℕ} {f : polynomial R} :
   f ∈ degree_le R n ↔ degree f ≤ n :=
 by simp only [degree_le, submodule.mem_infi, degree_le_iff_coeff_zero, linear_map.mem_ker]; refl
 
-theorem degree_le_mono {m n : with_bot ℕ} (H : m ≤ n) :
+@[mono] theorem degree_le_mono {m n : with_bot ℕ} (H : m ≤ n) :
   degree_le R m ≤ degree_le R n :=
 λ f hf, mem_degree_le.2 (le_trans (mem_degree_le.1 hf) H)
 
@@ -50,6 +54,33 @@ begin
   apply le_trans (degree_X_pow_le _) (with_bot.coe_le_coe.2 $ nat.le_of_lt_succ $ finset.mem_range.1 hk)
 end
 
+theorem mem_degree_lt {n : ℕ} {f : polynomial R} :
+  f ∈ degree_lt R n ↔ degree f < n :=
+by { simp_rw [degree_lt, submodule.mem_infi, linear_map.mem_ker, degree,
+    finset.sup_lt_iff (with_bot.bot_lt_coe n), finsupp.mem_support_iff, with_bot.some_eq_coe,
+    with_bot.coe_lt_coe, lt_iff_not_ge', ne, not_imp_not], refl }
+
+@[mono] theorem degree_lt_mono {m n : ℕ} (H : m ≤ n) :
+  degree_lt R m ≤ degree_lt R n :=
+λ f hf, mem_degree_lt.2 (lt_of_lt_of_le (mem_degree_lt.1 hf) $ with_bot.coe_le_coe.2 H)
+
+theorem degree_lt_eq_span_X_pow {n : ℕ} :
+  degree_lt R n = submodule.span R ↑((finset.range n).image (λ n, X^n) : finset (polynomial R)) :=
+begin
+  apply le_antisymm,
+  { intros p hp, replace hp := mem_degree_lt.1 hp,
+    rw [← finsupp.sum_single p, finsupp.sum],
+    refine submodule.sum_mem _ (λ k hk, _),
+    show monomial _ _ ∈ _,
+    have := with_bot.coe_lt_coe.1 ((finset.sup_lt_iff $ with_bot.bot_lt_coe n).1 hp k hk),
+    rw [single_eq_C_mul_X, C_mul'],
+    refine submodule.smul_mem _ _ (submodule.subset_span $ finset.mem_coe.2 $
+      finset.mem_image.2 ⟨_, finset.mem_range.2 this, rfl⟩) },
+  rw [submodule.span_le, finset.coe_image, set.image_subset_iff],
+  intros k hk, apply mem_degree_lt.2,
+  exact lt_of_le_of_lt (degree_X_pow_le _) (with_bot.coe_lt_coe.2 $ finset.mem_range.1 hk)
+end
+
 /-- Given a polynomial, return the polynomial whose coefficients are in
 the ring closure of the original coefficients. -/
 def restriction (p : polynomial R) : polynomial (ring.closure (↑p.frange : set R)) :=