Theory Rational_Function_Interpolation

section ‹Rational Function Interpolation›

theory Rational_Function_Interpolation
  imports
    "Poly_Lemmas"
    "Gauss_Jordan.System_Of_Equations"
    "Polynomial_Interpolation.Missing_Polynomial"
begin

subsection ‹Definitions›

text ‹General condition for rational functions interpolation›

definition interpolated_rational_function where
  "interpolated_rational_function pA pB E fA fB dA dB ≡
    (∀ e ∈ E. fA e * poly pB e = fB e * poly pA e) ∧
    degree pA ≤ (dA::real) ∧ degree pB ≤ (dB::real) ∧
    pA ≠ 0 ∧ pB ≠ 0"

text ‹Interpolation condition with given exact degrees›

definition monic_interpolated_rational_function where
 "monic_interpolated_rational_function pA pB E fA fB dA dB  ≡
    (∀ e ∈ E. fA e * poly pB e = fB e * poly pA e) ∧
    degree pA = ⌊dA::real⌋ ∧ degree pB = ⌊dB::real⌋ ∧
    monic pA ∧ monic pB"

lemma monic0: "¬ monic (0::'a::zero_neq_one poly)"
  by simp

lemma monic_interpolated_rational_function_interpolated_rational_function:
  "monic_interpolated_rational_function pA pB E fA fB dA dB
     ⟹ interpolated_rational_function pA pB E fA fB dA dB ∨ ¬(pA ≠ 0 ∧ pB ≠ 0)"
  unfolding monic_interpolated_rational_function_def interpolated_rational_function_def
  by linarith

definition rfi_coefficient_matrix :: "'a::field list ⇒ ('a ⇒ 'a) ⇒ nat ⇒ nat
    ⇒ nat ⇒ nat ⇒ 'a" where
  "rfi_coefficient_matrix E f dA dB i j = (
    if j < dA then
      (E ! i) ^ j
    else if j < dA + dB then
      - f (E ! i) * (E ! i) ^ (j-dA)
    else 0
  )"

definition rfi_constant_vector :: "'a::field list ⇒ ('a ⇒ 'a) ⇒ nat ⇒ nat ⇒ (nat ⇒ 'a)" where
   "rfi_constant_vector E f dA dB = (λ i. f (E ! i) * (E ! i) ^ dB - (E ! i) ^ dA)"

definition rational_function_interpolation :: "'a::field list ⇒ ('a ⇒ 'a) ⇒ nat ⇒ nat
    ⇒ 'm::mod_type itself ⇒ ('a,'m) vec" where
  "rational_function_interpolation E f dA dB m =
    (let solved = solve
      (χ (i::'m) (j::'m). rfi_coefficient_matrix E f dA dB (to_nat i) (to_nat j))
      (χ (i::'m). rfi_constant_vector E f dA dB (to_nat i))
    in fst (the solved))"

definition solution_to_poly :: "('a::finite_field, 'n::mod_type) vec ⇒
    nat ⇒ nat ⇒ 'a poly × 'a poly" where
  "solution_to_poly S dA dB = (let
    p = Abs_poly (λi. if i < dA then S $ (from_nat i) else 0) + monom 1 dA;
    q = Abs_poly (λi. if i < dB then S $ (from_nat (i+dA)) else 0) + monom 1 dB in
    (p, q))"

definition interpolate_rat_fun where
  "interpolate_rat_fun E f dA dB m =
  solution_to_poly (rational_function_interpolation E f dA dB m) dA dB"


subsection ‹Preliminary Results›

lemma consecutive_sum_combine:
  assumes "m ≥ n"
  shows "(∑i = 0..n. f i) + (∑i = Suc n ..m. f i) = (∑i = 0..m. f i)"
proof -
  from assms have "{0..n} ∪ {Suc n..m} = {0..m}"
    by auto
  moreover have "sum f ({0..n} ∪ {Suc n..m}) =
      sum f ({0..n} - {Suc n..m}) + sum f ({Suc n..m} - {0..n}) + sum f ({0..n} ∩ {Suc n..m})"
    using sum_Un2 finite_atLeastAtMost by fast
  ultimately show ?thesis
    by (simp add: Diff_triv)
qed

lemma poly_altdef_Abs_poly_le:
  fixes x :: "'a::{comm_semiring_0, semiring_1}"
  shows "poly (Abs_poly (λi. if i ≤ n then f i else 0)) x = (∑i = 0..n. f i * x ^ i)"
proof -
  let ?ifA0 = "(λi. if i ≤ n then f i else 0)"
  let ?p = "Abs_poly ?ifA0"

  have co: "coeff ?p = ?ifA0"
    using coeff_Abs_poly_If_le by blast

  then have "∀i>n. coeff ?p i = 0"
    by auto
  then have de: "degree ?p ≤ n"
    using degree_le by blast

  have "∀i>degree ?p. ?ifA0 i = 0"
    using co coeff_eq_0 by fastforce
  then have "∀i>degree ?p. ?ifA0 i * x ^ i = 0"
    by simp
  then have "∀ i ∈ {Suc (degree ?p)..n}. (?ifA0 i * x ^ i) = 0"
    using less_eq_Suc_le by fastforce
  then have db: "(∑i = Suc (degree ?p)..n. ?ifA0 i * x ^ i) = 0"
     by simp

  have "poly ?p x = (∑i≤degree ?p. coeff ?p i * x ^ i)"
    using poly_altdef by auto
  also have "… = (∑i≤degree ?p. ?ifA0 i * x ^ i) "
    using co by simp
  also have "… = (∑i = 0..degree ?p. ?ifA0 i * x ^ i) "
    using atMost_atLeast0 by simp
  also have "… = (∑i = 0..degree ?p. ?ifA0 i * x ^ i) +
      (∑i = Suc (degree ?p)..n. ?ifA0 i * x ^ i)"
    using db by simp
  also have  "… = (∑i = 0..n. ?ifA0 i * x ^ i)"
    using consecutive_sum_combine de by blast
  finally show ?thesis
    by simp
qed

lemma poly_altdef_Abs_poly_l:
  fixes x :: "'a::{comm_semiring_0,semiring_1}"
  shows "poly (Abs_poly (λi. if i < n then f i else 0)) x = (∑i<n. f i * x ^ i)"
proof (cases "n")
  case 0
  have p0: "Abs_poly (λi. 0) = 0"
    using zero_poly_def by fastforce
  show ?thesis
    using 0 by (simp add: p0)
next
  case (Suc m)
  have "poly (Abs_poly (λi. if i ≤ m then f i else 0)) x = (∑i = 0..m. f i * x ^ i)"
    using poly_altdef_Abs_poly_le by blast
  moreover have "poly (Abs_poly (λi. if i ≤ m then f i else 0)) x = poly (Abs_poly (λi. if i < n then f i else 0)) x "
    using Suc using less_Suc_eq_le by auto
  moreover have "(∑i = 0..m. f i * x ^ i) = (∑i<n. f i * x ^ i)"
    using Suc atLeast0AtMost lessThan_Suc_atMost by presburger
  ultimately show ?thesis by argo
qed

lemma degree_Abs_poly_If_l:
  assumes "n ≠ 0"
  shows "degree (Abs_poly (λi. if i < n then f i else 0)) < n"
proof -
  have "coeff (Abs_poly (λi. if i < n then f i else 0)) x = 0" if "x ≥ n" for x
    using coeff_Abs_poly [of n  "(λi. if i < n then f i else 0)"] using that by simp
  then show ?thesis
    using assms degree_lessI by blast
qed

lemma nth_less_length_in_set_eq:
  shows "(∀ i < length E. f (E ! i) = g (E ! i)) ⟷ (∀ e ∈ set E. f e = g e)"
proof standard
  show "∀i < length E. f (E ! i) = g (E ! i) ⟹ ∀e∈set E. f e = g e"
    using in_set_conv_nth by metis
next
  show "∀e∈set E. f e = g e ⟹ ∀i<length E. f (E ! i) = g (E ! i)"
    by simp
qed

lemma nat_leq_real_floor: "real (i::nat) ≤ (d::real) ⟷ real i ≤ ⌊d⌋" (is "?l = ?r")
proof
  assume ?l
  then show ?r
    using floor_mono by fastforce
next
  assume ?r
  then show ?l
    by linarith
qed

lemma mod_type_less_function_eq:
  fixes i :: "'a::mod_type"
  assumes "∀ i < CARD('a) . f i = g i"
  shows "f (to_nat i) = g (to_nat i)"
  using assms by (simp add: to_nat_less_card)


subsection ‹On @{term "solution_to_poly"}›

lemma fst_solution_to_poly_nz:
  "fst (solution_to_poly S dA dB) ≠ 0"
proof
  assume "fst (solution_to_poly S dA dB) = 0"
  hence "coeff (Abs_poly (λi. if i < dA then S $ (from_nat i) else 0) + monom 1 dA) dA = 0"
    unfolding solution_to_poly_def by simp
  hence "coeff (Abs_poly (λi. if i < dA then S $ (from_nat i) else 0)) dA + 1 = 0" by simp
  thus "False" by (subst (asm) coeff_Abs_poly[where n="dA"]) auto
qed

lemma snd_solution_to_poly_nz:
  "snd (solution_to_poly S dA dB) ≠ 0"
proof
  assume "snd (solution_to_poly S dA dB) = 0"
  hence "coeff (Abs_poly (λi. if i < dB then S $ (from_nat (i+dA)) else 0) + monom 1 dB) dB = 0"
    unfolding solution_to_poly_def by simp
  hence "coeff (Abs_poly (λi. if i < dB then S $ (from_nat (i+dA)) else 0)) dB + 1 = 0" by simp
  thus "False" by (subst (asm) coeff_Abs_poly[where n="dB"]) auto
qed

lemma degree_Abs0p1: "degree (Abs_poly (λi. 0) + 1) = 0"
    by (metis add_0 degree_1 zero_poly_def)

lemma degree_solution_to_poly_fst:
  "degree (fst (solution_to_poly S dA dB)) = dA"
proof (cases dA)
  case 0
  then show ?thesis unfolding solution_to_poly_def
    using degree_Abs0p1 by (simp add: one_pCons)
next
  case (Suc nat)
  then have "degree (Abs_poly (λi. if i < dA then S $ from_nat i else 0)) < dA"
    using degree_Abs_poly_If_l by fast
  moreover have "… = degree (monom (1::'a) dA)"
    by (simp add: degree_monom_eq)
  ultimately show ?thesis
    unfolding solution_to_poly_def
    by (simp add: degree_add_eq_right)
qed

lemma degree_solution_to_poly_snd:
  "degree (snd (solution_to_poly S dA dB)) = dB"
proof (cases dB)
  case 0
  then show ?thesis unfolding solution_to_poly_def
    using degree_Abs0p1 by (simp add: one_pCons)
next
  case (Suc nat)
  then have "degree (Abs_poly (λi. if i < dB then S $ from_nat (i + dA) else 0)) < dB"
    using degree_Abs_poly_If_l by fast
  moreover have "… = degree (monom (1::'a) dB)"
    by (simp add: degree_monom_eq)
  ultimately  show ?thesis
    unfolding solution_to_poly_def
    by (simp add: degree_add_eq_right)
qed

lemma monic_solution_to_poly_snd:
  "monic (snd (solution_to_poly S dA dB))"
proof (cases dB)
  case 0
  then show ?thesis unfolding solution_to_poly_def
    by (simp add: coeff_Abs_poly degree_Abs0p1)
next
  case (Suc x)
  have 1: "coeff (Abs_poly (λi. if i < Suc x then S $ from_nat (i + dA) else 0)) (Suc x) = 0"
    by (simp add: coeff_eq_0 degree_Abs_poly_If_l)
  have "degree (Abs_poly (λi. if i < dB then S $ from_nat (i + dA) else 0) + monom 1 dB) = dB"
    using degree_solution_to_poly_snd unfolding solution_to_poly_def by auto
  then show ?thesis
      unfolding solution_to_poly_def using 1 Suc by simp
  qed

lemma monic_solution_to_poly_fst:
  "monic (fst (solution_to_poly S dA dB))"
proof (cases dA)
  case 0
  then show ?thesis
    unfolding solution_to_poly_def by (simp add: coeff_Abs_poly degree_Abs0p1)
next
  case (Suc x)
  have 1: "coeff (Abs_poly (λi. if i < dA then S $ (from_nat i) else 0)) (Suc x) = 0"
    by (simp add: Suc coeff_eq_0 degree_Abs_poly_If_l)
  have "degree (Abs_poly (λi. if i < dA then S $ (from_nat i) else 0) + monom 1 dA) = dA"
    using degree_solution_to_poly_fst unfolding solution_to_poly_def by auto
  then show ?thesis
      unfolding solution_to_poly_def using 1 Suc by simp
qed

subsection "Correctness"

text ‹Needs the assumption that the system is consistent, because a solution exists.›
lemma rational_function_interpolation_correct_poly:
  assumes
    "∀ x ∈ set E. f x = fA x / fB x" "∀ x ∈ set E. fB x ≠ 0"
    "dA + dB ≤ length E"
    "CARD('m::mod_type) = length E"
    "consistent (χ (i::'m) (j::'m). rfi_coefficient_matrix E f dA dB (to_nat i) (to_nat j))
                (χ (i::'m). rfi_constant_vector E f dA dB (to_nat i))"
    "S = rational_function_interpolation E f dA dB TYPE('m)"
    "pA = fst (solution_to_poly S dA dB)"
    "pB = snd (solution_to_poly S dA dB)"
  shows
    "∀ e ∈ set E. fA e * poly pB e = fB e * poly pA e"
proof -

  let ?coeff = "rfi_coefficient_matrix E f dA dB"
  let ?const = "rfi_constant_vector E f dA dB"
  let ?coeff' = "(χ (i::'m) (j::'m). ?coeff (to_nat i) (to_nat j))"
  let ?const' = "(χ (i::'m). ?const (to_nat i))"

  have "is_solution S ?coeff' ?const'"
    by (simp add: assms(5,6) consistent_imp_is_solution_solve rational_function_interpolation_def)
  then have sol: "?coeff' *v S = ?const'"
    by (simp add: is_solution_def)

  have const: "?const i = ?const' $ from_nat i" if "i < length E" for i
    by (simp add: assms(4) that to_nat_from_nat_id)

  have coeff: "?coeff i j = ?coeff' $ from_nat i $ from_nat j"
    if "i < length E" "j < length E" for i j
  proof -
    have "to_nat (from_nat i ::'m) = i"
      using that assms(4)
      by (intro to_nat_from_nat_id) simp
    moreover have "to_nat (from_nat j ::'m) = j"
      using that assms(4,3)
      by (intro to_nat_from_nat_id) simp
    ultimately show ?thesis
      unfolding rfi_coefficient_matrix_def
      by (simp add: Let_def)
  qed

  have x: "(∑j < dA + dB. (?coeff i j) * S $ (from_nat j)) = ?const i"
    (is "?l = ?r") if "i < length E" for i
  proof -
    have "?l = (∑j < length E. ?coeff i j * S $ (from_nat j))"
      using assms(3) by (intro sum.mono_neutral_cong_left) (auto simp add:rfi_coefficient_matrix_def)
    also have "... = (∑j < length E. ?coeff' $ (from_nat i) $ (from_nat j) * S $ (from_nat j))"
      using coeff that by auto
    also have "… = (∑j ∈ {0..< length E}. ?coeff' $ (from_nat i) $ (from_nat j) * S $ (from_nat j))"
      by (intro sum.reindex_bij_betw [symmetric] bij_betwI [where g = id]) auto
    also have "… = (∑j∈(UNIV::'m set). ?coeff' $ (from_nat i) $ j * S $ j)"
      using bij_from_nat [where 'a = 'm] assms(3,4) by (intro sum.reindex_bij_betw) simp
    also have "… = (?coeff' *v S) $ (from_nat i)"
      unfolding matrix_vector_mult_def by simp
    also have "… = ?const' $ (from_nat i)"
      using sol by simp
    finally show "?l = ?r" using const that by simp
  qed

  let ?p_lam = "λi. if i < dA then S $ from_nat i else 0"
  let ?q_lam = "λi. if i < dB then S $ from_nat (i + dA) else 0"
  let ?p' = "Abs_poly ?p_lam + monom 1 dA"
  let ?q' = "Abs_poly ?q_lam + monom 1 dB"
  have pq: "pA = ?p'" "pB = ?q'"
    using assms(7,8) unfolding solution_to_poly_def by auto

  have "(∑j<dA. S $ from_nat j * E ! i ^ j) - f (E ! i) * (∑j<dB. S $ from_nat (j + dA) * E ! i ^ j)
    = f (E ! i) * E ! i ^ dB - E ! i ^ dA" if "i < length E" for i
  proof -
    let ?pq_lam = "(λj. (if j < dA then E ! i ^ j else
        if j < dA + dB then - f (E ! i) * E ! i ^ (j - dA) else 0) * S $ from_nat j)"

    have reindex: "(∑j ∈ {dA..< dA + dB}. - f (E ! i) * E ! i ^ (j - dA) * S $ from_nat j) =
      (∑j ∈ {0..< dB}. - f (E ! i) * E ! i ^ (j) * S $ from_nat (j+dA))"
      by (rule sum.reindex_bij_witness [of _ "λi. i + dA" "λi. i - dA"]) auto

    from x have "f (E ! i) * E ! i ^ dB - E ! i ^ dA = (∑j<dA + dB. ?pq_lam j )"
      unfolding rfi_coefficient_matrix_def rfi_constant_vector_def using that by simp
    also have "… = (∑j ∈ {0..< dA + dB}. ?pq_lam j)"
      using atLeast0LessThan by presburger
    also have "… = (∑j ∈ {0..< dA}. ?pq_lam j) + (∑j ∈ {dA..< dA + dB}. ?pq_lam j) "
      by (subst sum.atLeastLessThan_concat) auto
    also have "… = (∑j ∈ {0..< dA}. E ! i ^ j * S $ from_nat j) +
        (∑j ∈ {dA..< dA + dB}. - f (E ! i) * E ! i ^ (j - dA) * S $ from_nat j)"
      by auto
    also have "… = (∑j ∈ {0..< dA}. E ! i ^ j * S $ from_nat j) +
        (∑j ∈ {0..< dB}. - f (E ! i) * E ! i ^ (j) * S $ from_nat (j+dA))"
      using reindex by simp
    also have "… = (∑j ∈ {0..< dA}. E ! i ^ j * S $ from_nat j) +
        - f (E ! i) * (∑j ∈ {0..< dB}. E ! i ^ (j) * S $ from_nat (j+dA))"
      by (simp add: sum_distrib_left mult.commute mult.left_commute)
    finally have "f (E ! i) * E ! i ^ dB - E ! i ^ dA = …"
      by argo

    moreover have "(∑j ∈ {0..< dA}. E ! i ^ j * S $ from_nat j) =
        (∑j<dA. S $ from_nat j * E ! i ^ j)"
      by (subst atLeast0LessThan) (meson mult.commute)
    moreover have "(∑j ∈ {0..< dB}. E ! i ^ (j) * S $ from_nat (j+dA)) =
        (∑j<dB. S $ from_nat (j + dA) * E ! i ^ j)"
      by (subst atLeast0LessThan) (meson mult.commute)
    ultimately show ?thesis
      by simp
  qed

  then have "∀e ∈ set E. (∑j<dA. S $ from_nat j * e ^ j) - f e * (∑j<dB. S $ from_nat (j + dA) * e ^ j)
      = f e * e ^ dB - e ^ dA"
    by (subst nth_less_length_in_set_eq [symmetric]) auto

  then have "(∑i<dA. S $ from_nat i * e ^ i) - f e * (∑i<dB. S $ from_nat (i + dA) * e ^ i)
      = f e * e ^ dB - e ^ dA" if "e ∈ set E" for e
    using that by blast

  then have "(∑i<dA. S $ from_nat i * e ^ i) + e ^ dA
      = f e * e ^ dB + f e * (∑i<dB. S $ from_nat (i + dA) * e ^ i)" if "e ∈ set E" for e
    using that by (simp add:field_simps)

  then have "f e * ((∑i<dB. S $ from_nat (i + dA) * e ^ i) + e ^ dB) =
      (∑i<dA. S $ from_nat i * e ^ i) + e ^ dA" if "e ∈ set E" for e
    using that by (simp add: ring_class.ring_distribs(1))

  then have "f e * (poly (Abs_poly ?q_lam) e + poly (monom 1 dB) e) =
      poly (Abs_poly ?p_lam) e + poly (monom 1 dA) e" if "e ∈ set E" for e
    unfolding poly_altdef_Abs_poly_l poly_monom using that by auto

  then have "f e * (poly (Abs_poly ?q_lam) e + poly (monom 1 dB) e) =
      poly (Abs_poly ?p_lam) e + poly (monom 1 dA) e" if "e ∈ set E" for e
    using that by simp
  then have "f e * poly (Abs_poly ?q_lam + monom 1 dB) e =
      poly (Abs_poly ?p_lam + monom 1 dA) e" if "e ∈ set E" for e
    by (simp add: that)
  then have "(fA e / fB e) * poly ?q' e =  poly ?p' e" if "e ∈ set E" for e
    using that assms(1) by simp
  then have "fA e * poly ?q' e = fB e * poly ?p' e" if "e ∈ set E" for e
    using that by (simp add: assms(2) nonzero_divide_eq_eq)
  then have "∀e∈set E. fA e * poly (snd (solution_to_poly S dA dB)) e =
      fB e * poly (fst (solution_to_poly S dA dB)) e"
    unfolding solution_to_poly_def by auto
  then show "∀e∈set E. fA e * poly pB e = fB e * poly pA e"
    using assms(8,7) by simp
qed

lemma poly_lead_coeff_extract:
  "poly p x = (∑i<degree p. coeff p i * x ^ i) + lead_coeff p *  x ^ degree p"
    for x :: "'a::{comm_semiring_0,semiring_1}"
  unfolding poly_altdef using lessThan_Suc_atMost sum.lessThan_Suc by auto

lemma dA_dB_helper:
  assumes
    "finite A" "finite B"
    "int dA = ⌊(real (length E) + card A - card B)/2⌋"
    "int dB = ⌊(real (length E) + card B - card A)/2⌋"
    "card (sym_diff A B) ≤ length E"
  shows
    "dA + dB ≤ length E"
    "card (A - B) ≤ dA" "card (B - A) ≤ dB"
    "dB - card (B - A) = dA - card (A - B)"
proof -
  have a: "real dA = of_int ⌊(real (length E) + card A - card B)/2⌋"
    using assms(3) by simp
  have b: "real dB = of_int ⌊(real (length E) + card B - card A)/2⌋"
    using assms(4) by simp

  have "real dA + real dB ≤ (real (length E)+card A-card B)/2 + (real (length E)+card B-card A)/2"
    unfolding a b by (intro add_mono) linarith+
  also have "… = real (length E)" by argo
  finally have "real dA + real dB ≤ length E" by simp
  thus "dA + dB ≤ length E" by simp

  have "real (card (A - B)) = (real (card (sym_diff A B)) + real (card A) - real (card B))/2"
    unfolding card_sym_diff_finite[OF assms(1,2)] using card_sub_int_diff_finite[OF assms(1,2)]
    by simp
  also have "… ≤ (real (length E) + real (card A) - real (card B))/2"
    using assms(5) by simp
  finally have "real (card (A-B)) ≤ dA"
    unfolding a using nat_leq_real_floor by blast
  thus c:"card (A-B) ≤ dA" by auto

  have "real (card (B - A)) = (real (card (sym_diff A B)) + real (card B) - real (card A))/2"
    unfolding card_sym_diff_finite[OF assms(1,2)] using card_sub_int_diff_finite[OF assms(1,2)]
    by simp
  also have "… ≤ (real (length E) + real (card B) - real (card A))/2"
    using assms(5) by simp
  finally have "real (card (B-A)) ≤ dB"
    unfolding b using nat_leq_real_floor by blast
  thus d:"card (B-A) ≤ dB" by auto

  have "real dB - real dA =
    of_int ⌊(real (length E) - card B - card A)/2 + real (card B)⌋ -
    of_int ⌊(real (length E) - card A - card B)/2 + real (card A)⌋"
    unfolding a b by argo
  also have "… = real (card B) - real (card A)"
    by (simp add:algebra_simps)
  also have "… =  real (card (B - A)) - card (A - B)"
    using card_sub_int_diff_finite[OF assms(1,2)] by simp
  finally have "real dB - real dA = real (card (B - A)) - card (A - B)"
    by simp

  thus "dB - card (B - A) = dA - card (A - B)"
    using c d by simp
qed

text "Insert the solution we know that must exist to show it's consistent"
lemma rational_function_interpolation_consistent:
  fixes A B :: "'a::finite_field set"
  assumes
    "∀ x ∈ (set E). f x = fA x / fB x"
    "CARD('m::mod_type) = length E"
    "dA + dB ≤ length E"
    "card (A - B) ≤ dA"
    "card (B - A) ≤ dB"
    "dB - card (B - A) = dA - card (A - B)"
    "∀ x ∈ set E. x ∉ A" "∀ x ∈ set E. x ∉ B"
    "fA = (λ x ∈ set E. poly (set_to_poly A) x)"
    "fB = (λ x ∈ set E. poly (set_to_poly B) x)"
  shows
    "consistent (χ (i::'m) (j::'m). rfi_coefficient_matrix E f dA dB (to_nat i) (to_nat j))
                (χ (i::'m). rfi_constant_vector E f dA dB (to_nat i))"
proof -
  let ?coeff = "rfi_coefficient_matrix E f dA dB"
  let ?const = "rfi_constant_vector E f dA dB"
  let ?coeff' = "(χ (i::'m) (j::'m). ?coeff (to_nat i) (to_nat j))"
  let ?const' = "(χ (i::'m). ?const (to_nat i))"

  define sp where "sp = set_to_poly (A-B) * monom 1 (dA - card (A-B))"
  define sq where "sq = set_to_poly (B-A) * monom 1 (dB - card (B-A))"

  (* We show that ?x is a solution *)

  let ?x = "(χ (i::'m). if (to_nat i) < dA then coeff sp (to_nat i) else coeff sq (to_nat i - dA))"

  have poly_mul_eq: "fA x * poly sq x = fB x * poly sp x" if "x ∈ set E" for x
  proof -
    have "set_to_poly A * set_to_poly (B - A) = (set_to_poly B) * set_to_poly (A - B)"
      by (simp add: Un_commute set_to_poly_mult_distinct)
    then have "(set_to_poly A * set_to_poly (B - A) * monom 1 (dB - card (B - A)) =
       (set_to_poly B) * set_to_poly (A - B) * monom 1 (dA - card (A - B)))"
       using assms(6) by argo
    hence "poly (set_to_poly A) x * poly (set_to_poly (B - A) * monom 1 (dB - card (B - A))) x =
      poly (set_to_poly B) x * poly (set_to_poly (A - B) * monom 1 (dA - card (A - B))) x"
      by (metis (no_types, lifting) mult.commute mult.left_commute poly_mult)
    thus ?thesis
      using that unfolding assms sp_def sq_def by simp
  qed

  have x_sol_raw: "(∑j ∈ {0..<dA}. e ^ j * coeff sp j ) + (∑j ∈ {0..<dB}. - f e * e ^ j * (coeff sq j))
    = f e * e ^ dB - e ^ dA" if "e ∈ set E" for e
  proof -
    have fAz: "fA e ≠ 0"
      using assms (7,9) in_set_to_poly that by auto
    moreover have fBz: "fB e ≠ 0"
      using assms (8,10) in_set_to_poly that by auto
    ultimately have fz: "f e ≠ 0"
      using that assms(1) by simp
    have ffB: "f e = fA e / fB e"
       using that assms(1) by simp

    have "lead_coeff sp = 1"
      unfolding sp_def lead_coeff_mult using set_to_poly_lead_coeff lead_coeff_monom
      by (auto simp add: degree_monom_eq)
    moreover have "degree sp = dA"
      unfolding sp_def using assms(4)
      by (simp add: add_diff_inverse_nat degree_monom_eq degree_mult_eq order_less_imp_not_less set_to_poly_degree)
    ultimately have poly_sp: "poly sp e = (∑i<dA. coeff sp i * e ^ i) + e ^ dA" for e
      unfolding poly_lead_coeff_extract by simp

    have "lead_coeff sq = 1"
      unfolding sq_def lead_coeff_mult using set_to_poly_lead_coeff lead_coeff_monom
      by (auto simp add: degree_monom_eq)
    moreover have "degree sq = dB"
      using assms(5) unfolding sq_def
      by (simp add: degree_monom_eq degree_mult_eq le_eq_less_or_eq set_to_poly_degree)
    ultimately have poly_sq: "poly sq e = (∑i<dB. coeff sq i * e ^ i) + e ^ dB" for e
      unfolding poly_lead_coeff_extract by simp

    have "fB e * ((∑i = 0..<dA. coeff sp i * e ^ i) + e ^ dA) =
          fA e * ((∑i = 0..<dB. coeff sq i * e ^ i) + e ^ dB)"
      using that poly_mul_eq unfolding poly_sq poly_sp lessThan_atLeast0 by simp
    then have "fB e * ((∑j = 0..<dA. e ^ j * coeff sp j) + e ^ dA) =
               fA e * ((∑j = 0..<dB. e ^ j * coeff sq j) + e ^ dB)"
      by (metis (lifting) Finite_Cartesian_Product.sum_cong_aux mult.commute)
    then have "(∑j = 0..<dA. e ^ j * coeff sp j) + e ^ dA =
        f e * ((∑j = 0..<dB. e ^ j * coeff sq j) + e ^ dB)"
      unfolding ffB using fBz
      by (metis (no_types, lifting) fBz nonzero_mult_div_cancel_left times_divide_eq_left)
    also have "… = f e * (∑j = 0..<dB. e ^ j * coeff sq j) + f e * e ^ dB"
      by algebra
    also have "… = (∑j = 0..<dB. f e * e ^ j * coeff sq j) + f e * e ^ dB"
      by (metis (no_types, lifting) Finite_Cartesian_Product.sum_cong_aux mult.assoc sum_distrib_left)
    finally have "(∑j = 0..<dA. e ^ j * coeff sp j) + e ^ dA =
        (∑j = 0..<dB. f e * e ^ j * coeff sq j) + f e * e ^ dB"
      by argo
    then have "(∑j = 0..<dA. e ^ j * coeff sp j) =
        (∑j = 0..<dB. f e * e ^ j * coeff sq j) + f e * e ^ dB - e ^ dA"
      using add_implies_diff by blast
    then have "(∑j = 0..<dA. e ^ j * coeff sp j) + (- (∑j = 0..<dB. f e * e ^ j * coeff sq j)) =
        f e * e ^ dB - e ^ dA"
      by auto
    moreover have "- (∑j = 0..<dB. f e * e ^ j * (coeff sq j)) =
        (∑j = 0..<dB. - f e * e ^ j * coeff sq j)"
      using sum_negf [symmetric] by auto
    ultimately show ?thesis
      by argo
  qed

  (* Now we need to convert it to the matrix product representation *)

  let ?const_lam = "λe. f e * e ^ dB - e ^ dA"
  let ?const_lam' = "λi. ?const_lam (E ! i)"
  let ?coeff_lam = "λe j.(if j < dA then e ^ j
             else if j < dA + dB
                  then - f e * e ^ ( j - dA) else 0) *
            (if j < dA then coeff sp j else coeff sq (j- dA))"
  let ?coeff_lam' = "λi . ?coeff_lam (E ! i)"


  have "(∑j ∈ {0..<length E}. ?coeff_lam e j) = ?const_lam e" if "e ∈ set E" for e
  proof -
    have "(∑j ∈ {0..<length E}. ?coeff_lam e j) = (∑j ∈ {0..<dA + dB}. ?coeff_lam e j)"
      using assms(3) by (intro sum.mono_neutral_cong_right) auto
    also have "… = (∑j∈{0..<dA}. e^j * coeff sp j) + (∑j ∈ {0..<dB}. -f e * e^j*(coeff sq j))"
    proof -
      have "(∑j ∈ {0..<dA + dB}. ?coeff_lam e j) =
          (∑j ∈ {0..<dA}. ?coeff_lam e j) + (∑j ∈ {dA..<dA + dB}. ?coeff_lam e j) "
        by (intro sum.atLeastLessThan_concat [symmetric]) auto
      also have "… = (∑j ∈ {0..<dA}. e ^ j * coeff sp j ) +
         (∑j ∈ {dA..<dA + dB}. - f e * e ^ (j - dA) * (coeff sq (j-dA))) "
        by simp
      moreover have "(∑j ∈ {dA..<dA + dB}. - f e * e ^ (j - dA) * (coeff sq (j-dA))) =
          (∑j ∈ {0..<dB}. - f e * e ^ (j) * (coeff sq j))"
        by (rule sum.reindex_bij_witness [of _ "λi. i + dA" "λi. i - dA"]) auto
      ultimately show ?thesis
        by simp
    qed
    also have "… = ?const_lam e"
      using that x_sol_raw by simp
    finally show ?thesis by simp
  qed
  then have "(∑j ∈ {0..<length E}. ?coeff_lam' i j) = ?const_lam' i" if "i < length E" for i
    using that  by simp
  moreover have "(∑j∈(UNIV::'m set). ?coeff_lam i (to_nat j)) = (∑j ∈ {0..<CARD('m)}. ?coeff_lam i j)" for i
    using bij_to_nat by (intro sum.reindex_bij_betw) blast
  ultimately have "(∑j∈(UNIV::'m set). ?coeff_lam' i (to_nat j)) = ?const_lam' i" if "i < length E" for i
    using that assms using of_nat_eq_iff[of "card top" "length E"] assms(3) by force
  then have "(λi. ∑j∈(UNIV::'m set). ?coeff_lam' i (to_nat j)) (to_nat (i::'m)) = ?const_lam' (to_nat i)" for i
    using mod_type_less_function_eq [of "(λi. ∑j∈(UNIV::'m set). ?coeff_lam' i (to_nat j))" ?const_lam' i]
    using assms(2) assms(3) by auto
  then have eval: "(λi. ∑j∈(UNIV::'m set). ?coeff_lam' (to_nat (i::'m)) (to_nat j)) =
      (λ i. ?const_lam' (to_nat i))"
    by simp

  have "?coeff' *v ?x = ?const'"
    unfolding matrix_vector_mult_def
      rfi_coefficient_matrix_def
      rfi_constant_vector_def
    using eval by simp
  then show ?thesis
    unfolding consistent_def is_solution_def by auto
qed

subsection ‹Main lemma›

lemma rational_function_interpolation_correct:
  assumes
    "int dA = ⌊(real (length E) + card A - card B)/2⌋"
    "int dB = ⌊(real (length E) + card B - card A)/2⌋"
    "card (sym_diff A B) ≤ length E"

    "∀ x ∈ set E. x ∉ A" "∀ x ∈ set E. x ∉ B"
    "fA = (λ x ∈ set E. poly (set_to_poly A) x)"
    "fB = (λ x ∈ set E. poly (set_to_poly B) x)"
    "CARD('m::mod_type) = length E"
  defines
    "sol ≡ solution_to_poly (rational_function_interpolation E (λe. fA e / fB e) dA dB TYPE('m)) dA dB"
  shows
    "monic_interpolated_rational_function (fst sol) (snd sol) (set E) fA fB dA dB"
proof -
  let ?f = "(λe. fA e / fB e)"
  let ?S = "rational_function_interpolation E (λe. fA e / fB e) dA dB TYPE('m)"
  let ?p = "fst (solution_to_poly ?S dA dB)"
  let ?q = "snd (solution_to_poly ?S dA dB)"

  have f:"finite A" "finite B"
    using finite by blast+
  note pd_pq_props = dA_dB_helper[OF f assms(1-3)]

  have "consistent (χ (i::'m) (j::'m). rfi_coefficient_matrix E ?f dA dB (to_nat i) (to_nat j))
        (χ (i::'m). rfi_constant_vector E ?f dA dB (to_nat i))"
    using assms pd_pq_props
    by (intro rational_function_interpolation_consistent [where A = A and B = B and fA = fA and fB = fB])
      auto
  then have "∀e∈set E. fA e * poly ?q e = fB e * poly ?p e"
    using assms pd_pq_props(1) in_set_to_poly
    by (intro rational_function_interpolation_correct_poly [where f = ?f and dA = dA and dB = dB and S = ?S])
      auto
  moreover have "real (degree ?p) = real dA"
    using degree_solution_to_poly_fst by auto
  moreover have "real (degree ?q) = real dB"
    using degree_solution_to_poly_snd by auto
  moreover have "monic ?q"
    using monic_solution_to_poly_snd by auto
  moreover have "monic ?p"
    using monic_solution_to_poly_fst by auto
  ultimately show ?thesis using fst_solution_to_poly_nz snd_solution_to_poly_nz
    unfolding monic_interpolated_rational_function_def sol_def by force
qed

(*Can also interpolate poly with floor of degree*)
lemma interpolated_rational_function_floor_eq:
  "interpolated_rational_function pA pB E fA fB dA dB ⟷
   interpolated_rational_function pA pB E fA fB ⌊dA⌋ ⌊dB⌋"
  unfolding interpolated_rational_function_def using nat_leq_real_floor by simp

lemma sym_diff_bound_div2_ge0:
  fixes A B :: "'a :: finite set"
  assumes "card (sym_diff A B) ≤ length E"
  shows "(real (length E) + card A - card B)/2 ≥ 0"
proof -
  have *: "finite A" "finite B" using finite by auto

  have "0 ≤ real (card (sym_diff A B)) + real (card (A-B)) - (card (B-A))"
    unfolding card_sym_diff_finite[OF *] by simp
  also have "… ≤ real (length E) + real (card (A-B)) - (card (B-A))"
    using assms(1) by simp
  also have "… = (real (length E) + card A - card B)"
    using card_sub_int_diff_finite [OF *] by simp
  finally show ?thesis by simp
qed

text "If the degrees are reals we take the floor first"
lemma rational_function_interpolation_correct_real:
  fixes d'A d'B:: "real"
  assumes
    "card (sym_diff A B) ≤ length E"
    "∀ x ∈ set E. x ∉ A" "∀ x ∈ set E. x ∉ B"
    "fA = (λ x ∈ set E. poly (set_to_poly A) x)"
    "fB = (λ x ∈ set E. poly (set_to_poly B) x)"
    "CARD('m::mod_type) = length E"
  defines "d'A ≡ (real (length E) + card A - card B)/2"
  defines "d'B ≡ (real (length E) + card B - card A)/2"
  defines "dA ≡ nat ⌊d'A⌋"
  defines "dB ≡ nat ⌊d'B⌋"
  defines "sol_poly ≡ interpolate_rat_fun E (λe. fA e / fB e) dA dB TYPE('m)"
  shows
    "monic_interpolated_rational_function (fst sol_poly) (snd sol_poly) (set E) fA fB d'A d'B"
proof -
  have e: "d'A ≥ 0"
    unfolding d'A_def using sym_diff_bound_div2_ge0 assms(1) by auto

  hence a: "int dA = ⌊(real (length E) + real (card A) - real (card B)) / 2⌋"
    using d'A_def unfolding dA_def by simp

  have f: "d'B ≥ 0"
    unfolding d'B_def using sym_diff_bound_div2_ge0 assms(1) by (metis Un_commute)

  hence b: "int dB = ⌊(real (length E) + real (card B) - real (card A)) / 2⌋"
    using d'B_def unfolding dB_def by simp

  have c: "monic_interpolated_rational_function (fst sol_poly) (snd sol_poly) (set E) fA fB dA dB"
    unfolding sol_poly_def interpolate_rat_fun_def
    by (intro rational_function_interpolation_correct [OF a b assms(1-6)])
  moreover have "⌊d'A⌋ = real (nat ⌊d'A⌋)"
    using e by (intro of_nat_nat[symmetric]) simp
  moreover have "⌊d'B⌋ = real (nat ⌊d'B⌋)"
    using f by (intro of_nat_nat[symmetric]) simp
  ultimately have
    "monic_interpolated_rational_function (fst sol_poly) (snd sol_poly) (set E) fA fB (nat ⌊d'A⌋) (nat ⌊d'B⌋)"
    unfolding dA_def dB_def
    by simp
  thus ?thesis unfolding monic_interpolated_rational_function_def
    using assms(9,10) a b d'A_def d'B_def floor_of_nat by simp
qed

end