Theory Lagrange_Interpolation

section ‹Lagrange Interpolation›

text ‹This section introduces the function @{term "interpolate"}, which constructs the Lagrange
interpolation polynomials for a given set of points, followed by a theorem of its correctness.›

theory Lagrange_Interpolation
  imports "HOL-Algebra.Polynomial_Divisibility"
begin

text ‹A finite product in a domain is $0$ if and only if at least one factor is. This could be added
to @{theory "HOL-Algebra.FiniteProduct"} or @{theory "HOL-Algebra.Ring"}.›
lemma (in domain) finprod_zero_iff:
  assumes "finite A"
  assumes "a. a  A  f a  carrier R"
  shows "finprod R f A = 𝟬  (x  A. f x = 𝟬)"
  using assms
proof (induct A rule: finite_induct)
  case empty
  then show ?case by simp
next
  case (insert y F)
  moreover have "f  F  carrier R" using insert by blast
  ultimately show ?case by (simp add:integral_iff)
qed

lemma (in ring) poly_of_const_in_carrier:
  assumes "s  carrier R"
  shows "poly_of_const s  carrier (poly_ring R)"
  using poly_of_const_def assms
  by (simp add:univ_poly_carrier[symmetric] polynomial_def)

lemma (in ring) eval_poly_of_const:
  assumes "x  carrier R"
  shows "eval (poly_of_const x) y = x"
  using assms by (simp add:poly_of_const_def)

lemma (in ring) eval_in_carrier_2:
  assumes "x  carrier (poly_ring R)"
  assumes "y  carrier R"
  shows "eval x y  carrier R"
  using eval_in_carrier univ_poly_carrier polynomial_incl assms by blast

lemma (in domain) poly_mult_degree_le_1:
  assumes "x  carrier (poly_ring R)"
  assumes "y  carrier (poly_ring R)"
  shows "degree (x poly_ring Ry)  degree x + degree y"
proof -
  have "degree (x poly_ring Ry) = (if x = []  y = [] then 0 else degree x + degree y)"
    unfolding univ_poly_mult
    by (metis univ_poly_carrier assms(1,2) carrier_is_subring poly_mult_degree_eq)
  thus ?thesis by (metis nat_le_linear zero_le)
qed

lemma (in domain) poly_mult_degree_le:
  assumes "x  carrier (poly_ring R)"
  assumes "y  carrier (poly_ring R)"
  assumes "degree x  n"
  assumes "degree y  m"
  shows "degree (x poly_ring Ry)  n + m"
  using poly_mult_degree_le_1 assms add_mono by force

lemma (in domain) poly_add_degree_le:
  assumes "x  carrier (poly_ring R)" "degree x  n"
  assumes "y  carrier (poly_ring R)" "degree y  n"
  shows "degree (x poly_ring Ry)  n"
  using assms poly_add_degree
  by (metis dual_order.trans max.bounded_iff univ_poly_add)

lemma (in domain) poly_sub_degree_le:
  assumes "x  carrier (poly_ring R)" "degree x  n"
  assumes "y  carrier (poly_ring R)" "degree y  n"
  shows "degree (x poly_ring Ry)  n"
proof -
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto

  show ?thesis
    unfolding a_minus_def
    using assms univ_poly_a_inv_degree carrier_is_subring poly_add_degree_le x.a_inv_closed
    by simp
qed

lemma (in domain) poly_sum_degree_le:
  assumes "finite A"
  assumes "x. x  A  degree (f x)  n"
  assumes "x. x  A  f x  carrier (poly_ring R)"
  shows "degree (finsum (poly_ring R) f A)  n"
  using assms
proof (induct A rule:finite_induct)
  case empty
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto
  show ?case using empty by (simp add:univ_poly_zero)
next
  case (insert x F)
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto
  have a: "degree (f x poly_ring Rfinsum (poly_ring R) f F)  n"
    using insert poly_add_degree_le x.finsum_closed by auto
  show ?case using insert a by auto
qed

definition (in ring) lagrange_basis_polynomial_aux where
  "lagrange_basis_polynomial_aux S =
    (poly_ring Rs  S. X poly_ring R(poly_of_const s))"

lemma (in domain) lagrange_aux_eval:
  assumes "finite S"
  assumes "S  carrier R"
  assumes "x  carrier R"
  shows "(eval (lagrange_basis_polynomial_aux S) x) = (s  S. x  s)"
proof -
  interpret x:ring_hom_cring "poly_ring R" "R" "(λp. eval p x)"
    by (rule eval_cring_hom[OF carrier_is_subring assms(3)])

  have "a. a  S  X poly_ring Rpoly_of_const a  carrier (poly_ring R)"
    by (meson poly_of_const_in_carrier carrier_is_subring assms(2) cring.cring_simprules(4)
        domain_def subsetD univ_poly_is_domain var_closed(1))

  moreover have "s. s  S  eval (X poly_ring Rpoly_of_const s) x = x  s"
    using assms var_closed carrier_is_subring poly_of_const_in_carrier subsetD[OF assms(2)]
    by (simp add:eval_var eval_poly_of_const)

  moreover have "a_minus R x  S  carrier R"
    using assms by blast

  ultimately show ?thesis
    by (simp add:lagrange_basis_polynomial_aux_def x.hom_finprod cong:finprod_cong')
qed

lemma (in domain) lagrange_aux_poly:
  assumes "finite S"
  assumes "S  carrier R"
  shows "lagrange_basis_polynomial_aux S  carrier (poly_ring R)"
proof -
  have a:"subring (carrier R) R"
    using carrier_is_subring assms by blast

  have b: "a. a  S  X poly_ring Rpoly_of_const a  carrier (poly_ring R)"
    by (meson poly_of_const_in_carrier a assms(2) cring.cring_simprules(4) domain_def subsetD
        univ_poly_is_domain var_closed(1))

  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto

  show ?thesis
    using lagrange_basis_polynomial_aux_def b x.finprod_closed[OF Pi_I] by simp
qed

lemma (in domain) poly_prod_degree_le:
  assumes "finite A"
  assumes "x. x  A  f x  carrier (poly_ring R)"
  shows "degree (finprod (poly_ring R) f A)  (x  A. degree (f x))"
  using assms
proof (induct A rule:finite_induct)
  case empty
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto
  show ?case by (simp add:univ_poly_one)
next
  case (insert x F)
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto
  have a:"f  F  carrier (poly_ring R)"
    using insert by blast
  have b:"f x  carrier (poly_ring R)"
    using insert by blast
  have "degree (finprod (poly_ring R) f (insert x F)) = degree (f x poly_ring Rfinprod (poly_ring R) f F)"
    using a b insert by simp
  also have "...  degree (f x) + degree (finprod (poly_ring R) f F)"
    using poly_mult_degree_le x.finprod_closed[OF a] b by auto
  also have "...  degree (f x) + (y  F. degree (f y))"
    using insert(3) a add_mono by auto
  also have "... = (y  (insert x F). degree (f y))" using insert by simp
  finally show ?case by simp
qed

lemma (in domain) lagrange_aux_degree:
  assumes "finite S"
  assumes "S  carrier R"
  shows "degree (lagrange_basis_polynomial_aux S)  card S"
proof -
  interpret x:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto

  have "degree X  1" by (simp add:var_def)
  moreover have "y. y S  degree (poly_of_const y)  1" by (simp add:poly_of_const_def)
  ultimately have a: "y. y S  degree (X poly_ring Rpoly_of_const y)  1"
    by (meson assms(2) in_mono poly_of_const_in_carrier poly_sub_degree_le var_closed[OF carrier_is_subring])

  have b:"y. y  S  (X poly_ring Rpoly_of_const y)  carrier (poly_ring R)"
    by (meson subsetD x.minus_closed var_closed(1)[OF carrier_is_subring] poly_of_const_in_carrier assms(2))

  have "degree (lagrange_basis_polynomial_aux S)  (y  S. degree (X poly_ring Rpoly_of_const y))"
    using lagrange_basis_polynomial_aux_def b poly_prod_degree_le[OF assms(1)] by auto
  also have "...  (y  S. 1)"
    using sum_mono a by force
  also have "... = card S" by simp
  finally show ?thesis by simp
qed

definition (in ring) lagrange_basis_polynomial where
  "lagrange_basis_polynomial S x = lagrange_basis_polynomial_aux S
    poly_ring R(poly_of_const (invR(s  S. x  s)))"

lemma (in field)
  assumes "finite S"
  assumes "S  carrier R"
  assumes "x  carrier R - S"
  shows
    lagrange_one: "eval (lagrange_basis_polynomial S x) x = 𝟭" and
    lagrange_degree: "degree (lagrange_basis_polynomial S x)  card S" and
    lagrange_zero: "s. s  S  eval (lagrange_basis_polynomial S x) s = 𝟬" and
    lagrange_poly: "lagrange_basis_polynomial S x  carrier (poly_ring R)"
proof -
  interpret x:ring_hom_cring "poly_ring R" "R" "(λp. eval p x)"
    using assms carrier_is_subring eval_cring_hom by blast

  define p where "p = lagrange_basis_polynomial_aux S"
  have a:"eval p x = (s  S. x  s)"
    using assms by (simp add:p_def lagrange_aux_eval)

  have b:"p  carrier (poly_ring R)" using assms
    by (simp add:p_def lagrange_aux_poly)

  have "y. y  S  a_minus R x y  carrier R"
    using assms by blast

  hence c:"finprod R (a_minus R x) S  Units R"
    using finprod_closed[OF Pi_I] assms
    by (auto simp add:field_Units finprod_zero_iff)

  have "eval (lagrange_basis_polynomial S x) x =
    (s  S. x  s)  eval (poly_of_const (inv finprod R (a_minus R x) S)) x"
    using poly_of_const_in_carrier Units_inv_closed c p_def[symmetric]
    by (simp add: lagrange_basis_polynomial_def x.hom_mult[OF b] a)
  also have "... =  𝟭"
    using poly_of_const_in_carrier Units_inv_closed c eval_poly_of_const by simp
  finally show "eval (lagrange_basis_polynomial S x) x = 𝟭" by simp

  have "degree (lagrange_basis_polynomial S x)  degree p + degree (poly_of_const (inv finprod R (a_minus R x) S))"
    unfolding lagrange_basis_polynomial_def p_def[symmetric]
    using poly_mult_degree_le[OF b] poly_of_const_in_carrier Units_inv_closed c by auto
  also have "...  card S + 0"
    using add_mono lagrange_aux_degree[OF assms(1) assms(2)] p_def poly_of_const_def by auto
  finally show "degree (lagrange_basis_polynomial S x)  card S" by simp

  show "s. s  S  eval (lagrange_basis_polynomial S x) s = 𝟬"
  proof -
    fix s
    assume d:"s  S"

    interpret s:ring_hom_cring "poly_ring R" "R" "(λp. eval p s)"
      using eval_cring_hom carrier_is_subring assms d by blast

    have "eval p s = finprod R (a_minus R s) S"
      using subsetD[OF assms(2) d] assms
      by (simp add:p_def lagrange_aux_eval)
    also have "... = 𝟬"
      using subsetD[OF assms(2)] d assms by (simp add: finprod_zero_iff)
    finally have "eval p s = 𝟬R⇙" by simp

    moreover have "eval (poly_of_const (inv finprod R (a_minus R x) S)) s  carrier R"
      using s.hom_closed poly_of_const_in_carrier  Units_inv_closed c by blast

    ultimately show "eval (lagrange_basis_polynomial S x) s = 𝟬"
      using poly_of_const_in_carrier Units_inv_closed c
      by (simp add:lagrange_basis_polynomial_def Let_def p_def[symmetric] s.hom_mult[OF b])
  qed

  interpret r:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto

  show "lagrange_basis_polynomial S x  carrier (poly_ring R)"
    using lagrange_basis_polynomial_def p_def[symmetric] poly_of_const_in_carrier Units_inv_closed
      a b c by simp
qed

definition (in ring) interpolate where
  "interpolate S f =
   (poly_ring Rs  S. lagrange_basis_polynomial (S - {s}) s poly_ring R(poly_of_const (f s)))"

text ‹Let @{term "f"} be a function and @{term "S"} be a finite subset of the domain of the field.
Then @{term "interpolate S f"} will return a polynomial with degree less than @{term "card S"}
interpolating @{term "f"} on @{term "S"}.›

theorem (in field)
  assumes "finite S"
  assumes "S  carrier R"
  assumes "f ` S  carrier R"
  shows
    interpolate_poly: "interpolate S f  carrier (poly_ring R)" and
    interpolate_degree: "degree (interpolate S f)  card S - 1" and
    interpolate_eval: "s. s  S  eval (interpolate S f) s = f s"
proof -
  interpret r:cring "poly_ring R"
    using carrier_is_subring domain.univ_poly_is_cring domain_axioms by auto

  have a:"x. x  S  lagrange_basis_polynomial (S - {x}) x  carrier (poly_ring R)"
    by (meson lagrange_poly assms Diff_iff finite_Diff in_mono insertI1 subset_insertI2 subset_insert_iff)

  have b:"x. x  S  f x  carrier R" using assms by blast

  have c:"x. x  S  degree (lagrange_basis_polynomial (S - {x}) x)  card S - 1"
    by (metis (full_types) lagrange_degree  DiffI Diff_insert_absorb assms(1) assms(2)
        card_Diff_singleton finite_insert insert_subset mk_disjoint_insert)

  have d: "x. x  S 
    degree (lagrange_basis_polynomial (S - {x}) x poly_ring Rpoly_of_const (f x))  (card S - 1) + 0"
    using poly_of_const_in_carrier[OF b] poly_mult_degree_le[OF a] c poly_of_const_def by fastforce

  show "interpolate S f  carrier (poly_ring R)"
    using interpolate_def poly_of_const_in_carrier a b by simp

  show "degree (interpolate S f)  card S - 1"
    using poly_sum_degree_le[OF assms(1) d] poly_of_const_in_carrier[OF b] interpolate_def a by simp

  have e:"subring (carrier R) R"
    using carrier_is_subring assms by blast

  show "s. s  S  eval (interpolate S f) s = f s"
  proof -
    fix s
    assume f:"s  S"
    interpret s:ring_hom_cring "poly_ring R" "R" "(λp. eval p s)"
      using eval_cring_hom[OF e] assms f by blast
    have g:"i. i  S 
         eval (lagrange_basis_polynomial (S - {i}) i poly_ring Rpoly_of_const (f i)) s =
         (if s = i then f s else 𝟬)"
    proof -
      fix i
      assume i_in_S: "i  S"
      have "eval (lagrange_basis_polynomial (S - {i}) i poly_ring Rpoly_of_const (f i)) s =
        eval (lagrange_basis_polynomial (S - {i}) i) s  f i"
        using b i_in_S poly_of_const_in_carrier
        by (simp add: s.hom_mult[OF a] eval_poly_of_const)
      also have "... = (if s = i then f s else 𝟬)"
        using b i_in_S poly_of_const_in_carrier assms f
        apply (cases "s=i", simp, subst lagrange_one, auto)
        by (subst lagrange_zero, auto)
      finally show
        "eval (lagrange_basis_polynomial (S - {i}) i poly_ring Rpoly_of_const (f i)) s =
         (if s = i then f s else 𝟬)" by simp
    qed

    have "eval (interpolate S f) s =
    (xS. eval (lagrange_basis_polynomial (S - {x}) x poly_ring Rpoly_of_const (f x)) s)"
      using  poly_of_const_in_carrier[OF b] a e
      by (simp add: interpolate_def s.hom_finsum[OF Pi_I] comp_def)
    also have "... =  (xS. if s = x then f s else 𝟬)"
      using b g by (simp cong: finsum_cong)
    also have "... = f s"
      using finsum_singleton[OF f assms(1)] f assms by auto
    finally show "eval (interpolate S f) s = f s" by simp
  qed
qed

end