section ‹Polynomial Interpolation›
text ‹We combine Newton's, Lagrange's, and Neville-Aitken's interpolation algorithms
to a combined interpolation algorithm which is parametric. This parametric algorithm is then
further extend from fields to also perform interpolation of integer polynomials.
In experiments it is revealed that Newton's algorithm performs better than the one
of Lagrange. Moreover, on the integer numbers, only Newton's algorithm has been
optimized with fast failure capabilities.›
theory Polynomial_Interpolation
imports
Improved_Code_Equations
Newton_Interpolation
Lagrange_Interpolation
Neville_Aitken_Interpolation
begin
datatype interpolation_algorithm = Newton | Lagrange | Neville_Aitken
fun interpolation_poly :: "interpolation_algorithm ⇒ ('a :: field × 'a)list ⇒ 'a poly" where
"interpolation_poly Newton = newton_interpolation_poly"
| "interpolation_poly Lagrange = lagrange_interpolation_poly"
| "interpolation_poly Neville_Aitken = neville_aitken_interpolation_poly"
fun interpolation_poly_int :: "interpolation_algorithm ⇒ (int × int)list ⇒ int poly option" where
"interpolation_poly_int Newton xs_ys = newton_interpolation_poly_int xs_ys"
| "interpolation_poly_int alg xs_ys = (let
rxs_ys = map (λ (x,y). (of_int x, of_int y)) xs_ys;
rp = interpolation_poly alg rxs_ys
in if (∀ x ∈ set (coeffs rp). is_int_rat x) then
Some (map_poly int_of_rat rp) else None)"
lemma interpolation_poly_int_def: "distinct (map fst xs_ys) ⟹
interpolation_poly_int alg xs_ys = (let
rxs_ys = map (λ (x,y). (of_int x, of_int y)) xs_ys;
rp = interpolation_poly alg rxs_ys
in if (∀ x ∈ set (coeffs rp). is_int_rat x) then
Some (map_poly int_of_rat rp) else None)"
by (cases alg, auto simp: newton_interpolation_poly_int)
lemma interpolation_poly: assumes dist: "distinct (map fst xs_ys)"
and p: "p = interpolation_poly alg xs_ys"
and xy: "(x,y) ∈ set xs_ys"
shows "poly p x = y"
proof (cases alg)
case Newton
thus ?thesis using newton_interpolation_poly[OF dist _ xy] p by simp
next
case Lagrange
thus ?thesis using lagrange_interpolation_poly[OF dist _ xy] p by simp
next
case Neville_Aitken
thus ?thesis using neville_aitken_interpolation_poly[OF dist _ xy] p by simp
qed
lemma degree_interpolation_poly:
shows "degree (interpolation_poly alg xs_ys) ≤ length xs_ys - 1"
using degree_lagrange_interpolation_poly[of xs_ys]
degree_newton_interpolation_poly[of xs_ys]
degree_neville_aitken_interpolation_poly[of xs_ys]
by (cases alg, auto)
lemma uniqueness_of_interpolation: fixes p :: "'a :: idom poly"
assumes cS: "card S = Suc n"
and "degree p ≤ n" and "degree q ≤ n" and
id: "⋀ x. x ∈ S ⟹ poly p x = poly q x"
shows "p = q"
proof -
define f where "f = p - q"
let ?R = "{x. poly f x = 0}"
have sub: "S ⊆ ?R" unfolding f_def using id by auto
show ?thesis
proof (cases "f = 0")
case True thus ?thesis unfolding f_def by simp
next
case False note f = this
let ?R = "{x. poly f x = 0}"
from poly_roots_finite[OF f] have "finite ?R" .
from card_mono[OF this sub] poly_roots_degree[OF f]
have "Suc n ≤ degree f" unfolding cS by auto
also have "… ≤ n" unfolding f_def
by (rule degree_diff_le, insert assms, auto)
finally show ?thesis by auto
qed
qed
lemma uniqueness_of_interpolation_point_list: fixes p :: "'a :: idom poly"
assumes dist: "distinct (map fst xs_ys)"
and p: "⋀ x y. (x,y) ∈ set xs_ys ⟹ poly p x = y" "degree p < length xs_ys"
and q: "⋀ x y. (x,y) ∈ set xs_ys ⟹ poly q x = y" "degree q < length xs_ys"
shows "p = q"
proof -
let ?xs = "map fst xs_ys"
from q obtain n where len: "length xs_ys = Suc n" and dq: "degree q ≤ n" by (cases xs_ys, auto)
from p have dp: "degree p ≤ n" unfolding len by auto
from dist have card: "card (set ?xs) = Suc n" unfolding len[symmetric]
using distinct_card by fastforce
show "p = q"
proof (rule uniqueness_of_interpolation[OF card dp dq])
fix x
assume "x ∈ set ?xs"
then obtain y where "(x,y) ∈ set xs_ys" by auto
from p(1)[OF this] q(1)[OF this] show "poly p x = poly q x" by simp
qed
qed
lemma exactly_one_poly_interpolation: assumes xs: "xs_ys ≠ []" and dist: "distinct (map fst xs_ys)"
shows "∃! p. degree p < length xs_ys ∧ (∀ x y. (x,y) ∈ set xs_ys ⟶ poly p x = (y :: 'a :: field))"
proof -
let ?alg = "undefined"
let ?p = "interpolation_poly ?alg xs_ys"
note inter = interpolation_poly[OF dist refl]
show ?thesis
proof (rule ex1I[of _ ?p], intro conjI allI impI)
show dp: "degree ?p < length xs_ys" using degree_interpolation_poly[of ?alg xs_ys] xs by (cases xs_ys, auto)
show "⋀x y. (x, y) ∈ set xs_ys ⟹ poly (interpolation_poly ?alg xs_ys) x = y"
by (rule inter)
fix q
assume q: "degree q < length xs_ys ∧ (∀x y. (x, y) ∈ set xs_ys ⟶ poly q x = y)"
show "q = ?p"
by (rule uniqueness_of_interpolation_point_list[OF dist _ _ inter dp], insert q, auto)
qed
qed
lemma interpolation_poly_int_Some: assumes dist': "distinct (map fst xs_ys)"
and p: "interpolation_poly_int alg xs_ys = Some p"
shows "⋀ x y. (x,y) ∈ set xs_ys ⟹ poly p x = y" "degree p ≤ length xs_ys - 1"
proof -
let ?r = "rat_of_int"
define rxs_ys where "rxs_ys = map (λ(x, y). (?r x, ?r y)) xs_ys"
have dist: "distinct (map fst rxs_ys)" using dist' unfolding distinct_map rxs_ys_def inj_on_def by force
obtain rp where rp: "rp = interpolation_poly alg rxs_ys" by blast
from p[unfolded interpolation_poly_int_def[OF dist'] Let_def, folded rxs_ys_def rp]
have p: "p = map_poly int_of_rat rp" and ball: "Ball (set (coeffs rp)) is_int_rat"
by (auto split: if_splits)
have id: "rp = map_poly ?r p" unfolding p
by (rule sym, subst map_poly_map_poly, force, rule map_poly_idI, insert ball, auto)
note inter = interpolation_poly[OF dist rp]
{
fix x y
assume "(x,y) ∈ set xs_ys"
hence "(?r x, ?r y) ∈ set rxs_ys" unfolding rxs_ys_def by auto
from inter[OF this] have "poly rp (?r x) = ?r y" by auto
from this[unfolded id of_int_hom.poly_map_poly] show "poly p x = y" by auto
}
show "degree p ≤ length xs_ys - 1" using degree_interpolation_poly[of alg rxs_ys, folded rp]
unfolding id rxs_ys_def by simp
qed
lemma interpolation_poly_int_None: assumes dist: "distinct (map fst xs_ys)"
and p: "interpolation_poly_int alg xs_ys = None"
and q: "⋀ x y. (x,y) ∈ set xs_ys ⟹ poly q x = y"
and dq: "degree q < length xs_ys"
shows False
proof -
let ?r = "rat_of_int"
let ?rp = "map_poly ?r"
define rxs_ys where "rxs_ys = map (λ(x, y). (?r x, ?r y)) xs_ys"
have dist': "distinct (map fst rxs_ys)" using dist unfolding distinct_map rxs_ys_def inj_on_def by force
obtain rp where rp: "rp = interpolation_poly alg rxs_ys" by blast
note degrp = degree_interpolation_poly[of alg rxs_ys, folded rp]
from q have q': "⋀ x y. (x,y) ∈ set rxs_ys ⟹ poly (?rp q) x = y" unfolding rxs_ys_def
by auto
have [simp]: "degree (?rp q) = degree q" by simp
have id: "rp = ?rp q"
by (rule uniqueness_of_interpolation_point_list[OF dist' interpolation_poly[OF dist' rp]],
insert q' dq degrp, auto simp: rxs_ys_def)
from p[unfolded interpolation_poly_int_def[OF dist] Let_def, folded rxs_ys_def rp]
have "∃ c ∈ set (coeffs rp). c ∉ ℤ" by (auto split: if_splits)
from this[unfolded id] show False by auto
qed
lemmas newton_interpolation_poly_int_Some =
interpolation_poly_int_Some[where alg = Newton, unfolded interpolation_poly_int.simps]
lemmas newton_interpolation_poly_int_None =
interpolation_poly_int_None[where alg = Newton, unfolded interpolation_poly_int.simps]
text ‹We can also use Newton's improved algorithm for integer polynomials to show that
there is no polynomial $p$ over the integers such that $p(0) = 0$ and $p(2) = 1$.
The reason is that the intermediate result for computing the linear interpolant for these
two point fails, and so adding further points (which corresponds to increasing the degree)
will also fail. Of course, this can be generalized, showing that whenever you cannot
interpolate a set of $n$ points with an integer polynomial of degree $n-1$, then you
cannot interpolate this set of points with any integer polynomial. However, we did not
formally prove this more general fact.›
lemma impossible_p_0_is_0_and_p_2_is_1: "¬ (∃ p. poly p 0 = 0 ∧ poly p 2 = (1 :: int))"
proof
assume "∃ p. poly p 0 = 0 ∧ poly p 2 = (1 :: int)"
then obtain p where p: "poly p 0 = 0" "poly p 2 = (1 :: int)" by auto
define xs_ys where "xs_ys = map (λ i. (int i, poly p (int i))) [ 3 ..< 3 + degree p]"
let ?l = "λ xs. (0,0) # (2 :: int,1 :: int) # xs"
let ?xs_ys = "?l xs_ys"
define list where "list = map fst ?xs_ys"
have dist: "distinct (map fst ?xs_ys)" unfolding xs_ys_def by (auto simp: o_def distinct_map inj_on_def)
have p: "⋀ x y. (x,y) ∈ set ?xs_ys ⟹ poly p x = y" unfolding xs_ys_def using p by auto
have deg: "degree p < length ?xs_ys" unfolding xs_ys_def by simp
have "newton_coefficients_main_int list (rev (map snd ?xs_ys)) (rev (map fst ?xs_ys)) = None"
proof (induct xs_ys rule: rev_induct)
case Nil
show ?case unfolding list_def by (simp add: divmod_int_def)
next
case (snoc xy xs_ys) note IH = this
obtain x y where xy: "xy = (x,y)" by force
show ?case
proof (cases xs_ys rule: rev_cases)
case Nil
show ?thesis unfolding Nil xy
by (simp add: list_def divmod_int_def)
next
case (snoc xs_ys' xy')
obtain x' y' where xy': "xy' = (x',y')" by force
show ?thesis using IH unfolding xy' snoc xy by simp
qed
qed
hence newton: "newton_interpolation_poly_int ?xs_ys = None"
unfolding newton_interpolation_poly_int_def Let_def newton_poly_impl_int_def
Newton_Interpolation.newton_coefficients_int_def list_def by simp
from newton_interpolation_poly_int_None[OF dist newton p deg]
show False .
qed
end