Theory Newton_Interpolation

(*  
    Author:      René Thiemann 
                 Akihisa Yamada
    License:     BSD
*)
section ‹Newton Interpolation›

text ‹We proved the soundness of the Newton interpolation, i.e., a method to interpolate a polynomial $p$
  from a list of points $(x_1,p(x_1)), (x_2, p(x_2)), \ldots$. In experiments it performs
  much faster than the Lagrange interpolation.›
theory Newton_Interpolation
imports 
  "HOL-Library.Monad_Syntax"
  Ring_Hom_Poly
  Divmod_Int
  Is_Rat_To_Rat
begin

text ‹For the Newton interpolation, we start with an efficient implementation (which in prior examples
  we used as an uncertified oracle). Later on, a more abstract definition of the algorithm
  is described for which soundness is proven, and which is provably equivalent to the efficient
  implementation.

  The implementation is based on divided differences and the Horner schema.›

fun horner_composition :: "'a :: comm_ring_1 list ⇒ 'a list ⇒ 'a poly" where
  "horner_composition [cn] xis = [:cn:]"
| "horner_composition (ci # cs) (xi # xis) = horner_composition cs xis * [:- xi, 1:] + [:ci:]"
| "horner_composition _ _ = 0"

lemma (in map_poly_comm_ring_hom) horner_composition_hom:
  "horner_composition (map hom cs) (map hom xs) = map_poly hom (horner_composition cs xs)"
  by (induct cs xs rule: horner_composition.induct, auto simp: hom_distribs)

lemma horner_coeffs_ints: assumes len: "length cs ≤ Suc (length ys)"
  shows "(set (coeffs (horner_composition cs (map rat_of_int ys))) ⊆ ℤ) = (set cs ⊆ ℤ)"
proof -
  let ?ir = "int_of_rat"
  let ?ri = "rat_of_int"
  let ?mir = "map ?ir"
  let ?mri = "map ?ri"
  show ?thesis
  proof
    define ics where "ics = map ?ir cs"
    assume "set cs ⊆ ℤ"
    hence ics: "cs = ?mri ics" unfolding ics_def map_map o_def
      by (simp add: map_idI subset_code(1))      
    show "set (coeffs (horner_composition cs (?mri ys))) ⊆ ℤ"
      unfolding ics of_int_poly_hom.horner_composition_hom by auto
  next
    assume "set (coeffs (horner_composition cs (?mri ys))) ⊆ ℤ"
    thus "set cs ⊆ ℤ" using len
    proof (induct cs arbitrary: ys)
      case (Cons c cs xs)
      show ?case
      proof (cases "cs = [] ∨ xs = []")
        case True
        with Cons show ?thesis by (cases "c = 0"; cases cs, auto)
      next
        case False
        then obtain d ds and y ys where cs: "cs = d # ds" and xs: "xs = y # ys" 
          by (cases cs, auto, cases xs, auto)
        let ?q = "horner_composition cs (?mri ys)"
        define q where "q = ?q"
        define p where "p = q * [:- ?ri y, 1:] + [:c:]"
        have id: "horner_composition (c # cs) (?mri xs) = p" 
          unfolding cs xs q_def p_def by simp
        have coeff: "coeff p i ∈ ℤ" for i
        proof (cases "coeff p i ∈ set (coeffs p)")
          case True
          with Cons(2)[unfolded id] show ?thesis by blast
        next
          case False
          hence "coeff p i = 0" using range_coeff[of p] by blast
          thus ?thesis by simp
        qed
        {
          fix i
          let ?f = "λ j. coeff [:- ?ri y, 1:] j * coeff q (Suc i - j)"
          have "coeff p (Suc i) = coeff ([: -?ri y, 1 :] * q) (Suc i)" unfolding p_def by simp
          also have "… = (∑j≤Suc i. ?f j)" unfolding coeff_mult by simp
          also have "… = ?f 0 + ?f 1 + (∑j∈{..Suc i} - {0} - {Suc 0}. ?f j)"
            by (subst sum.remove[of _ 0], force+, subst sum.remove[of _ 1], force+)
          also have "(∑j∈{..Suc i} - {0} - {Suc 0}. ?f j) = 0"
          proof (rule sum.neutral, auto, goal_cases)
            case (1 x)
            thus ?case by (cases x, auto, cases "x - 1", auto)
          qed
          also have "?f 0 = - ?ri y * coeff q (Suc i)" by simp
          also have "?f 1 = coeff q i" by simp
          finally have int: "coeff q i - ?ri y * coeff q (Suc i) ∈ ℤ" using coeff[of "Suc i"] by auto
          assume "coeff q (Suc i) ∈ ℤ"
          hence "?ri y * coeff q (Suc i) ∈ ℤ" by simp
          hence "coeff q i ∈ ℤ" using int Ints_diff Ints_minus by force
        } note coeff_q = this
        {
          fix i
          assume "i ≤ degree q"
          hence "coeff q (degree q - i) ∈ ℤ"
          proof (induct i)
            case 0
            from coeff_q[of "degree q"] show ?case
              by (metis Ints_0 Suc_n_not_le_n diff_zero le_degree)
          next
            case (Suc i)
            with coeff_q[of i] show ?case
              by (metis Suc_diff_Suc Suc_leD Suc_n_not_le_n coeff_q le_less)
          qed
        } note coeff_q = this
        {
          fix i
          have "coeff q i ∈ ℤ"
          proof (cases "i ≤ degree q")
            case True
            with coeff_q[of "degree q - i"] show ?thesis by auto
          next
            case False
            hence "coeff q i = 0" using le_degree by blast
            thus ?thesis by simp
          qed
        } note coeff_q = this
        hence "set (coeffs q) ⊆ ℤ" by (auto simp: coeffs_def)
        
        from Cons(1)[OF this[unfolded q_def]] Cons(3) xs have IH: "set cs ⊆ ℤ" by auto
        define r where "r = coeff q 0 * (- ?ri y)"
        have r: "r ∈ ℤ" using coeff_q[of 0] unfolding r_def by auto
        have "coeff p 0 ∈ ℤ" by fact
        also have "coeff p 0 = r + c" unfolding p_def r_def by simp
        finally have c: "c ∈ ℤ" using r using Ints_diff by force
        with IH show ?thesis by auto
      qed
    qed simp
  qed
qed

context
fixes 
  ty :: "'a :: field itself"
  and xs :: "'a list"
  and fs :: "'a list"
begin


fun divided_differences_impl :: "'a list ⇒ 'a ⇒ 'a ⇒ 'a list ⇒ 'a list" where
  "divided_differences_impl (xi_j1 # x_j1s) fj xj (xi # xis) = (let 
    x_js = divided_differences_impl x_j1s fj xj xis;
    new = (hd x_js - xi_j1) / (xj - xi)
    in new # x_js)"
| "divided_differences_impl [] fj xj xis = [fj]"

fun newton_coefficients_main :: "'a list ⇒ 'a list ⇒ 'a list list" where
  "newton_coefficients_main [fj] xjs = [[fj]]"
| "newton_coefficients_main (fj # fjs) (xj # xjs) = (
    let rec = newton_coefficients_main fjs xjs; row = hd rec;
      new_row = divided_differences_impl row fj xj xs
    in new_row # rec)"
| "newton_coefficients_main _ _ = []"

definition newton_coefficients :: "'a list" where
  "newton_coefficients = map hd (newton_coefficients_main (rev fs) (rev xs))"

definition newton_poly_impl :: "'a poly" where
  "newton_poly_impl = horner_composition (rev newton_coefficients) xs"

qualified definition "x i = xs ! i"
qualified definition "f i = fs ! i"

private definition "xd i j = x i - x j"

lemma [simp]: "xd i i = 0" "xd i j + xd j k = xd i k" "xd i j + xd k i = xd k j" 
  unfolding xd_def by simp_all

(* divided differences [xi,..,xj]f *)
private function xij_f :: "nat ⇒ nat ⇒ 'a" where
  "xij_f i j = (if i < j then (xij_f (i + 1) j - xij_f i (j - 1)) / xd j i else f i)"
  by pat_completeness auto

termination by (relation "measure (λ (i,j). j - i)", auto)

private definition c :: "nat ⇒ 'a" where
  "c i = xij_f 0 i"

private definition "X j = [: - x j, 1:]"

private function b :: "nat ⇒ nat ⇒ 'a poly" where
  "b i n = (if i ≥ n then [:c n:] else b (Suc i) n * X i + [:c i:])"
  by pat_completeness auto

termination by (relation "measure (λ (i,n). Suc n - i)", auto)

declare b.simps[simp del]

definition newton_poly :: "nat ⇒ 'a poly" where
  "newton_poly n = b 0 n"

private definition "Xij i j = prod_list (map X [i ..< j])"

private definition "N i = Xij 0 i"

lemma Xii_1[simp]: "Xij i i = 1" unfolding Xij_def by simp
lemma smult_1[simp]: "smult d 1 = [:d:]"
  by (fact smult_one)

private lemma newton_poly_sum: 
  "newton_poly n = sum_list (map (λ i. smult (c i) (N i)) [0 ..< Suc n])"
  unfolding newton_poly_def N_def
proof -
  {
    fix j
    assume "j ≤ n"
    hence "b j n = (∑i←[j..<Suc n]. smult (c i) (Xij j i))"
    proof (induct j n rule: b.induct)
      case (1 j n)
      show ?case
      proof (cases "j ≥ n")
        case True
        with 1(2) have j: "j = n" by auto
        hence "b j n = [:c n:]" unfolding b.simps[of j n] by simp
        thus ?thesis unfolding j by simp
      next
        case False
        hence b: "b j n = b (Suc j) n * X j + [: c j:]" unfolding b.simps[of j n] by simp
        define nn where "nn = Suc n"
        from 1(2) have id: "[j..< nn] = j # [Suc j ..< nn]" unfolding nn_def by (simp add: upt_rec)
        from False have "Suc j ≤ n" by auto
        note IH = 1(1)[OF False this]
        have id2: "(∑x←[Suc j..< nn]. smult (c x) (Xij (Suc j) x * X j)) =
          (∑i←[Suc j..< nn]. smult (c i) (Xij j i))"
        proof (rule arg_cong[of _ _ sum_list], rule map_ext, intro impI, goal_cases)
          case (1 i)
          hence "Xij (Suc j) i * X j = Xij j i" by (simp add: Xij_def upt_conv_Cons)
          thus ?case by simp
        qed
        show ?thesis unfolding b IH sum_list_mult_const[symmetric]
          unfolding nn_def[symmetric] id
          by (simp add: id2)
      qed
    qed
  }
  from this[of 0] show "b 0 n = (∑i←[0..<Suc n]. smult (c i) (Xij 0 i))" by simp
qed

private lemma poly_newton_poly: "poly (newton_poly n) y = sum_list (map (λ i. c i * poly (N i) y) [0 ..< Suc n])"
  unfolding newton_poly_sum poly_sum_list map_map o_def by simp

private definition "pprod k i j = (∏l←[i..<j]. xd k l)"

private lemma poly_N_xi: "poly (N i) (x j) = pprod j 0 i"
proof -
  have "poly (N i) (x j) = (∏l←[0..<i]. xd j l)"
    unfolding N_def Xij_def poly_prod_list X_def[abs_def] map_map o_def xd_def by simp
  also have "… = pprod j 0 i" unfolding pprod_def ..
  finally show ?thesis .
qed

private lemma poly_N_xi_cond: "poly (N i) (x j) = (if j < i then 0 else pprod j 0 i)"
proof -
  show ?thesis
  proof (cases "j < i")
    case False
    thus ?thesis using poly_N_xi by simp
  next
    case True
    hence "j ∈ set [0 ..< i]" by auto
    from split_list[OF this] obtain bef aft where id2: "[0 ..< i] = bef @ j # aft" by auto
    have "(∏k←[0..<i]. xd j k) = 0" unfolding id2 by auto
    with True show ?thesis unfolding poly_N_xi pprod_def by auto
  qed
qed

private lemma poly_newton_poly_xj: assumes "j ≤ n"
  shows "poly (newton_poly n) (x j) = sum_list (map (λ i. c i * poly (N i) (x j)) [0 ..< Suc j])"
proof -
  from assms have id: "[0 ..< Suc n] = [0 ..< Suc j] @ [Suc j ..< Suc n]" 
    by (metis Suc_le_mono le_Suc_ex less_eq_nat.simps(1) upt_add_eq_append)
  have id2: "(∑i←[Suc j..< Suc n]. c i * poly (N i) (x j)) = 0"
    by (rule sum_list_neutral, unfold poly_N_xi_cond, auto)
  show ?thesis unfolding poly_newton_poly id map_append sum_list_append id2 by simp
qed

declare xij_f.simps[simp del]

context
  fixes n
  assumes dist: "⋀ i j. i < j ⟹ j ≤ n ⟹ x i ≠ x j"
begin
private lemma xd_diff: "i < j ⟹ j ≤ n ⟹ xd i j ≠ 0"
   "i < j ⟹ j ≤ n ⟹ xd j i ≠ 0" using dist[of i j] dist[of j i] unfolding xd_def by auto

text ‹This is the key technical lemma for soundness of Newton interpolation.›

private lemma divided_differences_main: assumes "k ≤ n" "i < k"
  shows "sum_list (map (λ j. xij_f i (i + j) * pprod k i (i + j)) [0..<Suc k - i]) = 
  sum_list (map (λ j. xij_f (Suc i) (Suc i + j) * pprod k (Suc i) (Suc i + j)) [0..<Suc k - Suc i])"
proof -
  let ?exp = "λ i j. xij_f i (i + j) * pprod k i (i + j)"
  define ei where "ei = ?exp i"
  define esi where "esi = ?exp (Suc i)"  
  let ?ki = "k - i"
  let ?sumi = "λ xs. sum_list (map ei xs)"
  let ?sumsi = "λ xs. sum_list (map esi xs)"
  let ?mid = "λ j. xij_f i (k - j) * pprod k (Suc i) (k - j) * xd (k - j) i"
  let ?sum = "λ j. ?sumi [0 ..< ?ki - j] + ?sumsi [?ki - j ..< ?ki] + ?mid j"
  define fin where "fin = ?ki - 1"
  have fin: "fin < ?ki" unfolding fin_def using assms by auto
  have id: "[ 0 ..< Suc k - i] = [0 ..< ?ki] @ [?ki]" and 
    id2: "[i..<k] = i # [Suc i ..< k]" and
    id3: "k - (i + (k - Suc i)) = 1" "k - (?ki - 1) = Suc i" using assms 
    by (auto simp: Suc_diff_le upt_conv_Cons)
  have neq: "xd (Suc i) i ≠ 0" using xd_diff[of i "Suc i"] assms by auto
  have "sum_list (map (λ j. xij_f i (i + j) * pprod k i (i + j)) [0..<Suc k - i])
    = ?sumi [0 ..< Suc k - i]" unfolding ei_def by simp
  also have "… = ?sumi [0 ..< ?ki] + ?sumsi [?ki ..< ?ki] + ei ?ki" 
    unfolding id by simp
  also have "… = ?sum 0"
    unfolding ei_def using assms by (simp add: pprod_def id2)
  also have "?sum 0 = ?sum fin" using fin
  proof (induct fin)
    case (Suc fin)
    from Suc(2) assms 
    have fki: "fin < ?ki" and ikf: "i < k - Suc fin" "i < k - fin" and kfn: "k - fin ≤ n" by auto
    from xd_diff[OF ikf(2) kfn] have nz: "xd (k - fin) i ≠ 0" by auto
    note IH = Suc(1)[OF fki]
    have id4: "[0 ..< ?ki - fin] = [0 ..< ?ki - Suc fin] @ [?ki - Suc fin]" 
      "i + (k - i - Suc fin) = k - Suc fin" 
      "Suc (k - Suc fin) = k - fin" using Suc(2) assms  ‹fin < ?ki›
      by (metis Suc_diff_Suc le0 upt_Suc) (insert Suc(2), auto)      
    from Suc(2) assms have id5: "[i..<k - Suc fin] = i # [Suc i ..< k - Suc fin]"
      "[Suc i..<k - fin] = [Suc i..<k - Suc fin] @ [k - Suc fin]"
      by (force simp: upt_rec) (metis Suc_leI id4(3) ikf(1) upt_Suc)
    have "?sum 0 = ?sum fin" by (rule IH)
    also have "… = ?sumi [0 ..< ?ki - Suc fin] + ?sumsi [?ki - fin ..< ?ki] + 
      (ei (?ki - Suc fin) + ?mid fin)"
      unfolding id4 by simp
    also have "?mid fin = (xij_f (Suc i) (k - fin) - xij_f i (k - Suc fin))
      * pprod k (Suc i) (k - fin)" unfolding xij_f.simps[of i "k - fin"]
      using ikf nz by simp
    also have "… = xij_f (Suc i) (k - fin) * pprod k (Suc i) (k - fin) -
      xij_f i (k - Suc fin) * pprod k (Suc i) (k - fin)" by algebra
    also have "xij_f (Suc i) (k - fin) * pprod k (Suc i) (k - fin) = esi (?ki - Suc fin)"
      unfolding esi_def using ikf by (simp add: id4)
    also have "ei (?ki - Suc fin) = xij_f i (k - Suc fin) * pprod k i (k - Suc fin)"       
      unfolding ei_def id4 using ikf by (simp add: ac_simps)
    finally have "?sum 0 = ?sumi [0 ..< ?ki - Suc fin] 
      + (esi (?ki - Suc fin) + ?sumsi [?ki - fin ..< ?ki])
      + (xij_f i (k - Suc fin) * (pprod k i (k - Suc fin) - pprod k (Suc i) (k - fin)))"
      by algebra
    also have "esi (?ki - Suc fin) + ?sumsi [?ki - fin ..< ?ki] 
      = ?sumsi ((?ki - Suc fin) # [?ki - fin ..< ?ki])" by simp
    also have "(?ki - Suc fin) # [?ki - fin ..< ?ki] = [?ki - Suc fin ..< ?ki]"
      using Suc(2) by (simp add: Suc_diff_Suc upt_rec)
    also have "pprod k i (k - Suc fin) - pprod k (Suc i) (k - fin) 
      = (xd k i) * pprod k (Suc i) (k - Suc fin) - (xd k (k - Suc fin)) * pprod k (Suc i) (k - Suc fin)"
      unfolding pprod_def id5 by simp
    also have "… = (xd k i - xd k (k - Suc fin)) * pprod k (Suc i) (k - Suc fin)"
      by algebra
    also have "… = (xd (k - Suc fin) i) * pprod k (Suc i) (k - Suc fin)" unfolding xd_def by simp
    also have "xij_f i (k - Suc fin) * … = ?mid (Suc fin)" by simp
    finally show ?case by simp
  qed simp
  also have "… = (ei 0 + ?mid (k - i - 1)) + ?sumsi [1 ..< k - i]" 
    unfolding fin_def by (simp add: id3)
  also have "ei 0 + ?mid (k - i - 1) = esi 0" unfolding id3
    unfolding ei_def esi_def xij_f.simps[of i i]  using neq assms
    by (simp add: field_simps xij_f.simps pprod_def)
  also have "esi 0 + ?sumsi [1 ..< k - i] = ?sumsi (0 # [1 ..< k - i])" by simp
  also have "0 # [1 ..< k - i] = [0 ..< Suc k - Suc i]" 
    using assms by (simp add: upt_rec)
  also have "?sumsi … = sum_list (map (λ j. xij_f (Suc i) (Suc i + j) * 
    pprod k (Suc i) (Suc i + j)) [0..<Suc k - Suc i])"
    unfolding esi_def using assms by simp
  finally show ?thesis .
qed

private lemma divided_differences: assumes kn: "k ≤ n" and ik: "i ≤ k"
  shows "sum_list (map (λ j. xij_f i (i + j) * pprod k i (i + j)) [0..<Suc k - i]) = f k"
proof -
  {
    fix ii
    assume "i + ii ≤ k"
    hence "sum_list (map (λ j. xij_f i (i + j) * pprod k i (i + j)) [0..<Suc k - i])
      = sum_list (map (λ j. xij_f (i + ii) (i + ii + j) * pprod k (i + ii) (i + ii + j)) [0..<Suc k - (i + ii)])"
    proof (induct ii)
      case (Suc ii)
      hence le1: "i + ii ≤ k" and le2: "i + ii < k" by simp_all
      show ?case unfolding Suc(1)[OF le1] unfolding divided_differences_main[OF kn le2]
        using Suc(2) by simp
    qed simp
  } note main = this
  have ik: "i + (k - i) ≤ k" and id: "i + (k - i) = k" using ik by simp_all
  show ?thesis unfolding main[OF ik] unfolding id
    by (simp add: xij_f.simps pprod_def)
qed      

lemma newton_poly_sound: assumes "k ≤ n"
  shows "poly (newton_poly n) (x k) = f k"
proof -
  have "poly (newton_poly n) (x k) = 
    sum_list (map (λ j. xij_f 0 (0 + j) * pprod k 0 (0 + j)) [0..<Suc k - 0])"
    unfolding poly_newton_poly_xj[OF assms] c_def poly_N_xi by simp
  also have "… = f k"
    by (rule divided_differences[OF assms], simp)
  finally show ?thesis by simp
qed
end

lemma newton_poly_degree: "degree (newton_poly n) ≤ n" 
proof -
  {
    fix i
    have "i ≤ n ⟹ degree (b i n) ≤ n - i"
    proof (induct i n rule: b.induct)
      case (1 i n)
      note b = b.simps[of i n]
      show ?case 
      proof (cases "n ≤ i")
        case True
        thus ?thesis unfolding b by auto
      next
        case False
        have "degree (b i n) = degree (b (Suc i) n * X i + [:c i:])" using False unfolding b by simp
        also have "… ≤ max (degree (b (Suc i) n * X i)) (degree [:c i:])"
          by (rule degree_add_le_max)
        also have "… = degree (b (Suc i) n * X i)" by simp
        also have "… ≤ degree (b (Suc i) n) + degree (X i)"
          by (rule degree_mult_le)
        also have "… ≤ n - Suc i + degree (X i)" 
          using 1(1)[OF False] 1(2) False add_le_mono1 not_less_eq_eq by blast
        also have "… = n - Suc i + 1" unfolding X_def by simp
        also have "… = n - i" using 1(2) False by auto
        finally show ?thesis .
      qed
    qed
  }
  from this[of 0] show ?thesis unfolding newton_poly_def by simp
qed

context
  fixes n
  assumes xs: "length xs = n"
    and fs: "length fs = n"
begin
lemma newton_coefficients_main: 
  "k < n ⟹ newton_coefficients_main (rev (map f [0..<Suc k])) (rev (map x [0..<Suc k]))
    = rev (map (λ i. map (λ j. xij_f j i) [0..<Suc i]) [0..<Suc k])"
proof (induct k)
  case 0
  show ?case
    by (simp add: xij_f.simps)
next
  case (Suc k)
  hence "k < n" by auto
  note IH = Suc(1)[OF this]
  have id: "⋀ f. rev (map f [0..<Suc (Suc k)]) = f (Suc k) # f k # rev (map f [0..< k])" 
    and id2: "⋀ f. f k # rev (map f [0..<k]) = rev (map f [0..< Suc k])" by simp_all
  show ?case unfolding id newton_coefficients_main.simps Let_def
    unfolding id2 IH
    unfolding list.simps id2[symmetric]
  proof (rule conjI, goal_cases)
    case 1
    have xs: "xs = map x [0 ..< n]" using xs unfolding x_def[abs_def] 
      by (intro nth_equalityI, auto)
    define nn where "nn = (0 :: nat)"
    define m where "m = Suc k - nn"
    have prems: "m = Suc k - nn" "nn < Suc (Suc k)" unfolding m_def nn_def by auto
    have "?case = (divided_differences_impl (map ((λj. xij_f j k)) [nn..< Suc k]) (f (Suc k)) (x (Suc k)) (map x [nn ..< n]) =
      map ((λj. xij_f j (Suc k))) [nn..<Suc (Suc k)])"
      unfolding nn_def xs[symmetric] by simp
    also have "…" using prems
    proof (induct m arbitrary: nn)
      case 0
      hence nn: "nn = Suc k" by auto
      show ?case unfolding nn by (simp add: xij_f.simps)
    next
      case (Suc m)
      with ‹Suc k < n› have "nn < n" and le: "nn < Suc k" by auto
      with Suc(2-) have id: 
        "[nn..<Suc k] = nn # [Suc nn..< Suc k]"
        "[nn..<n] = nn # [Suc nn..< n]"
      and id2: "[nn..<Suc (Suc k)] = nn # [Suc nn..<Suc (Suc k)]"
        "[Suc nn..<Suc (Suc k)] = Suc nn # [Suc (Suc nn)..<Suc (Suc k)]"
        by (auto simp: upt_rec)
      from Suc(2-) have "m = Suc k - Suc nn" "Suc nn < Suc (Suc k)" by auto
      note IH = Suc(1)[OF this]
      show ?case unfolding id list.simps divided_differences_impl.simps IH Let_def
        unfolding id2 list.simps 
        using le
        by (simp add: xij_f.simps[of nn "Suc k"] xd_def)
    qed
    finally show ?case by simp
  qed simp
qed

lemma newton_coefficients: "newton_coefficients = rev (map c [0 ..< n])"
proof (cases n)
  case 0
  hence xs: "xs = []" "fs = []" using xs fs by auto
  show ?thesis unfolding newton_coefficients_def 0 
    using newton_coefficients_main.simps
    unfolding xs by simp
next
  case (Suc nn)
  hence sn: "Suc nn = n" and nn: "nn < n" by auto
  from fs have fs: "map f [0..<Suc nn] = fs" unfolding sn
    by (intro nth_equalityI, auto simp: f_def)
  from xs have xs: "map x [0..<Suc nn] = xs" unfolding sn
    by (intro nth_equalityI, auto simp: x_def)
  show ?thesis
    unfolding newton_coefficients_def
      newton_coefficients_main[OF nn, unfolded fs xs]
    unfolding sn rev_map[symmetric] map_map o_def 
    by (rule arg_cong[of _ _ rev], subst upt_rec, intro nth_equalityI, auto simp: c_def)
qed

lemma newton_poly_impl: assumes "n = Suc nn"
  shows "newton_poly_impl = newton_poly nn"
proof -
  define i where "i = (0 :: nat)"
  have xs: "map x [0..<n] = xs" using xs
    by (intro nth_equalityI, auto simp: x_def)
  have "i ≤ nn" unfolding i_def by simp
  hence "horner_composition (map c [i..<Suc nn]) (map x [i..<Suc nn]) = b i nn"
  proof (induct i nn rule: b.induct)
    case (1 i n)
    show ?case
    proof (cases "n ≤ i")
      case True
      with 1(2) have i: "i = n" by simp
      show ?thesis unfolding i b.simps[of n n] by simp
    next
      case False
      hence "Suc i ≤ n" by simp
      note IH = 1(1)[OF False this]
      have bi: "b i n = b (Suc i) n * X i + [:c i:]" using False by (simp add: b.simps)    
      from False have id: "[i ..< Suc n] = i # [Suc i ..< Suc n]" by (simp add: upt_rec)
      from False have id2: "[Suc i ..< Suc n] = Suc i # [Suc (Suc i) ..< Suc n]" by (simp add: upt_rec)
      show ?thesis unfolding id bi list.simps horner_composition.simps id2
        unfolding IH[unfolded id2 list.simps] by (simp add: X_def)
    qed
  qed
  thus ?thesis
    unfolding newton_poly_impl_def newton_coefficients rev_rev_ident newton_poly_def i_def
      assms[symmetric] xs .
qed
end
end

context 
  fixes xs fs :: "int list"
begin

fun divided_differences_impl_int :: "int list ⇒ int ⇒ int ⇒ int list ⇒ int list option" where
  "divided_differences_impl_int (xi_j1 # x_j1s) fj xj (xi # xis) = (
     case divided_differences_impl_int x_j1s fj xj xis of None ⇒ None
   | Some x_js ⇒ let (new,m) = divmod_int (hd x_js - xi_j1) (xj - xi)
     in if m = 0 then Some (new # x_js) else None)"
| "divided_differences_impl_int [] fj xj xis = Some [fj]"
   
fun newton_coefficients_main_int :: "int list ⇒ int list ⇒ int list list option" where
  "newton_coefficients_main_int [fj] xjs = Some [[fj]]"
| "newton_coefficients_main_int (fj # fjs) (xj # xjs) = (do {
    rec ← newton_coefficients_main_int fjs xjs;
    let row = hd rec;
    new_row ← divided_differences_impl_int row fj xj xs;
    Some (new_row # rec)})"
| "newton_coefficients_main_int _ _ = Some []"

definition newton_coefficients_int :: "int list option" where
  "newton_coefficients_int = map_option (map hd) (newton_coefficients_main_int (rev fs) (rev xs))"

lemma divided_differences_impl_int_Some:
  "length gs ≤ length ys 
  ⟹ divided_differences_impl_int gs g x ys = Some res
  ⟹ divided_differences_impl (map rat_of_int gs) (rat_of_int g) (rat_of_int x) (map rat_of_int ys) = map rat_of_int res
    ∧ length res = Suc (length gs)"
proof (induct gs g x ys arbitrary: res rule: divided_differences_impl_int.induct)
  case (1 xi_j1 x_j1s fj xj xi xis)
  note some = 1(3)
  from 1(2) have len: "length x_j1s ≤ length xis" by auto
  from some obtain x_js where rec: "divided_differences_impl_int x_j1s fj xj xis = Some x_js"
    by (auto split: option.splits)
  note IH = 1(1)[OF len rec]
  have id: "hd (map rat_of_int x_js) = rat_of_int (hd x_js)" using IH by (cases x_js, auto)
  from some[simplified, unfolded rec divmod_int_def] have mod: "(hd x_js - xi_j1) mod (xj - xi) = 0"
    and res: "res = (hd x_js - xi_j1) div (xj - xi) # x_js" by (auto split: if_splits)
  have "rat_of_int ((hd x_js - xi_j1) div (xj - xi)) = rat_of_int (hd x_js - xi_j1) / rat_of_int (xj - xi)"
    using mod by force
  hence "(rat_of_int (hd x_js) - rat_of_int xi_j1) / (rat_of_int xj - rat_of_int xi) = 
    rat_of_int ((hd x_js - xi_j1) div (xj - xi))"
    by simp
  thus ?case by (simp add: IH Let_def res id)
next
  case (2 fj xj xis res)
  hence res: "res = [fj]" by simp
  thus ?case by simp
qed simp

lemma div_Ints_mod_0: assumes "rat_of_int a / rat_of_int b ∈ ℤ" "b ≠ 0"
  shows "a mod b = 0"
proof -
  define c where "c = int_of_rat (rat_of_int a / rat_of_int b)"
  have "rat_of_int a / rat_of_int b = rat_of_int c" unfolding c_def using assms(1) by simp
  hence "rat_of_int a = rat_of_int b * rat_of_int c" using assms(2)
    by (metis divide_cancel_right nonzero_mult_div_cancel_left of_int_eq_0_iff)
  hence a: "a = b * c" by (simp add: of_int_hom.injectivity)
  show "a mod b = 0" unfolding a by simp
qed

lemma divided_differences_impl_int_None:
  "length gs ≤ length ys 
  ⟹ divided_differences_impl_int gs g x ys = None
  ⟹ x ∉ set (take (length gs) ys)
  ⟹ hd (divided_differences_impl (map rat_of_int gs) (rat_of_int g) (rat_of_int x) (map rat_of_int ys)) ∉ ℤ"
proof (induct gs g x ys rule: divided_differences_impl_int.induct)
  case (1 xi_j1 x_j1s fj xj xi xis)
  note none = 1(3)
  from 1(2,4) have len: "length x_j1s ≤ length xis" and xj: "xj ∉ set (take (length x_j1s) xis)" and xji: "xj ≠ xi" by auto
  define d where "d = divided_differences_impl (map rat_of_int x_j1s) (rat_of_int fj) (rat_of_int xj) (map rat_of_int xis)"
  note IH = 1(1)[OF len _ xj]
  show ?case
  proof (cases "divided_differences_impl_int x_j1s fj xj xis")
    case None
    from IH[OF None] have d: "hd d ∉ ℤ" unfolding d_def by auto
    {
      let ?x = "(hd d - rat_of_int xi_j1) / (rat_of_int xj - rat_of_int xi)"
      assume "?x ∈ ℤ"
      hence "?x * (of_int (xj - xi)) + rat_of_int xi_j1 ∈ ℤ"
        using Ints_mult Ints_add Ints_of_int by blast
      also have "?x * (of_int (xj - xi)) = hd d - rat_of_int xi_j1" using xji by auto
      also have "… + rat_of_int xi_j1 = hd d" by simp
      finally have False using d by simp
    }
    thus ?thesis 
      by (auto simp: Let_def d_def[symmetric])
  next
    case (Some res)
    from divided_differences_impl_int_Some[OF len Some]
    have id: "divided_differences_impl (map rat_of_int x_j1s) (rat_of_int fj) (rat_of_int xj) (map rat_of_int xis) =
      map rat_of_int res" and res: "res ≠ []" by auto
    have hd: "hd (map rat_of_int res) = of_int (hd res)" using res by (cases res, auto)
    define a where "a = (hd res - xi_j1)"
    define b where "b = xj - xi"
    from none[simplified, unfolded Some divmod_int_def] 
    have mod: "a mod b ≠ 0"
      by (auto split: if_splits simp: a_def b_def)
    {
      assume "(rat_of_int (hd res) - rat_of_int xi_j1) / (rat_of_int xj - rat_of_int xi) ∈ ℤ"
      hence "rat_of_int a / rat_of_int b ∈ ℤ" unfolding a_def b_def by simp
      moreover have "b ≠ 0" using xji unfolding b_def by simp
      ultimately have False using mod div_Ints_mod_0 by auto
    }
    thus ?thesis 
      by (auto simp: id Let_def hd)
  qed
qed auto
  
lemma newton_coefficients_main_int_Some:
  "length gs = length ys ⟹ length ys ≤ length xs  
  ⟹ newton_coefficients_main_int gs ys = Some res
  ⟹ newton_coefficients_main (map rat_of_int xs) (map rat_of_int gs) (map rat_of_int ys) = map (map rat_of_int) res
    ∧ (∀ x ∈ set res. x ≠ [] ∧ length x ≤ length ys) ∧ length res = length gs"
proof (induct gs ys arbitrary: res rule: newton_coefficients_main_int.induct)
  case (2 fv v va xj xjs res)
  from 2(2,3) have len: "length (v # va) = length xjs" "length xjs ≤ length xs" by auto
  note some = 2(4)
  let ?n = "newton_coefficients_main_int (v # va) xjs"
  let ?ri = rat_of_int
  let ?mri = "map ?ri"
  from some obtain rec where n: "?n = Some rec"
    by (cases ?n, auto)
  note some = some[simplified, unfolded n]
  let ?d = "divided_differences_impl_int (hd rec) fv xj xs"
  from some obtain dd where d: "?d = Some dd" and res: "res = dd # rec" 
    by (cases ?d, auto)
  note IH = 2(1)[OF len n]
  from IH have lenn: "length (hd rec) ≤ length xjs" by (cases rec, auto)
  with len have "length (hd rec) ≤ length xs" by auto
  note dd = divided_differences_impl_int_Some[OF this d]
  have hd: "hd (map ?mri rec) = ?mri (hd rec)" using IH by (cases rec, auto)
  show ?case unfolding newton_coefficients_main.simps list.simps
    IH[THEN conjunct1, unfolded list.simps] Let_def hd
    dd[THEN conjunct1] res 
  proof (intro conjI)
    show "length (dd # rec) = length (fv # v # va)" using len
      IH[THEN conjunct2] dd[THEN conjunct2] by auto
    show "∀x∈insert dd (set rec). x ≠ [] ∧ length x ≤ length (xj # xjs)"
      using len IH[THEN conjunct2] dd[THEN conjunct2] lenn by auto
  qed auto
qed auto

lemma newton_coefficients_main_int_None: assumes dist: "distinct xs"
  shows "length gs = length ys ⟹ length ys ≤ length xs  
  ⟹ newton_coefficients_main_int gs ys = None
  ⟹ ys = drop (length xs - length ys) (rev xs)
  ⟹ ∃ row ∈ set (newton_coefficients_main (map rat_of_int xs) (map rat_of_int gs) (map rat_of_int ys)). hd row ∉ ℤ"
proof (induct gs ys rule: newton_coefficients_main_int.induct)
  case (2 fv v va xj xjs)
  from 2(2,3) have len: "length (v # va) = length xjs" "length xjs ≤ length xs" by auto
  from arg_cong[OF 2(5), of tl] 2(3)
  have xjs: "xjs = drop (length xs - length xjs) (rev xs)"
    by (metis 2(5) butlast_snoc butlast_take length_drop rev.simps(2) rev_drop rev_rev_ident rev_take)
  note none = 2(4)
  let ?n = "newton_coefficients_main_int (v # va) xjs"
  let ?n' = "newton_coefficients_main (map rat_of_int xs) (map rat_of_int (v # va)) (map rat_of_int xjs)"
  let ?ri = rat_of_int
  let ?mri = "map ?ri"
  show ?case
  proof (cases ?n)
    case None
    from 2(1)[OF len None xjs] obtain row where 
      row: "row∈set ?n'" and "hd row ∉ ℤ" by auto
    thus ?thesis by (intro bexI[of _ row], auto simp: Let_def)
  next
    case (Some rec)
    note some = newton_coefficients_main_int_Some[OF len this]
    hence len': "length (hd rec) ≤ length xjs" by (cases rec, auto)
    hence lenn: "length (hd rec) ≤ length xs" using len by auto
    have hd: "hd (map ?mri rec) = ?mri (hd rec)" using some by (cases rec, auto)
    let ?d = "divided_differences_impl_int (hd rec) fv xj xs"
    from none[simplified, unfolded Some]
    have none: "?d = None" by (cases ?d, auto)
    have "xj ∉ set (take (length (hd rec)) xs)"
    proof
      assume "xj ∈ set (take (length (hd rec)) xs)"
      then obtain i where "i < length (hd rec)" and xj: "xj = xs ! i" 
        unfolding in_set_conv_nth by auto
      with len' have i: "i < length xjs" by simp
      have "Suc (length xjs) ≤ length xs" using 2(3) by auto
      with i have i0: "i ≠ 0"
        by (metis 2(5) Suc_diff_Suc Suc_le_lessD diff_less dist distinct_conv_nth 
          hd_drop_conv_nth length_Cons length_drop length_greater_0_conv length_rev less_le_trans 
          list.sel(1) list.simps(3) nat_neq_iff rev_nth xj xjs)      
      have "xj ∈ set xjs" 
        by (subst xjs, unfold xj in_set_conv_nth, rule exI[of _ "length xjs - Suc i"], insert i 2(3) i0,
           auto simp: rev_nth)
      hence ndist: "¬ distinct (xj # xjs)" by auto
      from dist have "distinct (rev xs)" by simp
      from distinct_drop[OF this] have "distinct (xj # xjs)" using 2(5) by metis
      with ndist
      show False ..
    qed
    note dd = divided_differences_impl_int_None[OF lenn none this]
    show ?thesis
      by (rule bexI, rule dd, insert some hd, auto)
  qed
qed auto


lemma newton_coefficients_int: assumes dist: "distinct xs"
  and len: "length xs = length fs"
  shows "newton_coefficients_int = (let cs = newton_coefficients (map rat_of_int xs) (map of_int fs)
    in if set cs ⊆ ℤ then Some (map int_of_rat cs) else None)" 
proof -
  from len have len: "length (rev fs) = length (rev xs)" "length (rev xs) ≤ length xs" by auto
  show ?thesis
  proof (cases "newton_coefficients_main_int (rev fs) (rev xs)")
    case (Some res)
    have rev: "⋀ xs. map rat_of_int (rev xs) = rev (map of_int xs)" unfolding rev_map ..
    note n = newton_coefficients_main_int_Some[OF len Some, unfolded rev]
    {
      fix row
      assume "row ∈ set res"
      with n have "row ≠ []" by auto
      hence id: "hd (map rat_of_int row) = rat_of_int (hd row)" by (cases row, auto)
      also have "… ∈ ℤ" by auto
      finally have int: "hd (map rat_of_int row) ∈ ℤ" by auto
      have "hd row = int_of_rat (hd (map rat_of_int row))" unfolding id by simp
      note this int
    } 
    thus ?thesis unfolding newton_coefficients_int_def Some newton_coefficients_def n[THEN conjunct1] Let_def option.simps
      by (auto simp: o_def)
  next
    case None
    have "rev xs = drop (length xs - length (rev xs)) (rev xs)" by simp
    from newton_coefficients_main_int_None[OF dist len None this]
    show ?thesis unfolding newton_coefficients_int_def newton_coefficients_def None by (auto simp: Let_def rev_map)
  qed
qed

definition newton_poly_impl_int :: "int poly option" where
  "newton_poly_impl_int ≡ case newton_coefficients_int of None ⇒ None 
     | Some nc ⇒ Some (horner_composition (rev nc) xs)"

lemma newton_poly_impl_int: assumes len: "length xs = length fs" 
  and dist: "distinct xs"
  shows "newton_poly_impl_int = (let p = newton_poly_impl (map rat_of_int xs) (map of_int fs)
    in if set (coeffs p) ⊆ ℤ then Some (map_poly int_of_rat p) else None)"
proof -
  let ?ir = "int_of_rat"
  let ?ri = "rat_of_int"
  let ?mir = "map ?ir"
  let ?mri = "map ?ri"
  let ?nc = "newton_coefficients (?mri xs) (?mri fs)"
  have id: "newton_poly_impl_int = (if set ?nc ⊆ ℤ
    then Some (horner_composition (rev (?mir ?nc)) xs) else None)" 
    unfolding newton_poly_impl_int_def newton_coefficients_int[OF dist len] Let_def by simp
  have len: "length (rev ?nc) ≤ Suc (length xs)" 
    unfolding length_rev
    by (subst newton_coefficients[OF refl], insert len, auto) 
  show ?thesis unfolding id
    unfolding newton_poly_impl_def
    unfolding Let_def set_rev rev_map horner_coeffs_ints[OF len]
  proof (rule if_cong[OF refl _ refl], rule arg_cong[of _ _ Some])
    define cs where "cs = rev ?nc"
    define ics where "ics = map ?ir cs"
    assume "set ?nc ⊆ ℤ"
    hence "set cs ⊆ ℤ" unfolding cs_def by auto
    hence ics: "cs = ?mri ics" unfolding ics_def map_map o_def
      by (simp add: map_idI subset_code(1))      
    have id: "horner_composition (rev ?nc) (?mri xs) = map_poly ?ri (horner_composition ics xs)"
      unfolding cs_def[symmetric] ics
      by (rule of_int_poly_hom.horner_composition_hom)
    show "horner_composition (?mir (rev ?nc)) xs
      = map_poly ?ir (horner_composition (rev ?nc) (?mri xs))"
      unfolding id unfolding cs_def[symmetric] ics_def[symmetric]
      by (subst map_poly_map_poly, auto simp: o_def map_poly_idI)
  qed
qed
end

definition newton_interpolation_poly :: "('a :: field × 'a)list ⇒ 'a poly" where
  "newton_interpolation_poly x_fs = (let 
    xs = map fst x_fs; fs = map snd x_fs in
    newton_poly_impl xs fs)"

definition newton_interpolation_poly_int :: "(int × int)list ⇒ int poly option" where
  "newton_interpolation_poly_int x_fs = (let 
    xs = map fst x_fs; fs = map snd x_fs in
    newton_poly_impl_int xs fs)"   

lemma newton_interpolation_poly: assumes dist: "distinct (map fst xs_ys)"
  and p: "p = newton_interpolation_poly xs_ys"
  and xy: "(x,y) ∈ set xs_ys"
  shows "poly p x = y"
proof (cases "length xs_ys")
  case 0
  thus ?thesis using xy by (cases xs_ys, auto)
next
  case (Suc nn)
  let ?xs = "map fst xs_ys" let ?fs = "map snd xs_ys" let ?n = "Suc nn"
  from xy[unfolded set_conv_nth] obtain i where xy: "i ≤ nn" "x = ?xs ! i" "y = ?fs ! i"
    using Suc
    by (metis (no_types, lifting) fst_conv in_set_conv_nth less_Suc_eq_le nth_map snd_conv xy)
  have id: "newton_interpolation_poly xs_ys = newton_poly ?xs ?fs nn"
    unfolding newton_interpolation_poly_def Let_def
    by (rule newton_poly_impl[OF _ _ Suc], auto)
  show ?thesis
    unfolding p id
  proof (rule newton_poly_sound[of nn ?xs _ ?fs, unfolded 
      Newton_Interpolation.x_def Newton_Interpolation.f_def, OF _ xy(1), folded xy(2-)])
    fix i j
    show "i < j ⟹ j ≤ nn ⟹ ?xs ! i ≠ ?xs ! j" using dist Suc nth_eq_iff_index_eq by fastforce
  qed
qed

lemma degree_newton_interpolation_poly:  
  shows "degree (newton_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 newton_interpolation_poly_def Let_def list.simps newton_poly_impl_def
    Newton_Interpolation.newton_coefficients_def
    by simp
next
  case (Suc nn)
  let ?xs = "map fst xs_ys" let ?fs = "map snd xs_ys" let ?n = "Suc nn"
  have id: "newton_interpolation_poly xs_ys = newton_poly ?xs ?fs nn"
    unfolding newton_interpolation_poly_def Let_def
    by (rule newton_poly_impl[OF _ _ Suc], auto)
  show ?thesis unfolding id using newton_poly_degree[of ?xs ?fs nn] Suc by simp
qed

text ‹For @{const newton_interpolation_poly_int} at this point we just prove that it 
  is equivalent to perfom an interpolation on the rational numbers, and then check
  whether all resulting coefficients are integers. That this corresponds to a 
  sound and complete interpolation algorithm on the integers is proven in the theory
  Polynomial-Interpolation, cf.\ lemmas newton-interpolation-poly-int-Some/None.›

lemma newton_interpolation_poly_int: assumes dist: "distinct (map fst xs_ys)"
  shows "newton_interpolation_poly_int xs_ys = (let 
     rxs_ys = map (λ (x,y). (rat_of_int x, rat_of_int y)) xs_ys;
     rp = newton_interpolation_poly rxs_ys
     in if (∀ x ∈ set (coeffs rp). is_int_rat x) then
       Some (map_poly int_of_rat rp) else None)"
proof -
  have id1: "map fst (map (λ(x, y). (rat_of_int x, rat_of_int y)) xs_ys) = map rat_of_int (map fst xs_ys)"
    by (induct xs_ys, auto)
  have id2: "map snd (map (λ(x, y). (rat_of_int x, rat_of_int y)) xs_ys) = map rat_of_int (map snd xs_ys)"
    by (induct xs_ys, auto)
  have id3: "length (map fst xs_ys) = length (map snd xs_ys)" by auto
  show ?thesis
    unfolding newton_interpolation_poly_def  newton_interpolation_poly_int_def Let_def newton_poly_impl_int[OF id3 dist]
    unfolding id1 id2
    by (rule sym, rule if_cong, auto simp: is_int_rat[abs_def])
qed


hide_const 
  Newton_Interpolation.x
  Newton_Interpolation.f
end