section ‹Neville Aitken Interpolation›
text ‹We prove soundness of Neville-Aitken's polynomial interpolation algorithm
using the recursive formula directly. We further provide an implementation
which avoids the exponential branching in the recursion.›
theory Neville_Aitken_Interpolation
imports
"HOL-Computational_Algebra.Polynomial"
begin
context
fixes x :: "nat ⇒ 'a :: field"
and f :: "nat ⇒ 'a"
begin
private definition X :: "nat ⇒ 'a poly" where [code_unfold]: "X i = [:-x i, 1:]"
function neville_aitken_main :: "nat ⇒ nat ⇒ 'a poly" where
"neville_aitken_main i j = (if i < j then
(smult (inverse (x j - x i)) (X i * neville_aitken_main (i + 1) j -
X j * neville_aitken_main i (j - 1)))
else [:f i:])"
by pat_completeness auto
termination by (relation "measure (λ (i,j). j - i)", auto)
definition neville_aitken :: "nat ⇒ 'a poly" where
"neville_aitken = neville_aitken_main 0"
declare neville_aitken_main.simps[simp del]
lemma neville_aitken_main: assumes dist: "⋀ i j. i < j ⟹ j ≤ n ⟹ x i ≠ x j"
shows "i ≤ k ⟹ k ≤ j ⟹ j ≤ n ⟹ poly (neville_aitken_main i j) (x k) = (f k)"
proof (induct i j arbitrary: k rule: neville_aitken_main.induct)
case (1 i j k)
note neville_aitken_main.simps[of i j, simp]
show ?case
proof (cases "i < j")
case False
with 1(3-) have "k = i" by auto
with False show ?thesis by auto
next
case True note ij = this
from dist[OF True 1(5)] have diff: "x i ≠ x j" by auto
from True have id: "neville_aitken_main i j =
(smult (inverse (x j - x i)) (X i * neville_aitken_main (i + 1) j - X j
* neville_aitken_main i (j - 1)))" by simp
note IH = 1(1-2)[OF True]
show ?thesis
proof (cases "k = i")
case True
show ?thesis unfolding id True poly_smult using IH(2)[of i] ij 1(3-) diff
by (simp add: X_def field_simps)
next
case False note ki = this
show ?thesis
proof (cases "k = j")
case True
show ?thesis unfolding id True poly_smult using IH(1)[of j] ij 1(3-) diff
by (simp add: X_def field_simps)
next
case False
with ki show ?thesis unfolding id poly_smult using IH(1-2)[of k] ij 1(3-) diff
by (simp add: X_def field_simps)
qed
qed
qed
qed
lemma degree_neville_aitken_main: "degree (neville_aitken_main i j) ≤ j - i"
proof (induct i j rule: neville_aitken_main.induct)
case (1 i j)
note simp = neville_aitken_main.simps[of i j]
show ?case
proof (cases "i < j")
case False
thus ?thesis unfolding simp by simp
next
case True
note IH = 1[OF this]
let ?n = neville_aitken_main
have X: "⋀ i. degree (X i) = Suc 0" unfolding X_def by auto
have "degree (X i * ?n (i + 1) j) ≤ Suc (degree (?n (i+1) j))"
by (rule order.trans[OF degree_mult_le], simp add: X)
also have "… ≤ Suc (j - (i+1))" using IH(1) by simp
finally have 1: "degree (X i * ?n (i + 1) j) ≤ j - i" using True by auto
have "degree (X j * ?n i (j - 1)) ≤ Suc (degree (?n i (j - 1)))"
by (rule order.trans[OF degree_mult_le], simp add: X)
also have "… ≤ Suc ((j - 1) - i)" using IH(2) by simp
finally have 2: "degree (X j * ?n i (j - 1)) ≤ j - i" using True by auto
have id: "?n i j = smult (inverse (x j - x i))
(X i * ?n (i + 1) j - X j * ?n i (j - 1))" unfolding simp using True by simp
have "degree (?n i j) ≤ degree (X i * ?n (i + 1) j - X j * ?n i (j - 1))"
unfolding id by simp
also have "… ≤ max (degree (X i * ?n (i + 1) j)) (degree (X j * ?n i (j - 1)))"
by (rule degree_diff_le_max)
also have "… ≤ j - i" using 1 2 by auto
finally show ?thesis .
qed
qed
lemma degree_neville_aitken: "degree (neville_aitken n) ≤ n"
unfolding neville_aitken_def using degree_neville_aitken_main[of 0 n] by simp
fun neville_aitken_merge :: "('a × 'a × 'a poly) list ⇒ ('a × 'a × 'a poly) list" where
"neville_aitken_merge ((xi,xj,p_ij) # (xsi,xsj,p_sisj) # rest) =
(xi,xsj, smult (inverse (xsj - xi)) ([:-xi,1:] * p_sisj
+ [:xsj,-1:] * p_ij)) # neville_aitken_merge ((xsi,xsj,p_sisj) # rest)"
| "neville_aitken_merge [_] = []"
| "neville_aitken_merge [] = []"
lemma length_neville_aitken_merge[termination_simp]: "length (neville_aitken_merge xs) = length xs - 1"
by (induct xs rule: neville_aitken_merge.induct, auto)
fun neville_aitken_impl_main :: "('a × 'a × 'a poly) list ⇒ 'a poly" where
"neville_aitken_impl_main (e1 # e2 # es) =
neville_aitken_impl_main (neville_aitken_merge (e1 # e2 # es))"
| "neville_aitken_impl_main [(_,_,p)] = p"
| "neville_aitken_impl_main [] = 0"
lemma neville_aitken_merge:
"xs = map (λ i. (x i, x (i + j), neville_aitken_main i (i + j))) [l ..< Suc (l + k)]
⟹ neville_aitken_merge xs
= (map (λ i. (x i, x (i + Suc j), neville_aitken_main i (i + Suc j))) [l ..< l + k])"
proof (induct xs arbitrary: l k rule: neville_aitken_merge.induct)
case (1 xi xj p_ij xsi xsj p_sisj rest l k)
let ?n = neville_aitken_main
let ?f = "λ j i. (x i, x (i + j), ?n i (i + j))"
define f where "f = ?f"
let ?map = "λ j. map (?f j)"
note res = 1(2)
from arg_cong[OF res, of length] obtain kk where k: "k = Suc kk" by (cases k, auto)
hence id: "[l..<Suc (l + k)] = l # [Suc l ..< Suc (Suc l + kk)]"
by (simp add: upt_rec)
from res[unfolded id] have id2: "(xsi, xsj, p_sisj) # rest =
?map j [Suc l..< Suc (Suc l + kk)]"
and id3: "xi = x l" "xj = x (l + j)" "p_ij = ?n l (l + j)"
"xsi = x (Suc l)" "xsj = x (Suc (l + j))" "p_sisj = ?n (Suc l) (Suc (l + j))"
by (auto simp: upt_rec)
note IH = 1(1)[OF id2]
have X: "[:x (Suc (l + j)), - 1:] = - X (Suc l + j)" unfolding X_def by simp
have id4: "(xi, xsj, smult (inverse (xsj - xi)) ([:- xi, 1:] * p_sisj +
[:xsj, - 1:] * p_ij)) = (x l, x (l + Suc j), ?n l (l + Suc j))"
unfolding id3 neville_aitken_main.simps[of l "l + Suc j"]
X_def[symmetric] X by simp
have id5: "[l..<l + k] = l # [Suc l ..< Suc l + kk]" unfolding k
by (simp add: upt_rec)
show ?case unfolding neville_aitken_merge.simps IH id4
unfolding id5 by simp
qed auto
lemma neville_aitken_impl_main:
"xs = map (λ i. (x i, x (i + j), neville_aitken_main i (i + j))) [l ..< Suc (l + k)]
⟹ neville_aitken_impl_main xs = neville_aitken_main l (l + j + k)"
proof (induct xs arbitrary: l k j rule: neville_aitken_impl_main.induct)
case (1 e1 e2 es l k j)
note res = 1(2)
from res obtain kk where k: "k = Suc kk" by (cases k, auto)
hence id1: "l + k = Suc (l + kk)" by auto
show ?case unfolding neville_aitken_impl_main.simps 1(1)[OF neville_aitken_merge[OF 1(2), unfolded id1]]
by (simp add: k)
qed auto
lemma neville_aitken_impl:
"xs = map (λ i. (x i, x i, [:f i:])) [0 ..< Suc k]
⟹ neville_aitken_impl_main xs = neville_aitken k"
unfolding neville_aitken_def using neville_aitken_impl_main[of xs 0 0 k]
by (simp add: neville_aitken_main.simps)
end
lemma neville_aitken: assumes "⋀ i j. i < j ⟹ j ≤ n ⟹ x i ≠ x j"
shows "j ≤ n ⟹ poly (neville_aitken x f n) (x j) = (f j)"
unfolding neville_aitken_def
by (rule neville_aitken_main[OF assms, of n], auto)
definition neville_aitken_interpolation_poly :: "('a :: field × 'a)list ⇒ 'a poly" where
"neville_aitken_interpolation_poly x_fs = (let
start = map (λ (xi,fi). (xi,xi,[:fi:])) x_fs in
neville_aitken_impl_main start)"
lemma neville_aitken_interpolation_impl: assumes "x_fs ≠ []"
shows "neville_aitken_interpolation_poly x_fs =
neville_aitken (λ i. fst (x_fs ! i)) (λ i. snd (x_fs ! i)) (length x_fs - 1)"
proof -
from assms have id: "Suc (length x_fs - 1) = length x_fs" by auto
show ?thesis
unfolding neville_aitken_interpolation_poly_def Let_def
by (rule neville_aitken_impl, unfold id, rule nth_equalityI, auto split: prod.splits)
qed
lemma neville_aitken_interpolation_poly: assumes dist: "distinct (map fst xs_ys)"
and p: "p = neville_aitken_interpolation_poly xs_ys"
and xy: "(x,y) ∈ set xs_ys"
shows "poly p x = y"
proof -
have p: "p = neville_aitken (λ i. fst (xs_ys ! i)) (λ i. snd (xs_ys ! i)) (length xs_ys - 1)"
unfolding p
by (rule neville_aitken_interpolation_impl, insert xy, auto)
from xy obtain i where i: "i < length xs_ys" and x: "x = fst (xs_ys ! i)" and y: "y = snd (xs_ys ! i)"
unfolding set_conv_nth by (metis fst_conv in_set_conv_nth snd_conv xy)
show ?thesis unfolding p x y
proof (rule neville_aitken)
fix i j
show "i < j ⟹ j ≤ length xs_ys - 1 ⟹ fst (xs_ys ! i) ≠ fst (xs_ys ! j)" using dist
by (metis (mono_tags, lifting) One_nat_def diff_less dual_order.strict_trans2 length_map
length_pos_if_in_set lessI less_or_eq_imp_le neq_iff nth_eq_iff_index_eq nth_map xy)
qed (insert i, auto)
qed
lemma degree_neville_aitken_interpolation_poly:
shows "degree (neville_aitken_interpolation_poly xs_ys) ≤ length xs_ys - 1"
proof (cases "length xs_ys")
case 0
hence id: "xs_ys = []" by (cases xs_ys, auto)
show ?thesis unfolding id neville_aitken_interpolation_poly_def Let_def by simp
next
case (Suc nn)
have id: "neville_aitken_interpolation_poly xs_ys =
neville_aitken (λ i. fst (xs_ys ! i)) (λ i. snd (xs_ys ! i)) (length xs_ys - 1)"
by (rule neville_aitken_interpolation_impl, insert Suc, auto)
show ?thesis unfolding id by (rule degree_neville_aitken)
qed
end