Theory Roots_via_IA

section ‹Root Filter via Interval Arithmetic›

subsection ‹Generic Framework›

text ‹We provide algorithms for finding all real or complex roots of a polynomial
  from a superset of the roots via interval arithmetic. 
  These algorithms are much faster than just
  evaluating the polynomial via algebraic number computations.›

theory Roots_via_IA
  imports 
    Algebraic_Numbers.Interval_Arithmetic
begin

definition interval_of_real :: "nat ⇒ real ⇒ real interval" where
  "interval_of_real prec x =
      (if is_rat x then Interval x x
       else let n = 2 ^ prec; x' = x * of_int n
            in  Interval (of_rat (Rat.Fract ⌊x'⌋ n)) (of_rat (Rat.Fract ⌈x'⌉ n)))"

definition interval_of_complex :: "nat ⇒ complex ⇒ complex_interval" where
  "interval_of_complex prec z =
     Complex_Interval (interval_of_real prec (Re z)) (interval_of_real prec (Im z))"

fun poly_interval :: "'a :: {plus,times,zero} list ⇒ 'a ⇒ 'a" where
  "poly_interval [] _ = 0"
| "poly_interval [c] _ = c"
| "poly_interval (c # cs) x = c + x * poly_interval cs x"

definition filter_fun_complex :: "complex poly ⇒ nat ⇒ complex ⇒ bool" where
  "filter_fun_complex p = (let c = coeffs p in
      (λ prec. let cs = map (interval_of_complex prec) c
      in (λ x. 0 ∈⇩c poly_interval cs (interval_of_complex prec x))))" 

definition filter_fun_real :: "real poly ⇒ nat ⇒ real ⇒ bool" where
  "filter_fun_real p = (let c = coeffs p in
      (λ prec. let cs = map (interval_of_real prec) c
      in (λ x. 0 ∈⇩i poly_interval cs (interval_of_real prec x))))" 

definition genuine_roots :: "_ poly ⇒ _ list ⇒ _ list" where
  "genuine_roots p xs = filter (λx. poly p x = 0) xs"

lemma zero_in_interval_0 [simp, intro]: "0 ∈⇩i 0"
  unfolding zero_interval_def by auto

lemma zero_in_complex_interval_0 [simp, intro]: "0 ∈⇩c 0"
  unfolding zero_complex_interval_def by (auto simp: in_complex_interval_def)

lemma length_coeffs_degree':
  "length (coeffs p) = (if p = 0 then 0 else Suc (degree p))"
  by (cases "p = 0") (auto simp: length_coeffs_degree)

lemma poly_in_poly_interval_complex:
  assumes "list_all2 (λc ivl. c ∈⇩c ivl) (coeffs p) cs" "x ∈⇩c ivl"
  shows   "poly p x ∈⇩c poly_interval cs ivl"
proof -
  have len_eq: "length (coeffs p) = length cs"
    using assms(1) list_all2_lengthD by blast
  have "coeffs p = map (λi. coeffs p ! i) [0..<length cs]"
    by (subst len_eq [symmetric], rule map_nth [symmetric])
  also have "… = map (poly.coeff p) [0..<length cs]"
    by (intro map_cong) (auto simp: nth_coeffs_coeff len_eq)
  finally have "list_all2 (λc ivl. c ∈⇩c ivl) (map (poly.coeff p) [0..<length cs]) cs"
    using assms by simp
  moreover have "length cs ≥ length (coeffs p)"
    using len_eq by simp
  ultimately show ?thesis using assms(2)
  proof (induction cs ivl arbitrary: p x rule: poly_interval.induct)
    case (1 ivl p x)
    thus ?case by auto
  next
    case (2 c ivl p x)
    have "degree p = 0"
      using 2 by (auto simp: degree_eq_length_coeffs)
    then obtain c' where [simp]: "p = [:c':]"
      by (meson degree_eq_zeroE)
    show ?case using 2 by auto
  next
    case (3 c1 c2 cs ivl p x)
    obtain q c where [simp]: "p = pCons c q"
      by (cases p rule: pCons_cases)
    have "list_all2 in_complex_interval (map (poly.coeff p) [0..<length (c1 # c2 # cs)])
                  (c1 # c2 # cs)"
      using "3.prems"(1) by simp
    also have "[0..<length (c1 # c2 # cs)] = 0 # map Suc [0..<length (c2 # cs)]"
      by (metis length_Cons map_Suc_upt upt_conv_Cons zero_less_Suc)
    also have "map (poly.coeff p) … = c # map (poly.coeff q) [0..<length (c2 # cs)]"
      by auto
    finally have "c ∈⇩c c1" and
        "list_all2 in_complex_interval (map (poly.coeff q) [0..<length (c2 # cs)]) (c2 # cs)"
      using "3.prems" by (simp_all del: upt_Suc)

    have IH: "poly q x ∈⇩c poly_interval (c2 # cs) ivl"
    proof (rule "3.IH")
      show "length (coeffs q) ≤ length (c2 # cs)"
        using "3.prems"(2) unfolding length_coeffs_degree' by auto
    qed fact+

    show ?case
      using IH "3.prems" ‹c ∈⇩c c1›
      by (auto intro!: plus_complex_interval times_complex_interval)
  qed
qed

lemma poly_in_poly_interval_real: fixes x :: real 
  assumes "list_all2 (λc ivl. c ∈⇩i ivl) (coeffs p) cs" "x ∈⇩i ivl"
  shows   "poly p x ∈⇩i poly_interval cs ivl"
proof -
  have len_eq: "length (coeffs p) = length cs"
    using assms(1) list_all2_lengthD by blast
  have "coeffs p = map (λi. coeffs p ! i) [0..<length cs]"
    by (subst len_eq [symmetric], rule map_nth [symmetric])
  also have "… = map (poly.coeff p) [0..<length cs]"
    by (intro map_cong) (auto simp: nth_coeffs_coeff len_eq)
  finally have "list_all2 (λc ivl. c ∈⇩i ivl) (map (poly.coeff p) [0..<length cs]) cs"
    using assms by simp
  moreover have "length cs ≥ length (coeffs p)"
    using len_eq by simp
  ultimately show ?thesis using assms(2)
  proof (induction cs ivl arbitrary: p x rule: poly_interval.induct)
    case (1 ivl p x)
    thus ?case by auto
  next
    case (2 c ivl p x)
    have "degree p = 0"
      using 2 by (auto simp: degree_eq_length_coeffs)
    then obtain c' where [simp]: "p = [:c':]"
      by (meson degree_eq_zeroE)
    show ?case using 2 by auto
  next
    case (3 c1 c2 cs ivl p x)
    obtain q c where [simp]: "p = pCons c q"
      by (cases p rule: pCons_cases)
    have "list_all2 in_interval (map (poly.coeff p) [0..<length (c1 # c2 # cs)])
                  (c1 # c2 # cs)"
      using "3.prems"(1) by simp
    also have "[0..<length (c1 # c2 # cs)] = 0 # map Suc [0..<length (c2 # cs)]"
      by (metis length_Cons map_Suc_upt upt_conv_Cons zero_less_Suc)
    also have "map (poly.coeff p) … = c # map (poly.coeff q) [0..<length (c2 # cs)]"
      by auto
    finally have "c ∈⇩i c1" and
        "list_all2 in_interval (map (poly.coeff q) [0..<length (c2 # cs)]) (c2 # cs)"
      using "3.prems" by (simp_all del: upt_Suc)

    have IH: "poly q x ∈⇩i poly_interval (c2 # cs) ivl"
    proof (rule "3.IH")
      show "length (coeffs q) ≤ length (c2 # cs)"
        using "3.prems"(2) unfolding length_coeffs_degree' by auto
    qed fact+

    show ?case
      using IH "3.prems" ‹c ∈⇩i c1›
      by (auto intro!: plus_in_interval times_in_interval)
  qed
qed


lemma in_interval_of_real [simp, intro]: "x ∈⇩i interval_of_real prec x"
  unfolding interval_of_real_def by (auto simp: Let_def of_rat_rat field_simps)

lemma in_interval_of_complex [simp, intro]: "z ∈⇩c interval_of_complex prec z"
  unfolding interval_of_complex_def in_complex_interval_def by auto

lemma distinct_genuine_roots [simp, intro]: 
  "distinct xs ⟹ distinct (genuine_roots p xs)"
  by (simp add: genuine_roots_def)

definition filter_fun :: "'a poly ⇒ (nat ⇒ 'a :: comm_ring ⇒ bool) ⇒ bool" where
  "filter_fun p f = (∀ n x. poly p x = 0 ⟶ f n x)" 

lemma filter_fun_complex: "filter_fun p (filter_fun_complex p)"
  unfolding filter_fun_def
proof (intro impI allI)
  fix prec x
  assume root: "poly p x = 0" 
  define cs where "cs = map (interval_of_complex prec) (coeffs p)"
  have cs: "list_all2 in_complex_interval (coeffs p) cs"
    unfolding cs_def list_all2_map2 by (intro list_all2_refl in_interval_of_complex)
  define P where "P = (λx. 0 ∈⇩c poly_interval cs (interval_of_complex prec x))"
  have "P x" 
  proof -
    have "poly p x ∈⇩c poly_interval cs (interval_of_complex prec x)"
      by (intro poly_in_poly_interval_complex in_interval_of_complex cs)
    with root show ?thesis
      by (simp add: P_def)
  qed
  thus "filter_fun_complex p prec x" unfolding filter_fun_complex_def Let_def P_def
    using cs_def by blast
qed

lemma filter_fun_real: "filter_fun p (filter_fun_real p)"
  unfolding filter_fun_def
proof (intro impI allI)
  fix prec x
  assume root: "poly p x = 0" 
  define cs where "cs = map (interval_of_real prec) (coeffs p)"
  have cs: "list_all2 in_interval (coeffs p) cs"
    unfolding cs_def list_all2_map2 by (intro list_all2_refl in_interval_of_real)
  define P where "P = (λx. 0 ∈⇩i poly_interval cs (interval_of_real prec x))"
  have "P x" 
  proof -
    have "poly p x ∈⇩i poly_interval cs (interval_of_real prec x)"
      by (intro poly_in_poly_interval_real in_interval_of_real cs)
    with root show ?thesis
      by (simp add: P_def)
  qed
  thus "filter_fun_real p prec x" unfolding filter_fun_real_def Let_def P_def
    using cs_def by blast
qed

context
  fixes p :: "'a :: comm_ring poly" and f
  assumes ff: "filter_fun p f" 
begin

lemma genuine_roots_step:
  "genuine_roots p xs = genuine_roots p (filter (f prec) xs)"
  unfolding genuine_roots_def filter_filter 
  using ff[unfolded filter_fun_def, rule_format, of _ prec] by metis 

lemma genuine_roots_step_preserve_invar:
  assumes "{z. poly p z = 0} ⊆ set xs"
  shows   "{z. poly p z = 0} ⊆ set (filter (f prec) xs)"
proof -
  have "{z. poly p z = 0} = set (genuine_roots p xs)"
    using assms by (auto simp: genuine_roots_def)
  also have "… = set (genuine_roots p (filter (f prec) xs))"
    using genuine_roots_step[of _ prec] by simp
  also have "… ⊆ set (filter (f prec) xs)"
    by (auto simp: genuine_roots_def)
  finally show ?thesis .
qed
end

lemma genuine_roots_finish:
  fixes p :: "'a :: field_char_0 poly" 
  assumes "{z. poly p z = 0} ⊆ set xs" "distinct xs"
  assumes "length xs = card {z. poly p z = 0}"
  shows   "genuine_roots p xs = xs"
proof -
  have [simp]: "p ≠ 0"
    using finite_subset[OF assms(1) finite_set] infinite_UNIV_char_0 by auto
  have "length (genuine_roots p xs) = length xs"
    unfolding genuine_roots_def using assms 
    by (simp add: Int_absorb2 distinct_length_filter)
  thus ?thesis
    unfolding genuine_roots_def
    by (metis filter_True length_filter_less linorder_not_less order_eq_iff)
qed

text ‹This is type of the initial search problem. It consists of a polynomial $p$, 
  a list $xs$ of candidate roots, the cardinality of the set of roots of $p$ and a filter function to
  drop non-roots that is parametric in a precision parameter.›
typedef (overloaded) 'a genuine_roots_aux =
  "{(p :: 'a :: field_char_0 poly, xs, n, ff). 
    distinct xs ∧ 
    {z. poly p z = 0} ⊆ set xs ∧ 
    card {z. poly p z = 0} = n ∧
    filter_fun p ff}"
  by (rule exI[of _ "(1, [], 0, λ _ _. False)"], auto simp: filter_fun_def)

setup_lifting type_definition_genuine_roots_aux

lift_definition genuine_roots' :: "nat ⇒ 'a :: field_char_0 genuine_roots_aux ⇒ 'a list" is
  "λprec (p, xs, n, ff). genuine_roots p xs" .

lift_definition genuine_roots_impl_step' :: "nat ⇒ 'a :: field_char_0 genuine_roots_aux ⇒ 'a genuine_roots_aux" is
  "λprec (p, xs, n, ff). (p, filter (ff prec) xs, n, ff)"
  by (safe, intro distinct_filter, auto simp: filter_fun_def)

lift_definition gr_poly :: "'a :: field_char_0 genuine_roots_aux ⇒ 'a poly" is
  "λ(p :: 'a poly, _, _, _). p" .

lift_definition gr_list :: "'a :: field_char_0 genuine_roots_aux ⇒ 'a list" is
  "λ(_, xs :: 'a list, _, _). xs" .

lift_definition gr_numroots :: "'a :: field_char_0 genuine_roots_aux ⇒ nat" is
  "λ(_, _, n, _). n" .

lemma genuine_roots'_code [code]:
  "genuine_roots' prec gr =
     (if length (gr_list gr) = gr_numroots gr then gr_list gr
      else genuine_roots' (2 * prec) (genuine_roots_impl_step' prec gr))"
proof (transfer, clarify)
  fix prec :: nat and p :: "'a poly" and xs :: "'a list" and ff
  assume *: "{z. poly p z = 0} ⊆ set xs" "distinct xs" "filter_fun p ff" 
  show "genuine_roots p xs =
          (if length xs = card {z. poly p z = 0} then xs
           else genuine_roots p (filter (ff prec) xs))"
    using genuine_roots_finish[of p xs] genuine_roots_step[of p] * by auto
qed

definition initial_precision :: nat where "initial_precision = 10" 

definition genuine_roots_impl :: "'a genuine_roots_aux ⇒ 'a :: field_char_0 list" where
  "genuine_roots_impl = genuine_roots' initial_precision" 

lemma genuine_roots_impl: "set (genuine_roots_impl p) = {z. poly (gr_poly p) z = 0}" 
  "distinct (genuine_roots_impl p)" 
  unfolding genuine_roots_impl_def
  by (transfer, auto simp: genuine_roots_def, transfer, auto)

end