Theory BenOr_Kozen_Reif.BKR_Decision

theory BKR_Decision
  imports BKR_Algorithm
    "Berlekamp_Zassenhaus.Factorize_Rat_Poly"
    "Algebraic_Numbers.Real_Roots"
     BKR_Proofs
    "HOL.Deriv"
begin

section "Algorithm"

subsection "Parsing"
  (* Formula type *)
datatype 'a fml =
  And "'a fml" "'a fml"
  | Or "'a fml" "'a fml"
  | Gt 'a   (* 'a > 0 *)
  | Geq 'a  (* 'a ≥ 0 *)
  | Lt 'a   (* 'a < 0 *)
  | Leq 'a  (* 'a ≤ 0 *)
  | Eq 'a   (* 'a = 0 *)
  | Neq 'a  (* 'a ≠ 0 *)

(* Evaluating a formula over a lookup semantics where 'a is nat *)
primrec lookup_sem :: "nat fml  ('a::linordered_field list)  bool"
  where
    "lookup_sem (And l r) ls = (lookup_sem l ls  lookup_sem r ls)"
  | "lookup_sem (Or l r) ls = (lookup_sem l ls  lookup_sem r ls)"
  | "lookup_sem (Gt p) ls = (ls ! p > 0)"  
  | "lookup_sem (Geq p) ls = (ls ! p  0)" 
  | "lookup_sem (Lt p) ls = (ls ! p < 0)" 
  | "lookup_sem (Leq p) ls = (ls ! p  0)" 
  | "lookup_sem (Eq p) ls = (ls ! p = 0)" 
  | "lookup_sem (Neq p) ls = (ls ! p  0)" 

(* (compute) all polynomials mentioned in a formula *)
primrec poly_list :: "'a fml  'a list"
  where
    "poly_list (And l r) = poly_list l @ poly_list r"  
  | "poly_list (Or l r) = poly_list l @ poly_list r"  
  | "poly_list (Gt p) = [p]"
  | "poly_list (Geq p) = [p]"
  | "poly_list (Lt p) = [p]"
  | "poly_list (Leq p) = [p]"
  | "poly_list (Eq p) = [p]"
  | "poly_list (Neq p) = [p]"

primrec index_of_aux :: "'a list  'a  nat  nat" where
  "index_of_aux [] y n = n"
| "index_of_aux (x#xs) y n =
  (if x = y then n else index_of_aux xs y (n+1))"

definition index_of :: "'a list  'a  nat" where
  "index_of xs y = index_of_aux xs y 0"

definition convert :: "'a fml  (nat fml × 'a list)"
  where
    "convert fml = (
    let ps = remdups (poly_list fml)
    in
    (map_fml (index_of ps) fml, ps)
  )"


subsection "Factoring"

(* Makes sure the result of factorize_rat_poly is monic *)
definition factorize_rat_poly_monic :: "rat poly  (rat × (rat poly × nat) list)"
  where
    "factorize_rat_poly_monic p = (
    let (c,fs) =  factorize_rat_poly p ;
        lcs = prod_list (map (λ(f,i). (lead_coeff f) ^ i) fs) ;
        fs = map (λ(f,i). (normalize f, i)) fs
    in
    (c * lcs,fs)
  )"

(* Factoring an input list of polynomials *)
definition factorize_polys :: "rat poly list  (rat poly list × (rat × (nat × nat) list) list)"
  where
    "factorize_polys ps = (
    let fact_ps = map factorize_rat_poly_monic ps;
        factors = remdups (map fst (concat (map snd fact_ps))) ;
        data = map (λ(c,fs). (c, map (λ(f,pow). (index_of factors f, pow) ) fs)) fact_ps
    in
      (factors,data)
  )"

(* After turning a polynomial into factors,
  this turns a sign condition on the factors
  into a sign condition for the polynomial *)
definition undo_factorize :: "rat × (nat × nat) list  rat list  rat"
  where
    "undo_factorize cfs signs =
   squash
    (case cfs of (c,fs) 
    (c * prod_list (map (λ(f,pow). (signs ! f) ^ pow) fs)))
  "

definition undo_factorize_polys :: "(rat × (nat × nat) list) list  rat list  rat list"
  where
    "undo_factorize_polys ls signs = map (λl. undo_factorize l signs) ls"

subsection "Auxiliary Polynomial"
definition crb:: "real poly  int" where
  "crb p = ceiling (2 + max_list_non_empty (map (λ i. norm (coeff p i)) [0 ..< degree p]) 
    / norm (lead_coeff p))"

(* Because we are using prod_list instead of lcm, it's important that this is called 
    when ps is pairwise coprime. *)
definition coprime_r :: "real poly list  real poly"
  where
    "coprime_r ps = pderiv (prod_list ps) * ([:-(crb (prod_list ps)),1:]) * ([:(crb (prod_list ps)),1:])"


subsection "Setting Up the Procedure"
  (* 0 indexed *)
definition insertAt :: "nat  'a   'a list  'a list" where
  "insertAt n x ls = take n ls @ x # (drop n ls)"

(* 0 indexed *)
definition removeAt :: "nat  'a list  'a list" where
  "removeAt n ls = take n ls @ (drop (n+1) ls)"

definition find_sgas_aux:: "real poly list  rat list list"
  where "find_sgas_aux in_list =
  concat (map (λi.
    map (λv. insertAt i 0 v) (find_consistent_signs_at_roots (in_list ! i) (removeAt i in_list))
  ) [0..<length in_list])"

(* For an input list of real polynomials, apply BKR to all positions *)
definition find_sgas :: "real poly list  rat list list"
  where
    "find_sgas ps = ( 
    let r = coprime_r ps in
    find_consistent_signs_at_roots r ps @ find_sgas_aux ps
  )"

(* Putting the sign condition preprocessing together with BKR *)
definition find_consistent_signs :: "rat poly list  rat list list"
  where
    "find_consistent_signs ps = (
     let (fs,data) = factorize_polys ps;
         sgas = find_sgas (map (map_poly of_rat) fs);
         rsgas = map (undo_factorize_polys data) sgas
    in
    (if fs = [] then [(map (λx. if poly x 0 < 0 then -1 else if poly x 0 = 0 then 0 else 1) ps)] else rsgas)
  )"


subsection "Deciding Univariate Problems"
definition decide_universal :: "rat poly fml  bool"
  where [code]:
    "decide_universal fml = (
    let (fml_struct,polys) = convert fml;
        conds = find_consistent_signs polys
    in
     list_all (lookup_sem fml_struct) conds
  )"

definition decide_existential :: "rat poly fml  bool"
  where [code]:
    "decide_existential fml = (
    let (fml_struct,polys) = convert fml;
        conds = find_consistent_signs polys
    in
      find (lookup_sem fml_struct) conds  None
  )"

section "Proofs"
subsection "Parsing and Semantics"
  (* Evaluating a formula where 'a is a real poly *)
primrec real_sem :: "real poly fml  real  bool"
  where
    "real_sem (And l r) x = (real_sem l x  real_sem r x)"
  | "real_sem (Or l r) x = (real_sem l x  real_sem r x)"
  | "real_sem (Gt p) x = (poly p x > 0)" 
  | "real_sem (Geq p) x = (poly p x  0)" 
  | "real_sem (Lt p) x = (poly p x < 0)" 
  | "real_sem (Leq p) x = (poly p x  0)" 
  | "real_sem (Eq p) x = (poly p x = 0)" 
  | "real_sem (Neq p) x = (poly p x  0)" 

(* Evaluating a formula where 'a is a rat poly *)
primrec fml_sem :: "rat poly fml  real  bool"
  where
    "fml_sem (And l r) x = (fml_sem l x  fml_sem r x)"
  | "fml_sem (Or l r) x = (fml_sem l x  fml_sem r x)"
  | "fml_sem (Gt p) x = (rpoly p x > 0)" 
  | "fml_sem (Geq p) x = (rpoly p x  0)" 
  | "fml_sem (Lt p) x = (rpoly p x < 0)" 
  | "fml_sem (Leq p) x = (rpoly p x  0)" 
  | "fml_sem (Eq p) x = (rpoly p x = 0)" 
  | "fml_sem (Neq p) x = (rpoly p x  0)" 

lemma poly_list_set_fml:
  shows "set (poly_list fml) = set_fml fml"
  apply (induction) by auto

lemma convert_semantics_lem:
  assumes "p. p  set (poly_list fml) 
    ls ! (index_of ps p) = rpoly p x"
  shows "fml_sem fml x = lookup_sem (map_fml (index_of ps) fml) ls"
  using assms apply (induct fml)
  by auto

lemma index_of_aux_more:
  shows "index_of_aux ls p n  n"
  apply (induct ls arbitrary: n)
  apply auto
  using Suc_leD by blast

lemma index_of_aux_lookup:
  assumes "p  set ls"
  shows "(index_of_aux ls p n) - n < length ls"
    "ls ! ((index_of_aux ls p n) - n) = p"
  using assms apply (induct ls arbitrary: n)
  apply auto
  apply (metis Suc_diff_Suc index_of_aux_more lessI less_Suc_eq_0_disj less_le_trans)
  by (metis Suc_diff_Suc index_of_aux_more lessI less_le_trans nth_Cons_Suc)

lemma index_of_lookup:
  assumes "p  set ls"
  shows "index_of ls p < length ls"
    "ls ! (index_of ls p) = p"
  apply (metis assms index_of_aux_lookup(1) index_of_def minus_nat.diff_0)
  by (metis assms index_of_aux_lookup(2) index_of_def minus_nat.diff_0)

lemma convert_semantics:
  shows "fml_sem fml x = lookup_sem (fst (convert fml)) (map (λp. rpoly p x) (snd (convert fml)))"
  unfolding convert_def Let_def apply simp
  apply (intro convert_semantics_lem)
  by (simp add: index_of_lookup(1) index_of_lookup(2))

lemma convert_closed:
  shows "i. i  set_fml (fst (convert fml))  i < length (snd (convert fml))"
  unfolding convert_def Let_def
  apply (auto simp add: fml.set_map)
  by (simp add: index_of_lookup(1) poly_list_set_fml)

(* Rational sign vector of polynomials qs with rational coefficients at x *)
definition sign_vec::"rat poly list  real  rat list"
  where "sign_vec qs x 
    map (squash  (λp. rpoly p x)) qs"

(* The set of all rational sign vectors for qs wrt the set S
  When S = UNIV, then this quantifies over all reals *)
definition consistent_sign_vectors::"rat poly list  real set  rat list set"
  where "consistent_sign_vectors qs S = (sign_vec qs) ` S"

lemma sign_vec_semantics:
  assumes "i. i  set_fml fml  i < length ls"
  shows "lookup_sem fml (map (λp. rpoly p x) ls) = lookup_sem fml (sign_vec ls x)"
  using assms apply (induction)
  by (auto simp add: sign_vec_def squash_def)

(* The universal and existential decision procedure is easy if we know the consistent sign vectors *)
lemma universal_lookup_sem:
  assumes "i. i  set_fml fml  i < length qs"
  assumes "set signs = consistent_sign_vectors qs UNIV"
  shows "(x::real. lookup_sem fml (map (λp. rpoly p x) qs)) 
    list_all (lookup_sem fml) signs"
  using assms(2) unfolding consistent_sign_vectors_def list_all_iff
  by (simp add: assms(1) sign_vec_semantics)

lemma existential_lookup_sem:
  assumes "i. i  set_fml fml  i < length qs"
  assumes "set signs = consistent_sign_vectors qs UNIV"
  shows "(x::real. lookup_sem fml (map (λp. rpoly p x) qs)) 
    find (lookup_sem fml) signs  None"
  using assms(2) unfolding consistent_sign_vectors_def find_None_iff
  by (simp add: assms(1) sign_vec_semantics)

subsection "Factoring Lemmas"
  (*definition real_factorize_list:: "rat poly list ⇒ real poly list"
where "real_factorize_list qs = map (map_poly of_rat) (fst(factorize_polys qs))"
*)
interpretation of_rat_poly_hom: map_poly_comm_semiring_hom of_rat..
interpretation of_rat_poly_hom: map_poly_comm_ring_hom of_rat..
interpretation of_rat_poly_hom: map_poly_idom_hom of_rat..

lemma finite_prod_map_of_rat_poly_hom:
  shows "poly (real_of_rat_poly ((a,b)s. f a b)) y = ((a,b)s. poly (real_of_rat_poly (f a b)) y)"
  apply (simp add: of_rat_poly_hom.hom_prod poly_prod)
  by (simp add: case_prod_app prod.case_distrib)

lemma sign_vec_index_of:
  assumes "f  set ftrs"
  shows "sign_vec ftrs x ! (index_of ftrs f) = squash (rpoly f x)"
  by (simp add: assms index_of_lookup(1) index_of_lookup(2) sign_vec_def)

lemma squash_idem:
  shows "squash (squash x) = squash x"
  unfolding squash_def by auto

lemma squash_mult:
  shows "squash ((a::real) * b) = squash a * squash b"
  unfolding squash_def apply auto
  using less_not_sym mult_neg_neg apply blast
  using mult_less_0_iff by blast

lemma squash_prod_list:
  shows "squash (prod_list (ls::real list)) = prod_list (map squash ls)"
  apply (induction ls)
  unfolding squash_def apply auto
  apply (simp add: mult_less_0_iff)
  by (simp add: zero_less_mult_iff)

lemma squash_pow:
  shows "squash ((x::real) ^ (y::nat)) = (squash x) ^ y"
  unfolding squash_def apply auto
  by (auto simp add: zero_less_power_eq)

lemma squash_real_of_rat[simp]:
  shows "squash (real_of_rat x) = squash x"
  unfolding squash_def by auto

lemma factorize_rat_poly_monic_irreducible_monic:
  assumes "factorize_rat_poly_monic f = (c,fs)"
  assumes "(fi,i)  set fs"
  shows "irreducible fi  monic fi"
proof -
  obtain c' fs' where cfs: "factorize_rat_poly f = (c',fs')"
    by (meson surj_pair)
  then have fs: "fs = map (λ(f,i). (normalize f, i)) fs'"
    using factorize_rat_poly_monic_def assms by auto
  obtain "fi'" where "(fi',i)  set fs'" "fi = normalize fi'"
    using assms(2) unfolding fs by auto
  thus ?thesis using factorize_rat_poly irreducible_normalize_iff
    by (metis cfs monic_normalize not_irreducible_zero)
qed

lemma square_free_normalize:
  assumes "square_free p"
  shows "square_free (normalize p)"
  by (metis assms square_free_multD(3) unit_factor_mult_normalize)

lemma coprime_normalize:
  assumes "coprime a b"
  shows "coprime (normalize a) b"
  using assms by auto

lemma undo_normalize:
  shows "a = Polynomial.smult (unit_factor (lead_coeff a)) (normalize a)"
  by (metis add.right_neutral mult_pCons_right mult_zero_right normalize_mult_unit_factor pCons_0_hom.hom_zero unit_factor_poly_def)

lemma finite_smult_distr:
  assumes "distinct fs"
  shows "((x,y)set fs. Polynomial.smult ((f x y)::rat) (g x y)) = 
    Polynomial.smult ((x,y)set fs. f x y) ((x,y)set fs. g x y)"
  using assms
proof (induction fs)
  case Nil
  then show ?case by auto
next
  case (Cons a fs)
  then show ?case apply auto
    using mult.commute mult_smult_right prod.case_distrib smult_smult split_cong split_conv
    by (simp add: Groups.mult_ac(2) split_beta)
qed

lemma normalize_coprime_degree:
  assumes "normalize (f::rat poly) = normalize g"
  assumes "coprime f g"
  shows "degree f = 0"
proof -
  have "f dvd g" by (simp add: assms(1) associatedD2) 
  then have "f dvd 1"
    using assms(2) associatedD1 by auto
  thus ?thesis
    using Missing_Polynomial_Factorial.is_unit_field_poly by blast
qed

lemma factorize_rat_poly_monic_square_free_factorization:
  assumes res: "factorize_rat_poly_monic f = (c,fs)"
  shows "square_free_factorization f (c,fs)"
proof (unfold square_free_factorization_def split, intro conjI impI allI)
  obtain c' fs' where cfs: "factorize_rat_poly f = (c',fs')"
    by (meson surj_pair)
  then have fs: "fs = map (λ(f,i). (normalize f, i)) fs'"
    using factorize_rat_poly_monic_def assms by auto
  have sq: "square_free_factorization f (c',fs')"
    using cfs factorize_rat_poly(1) by blast
  obtain lcs where lcs: "lcs = prod_list (map (λ(f,i). lead_coeff f ^ i) fs')" by force
  have c: "c = c' * lcs"  using assms unfolding factorize_rat_poly_monic_def cfs Let_def lcs by auto
  show "f = 0  c = 0" using c cfs by auto                       
  show "f = 0  fs = []" using fs cfs by auto
  have dist: "distinct fs'" using sq square_free_factorizationD(5) by blast
  show dist2: "distinct fs" unfolding fs
    unfolding distinct_conv_nth apply auto
  proof -
    fix i j
    assume ij: "i < length fs'" "j < length fs'" "i  j"
    assume eq: "(case fs' ! i of
            (f, x)  (normalize f, x)) =
           (case fs' ! j of
            (f, x)  (normalize f, x))"
    obtain f a where fa: "fs' ! i = (f,a)" by force
    obtain g where g: "fs' ! j = (g,a)" "normalize f = normalize g"
      using eq fa  apply auto
      by (metis case_prod_conv prod.collapse prod.inject)
    have "f  g" using dist ij fa g
      using nth_eq_iff_index_eq by fastforce
    then have "coprime f g"
      using square_free_factorizationD(3)[OF sq, of f a g a] fa g ij
      apply auto
      using nth_mem by force
    then have "degree f = 0"
      by (simp add: g(2) normalize_coprime_degree)
    thus False
      using fa ij(1) nth_mem sq square_free_factorizationD'(3) by fastforce 
  qed
  have ceq: "c = c' * ((a, i)set fs'. (lead_coeff a) ^ i)" using c lcs
    by (simp add: dist prod.distinct_set_conv_list) 
  have fseq: " ((a, i)set fs. a ^ i) =  ((a, i)set fs'. (normalize a) ^ i)"
    apply (subst prod.distinct_set_conv_list[OF dist])
    apply (subst prod.distinct_set_conv_list[OF dist2])
    unfolding fs apply (auto simp add: o_def )
    by (metis (no_types, lifting) case_prod_conv old.prod.exhaust)

  have "f = Polynomial.smult c' ((a, i)set fs'. a ^ i)"  using sq square_free_factorizationD(1) by blast
  moreover have "... =  Polynomial.smult c' ((a, i)set fs'. (Polynomial.smult ((unit_factor (lead_coeff a))) (normalize a)) ^ i)"
    apply (subst undo_normalize[symmetric]) by auto
  moreover have "... =  Polynomial.smult c'
    ((a, i)set fs'. (Polynomial.smult ((lead_coeff a) ^ i) ((normalize a) ^ i)))"
    apply (subst smult_power) by auto
  moreover have "... =  Polynomial.smult c'
    (Polynomial.smult ((a, i)set fs'. ((lead_coeff a) ^ i))
        ((a, i)set fs'. (normalize a) ^ i))"
    apply (subst finite_smult_distr) by (auto simp add: dist)
  moreover have "... =  Polynomial.smult (c' * ((a, i)set fs'. (lead_coeff a) ^ i))
        ((a, i)set fs'. (normalize a) ^ i)"
    using smult_smult by blast
  moreover have "... = Polynomial.smult c  ((a, i)set fs. a ^ i)"
    unfolding ceq fseq by auto
  ultimately show "f =  Polynomial.smult c  ((a, i)set fs. a ^ i)" by auto
  fix a i
  assume ai: "(a,i)  set fs"
  obtain a' where a': "(a',i)  set fs'" "a = normalize a'" using ai unfolding fs by auto
  show "square_free a" using square_free_normalize a'
    using sq square_free_factorizationD(2) by blast
  show "0 < degree a" "0 < i" using degree_normalize a'
    using sq square_free_factorizationD'(3) by fastforce+
  fix  b j
  assume bj: "(b,j)  set fs" "(a,i)  (b,j)"
  obtain b' where b': "(b',j)  set fs'" "b = normalize b'" using bj unfolding fs by auto
  show "algebraic_semidom_class.coprime a b" using a' b' apply auto
    using bj(2) sq square_free_factorizationD(3) by fastforce
qed  

lemma undo_factorize_correct:
  assumes "factorize_rat_poly_monic p = (c,fs)"
  assumes "f p. (f,p)  set fs  f  set ftrs"
  shows "undo_factorize (c,map (λ(f,pow). (index_of ftrs f, pow)) fs) (sign_vec ftrs x) = squash (rpoly p x)"
proof -
  have p: "p = smult c ((a, i) set fs. a ^ i)"
    using assms(1) factorize_rat_poly_monic_square_free_factorization square_free_factorizationD(1) by blast
  have fs: "distinct fs"
    using assms(1) factorize_rat_poly_monic_square_free_factorization square_free_factorizationD(5) by blast
  have "rpoly p x = ((real_of_rat c) * rpoly ((a, i) set fs. a ^ i) x)"
    using p by (simp add: of_rat_hom.map_poly_hom_smult)
  moreover have "... = ((real_of_rat c) * rpoly (ai set fs. case ai of (a,i)  a ^ i) x)"
    by blast
  moreover have "... = ((real_of_rat c) * (ai set fs. case ai of (a,i)  rpoly (a ^ i) x))"
    by (simp add: finite_prod_map_of_rat_poly_hom)
  moreover have "... = ((real_of_rat c) * (ai set fs. case ai of (a,i)  (rpoly a x) ^ i))"
    by (metis (mono_tags, lifting) of_rat_poly_hom.hom_power poly_hom.hom_power split_cong)
  moreover have "...  = ((real_of_rat c) * (prod_list (map (λai. case ai of (a,i)  (rpoly a x) ^ i) fs)))"
    by (simp add: fs prod.distinct_set_conv_list)
  ultimately have "rpoly p x = ((real_of_rat c) * (prod_list (map (λai. case ai of (a,i)  (rpoly a x) ^ i) fs)))" by auto

  then have "squash (rpoly p x) = squash c * prod_list (map squash (map (λai. case ai of (a,i)  (rpoly a x) ^ i) fs))"
    by (auto simp add: squash_mult squash_prod_list o_def)
  moreover have "... = squash c * prod_list (map (λai. case ai of (a,i)  squash ((rpoly a x) ^ i)) fs)"
    apply (simp add: o_def)
    by (simp add: prod.case_distrib)
  ultimately have rp:"squash(rpoly p x) = squash c * prod_list (map (λai. case ai of (a,i)  squash (rpoly a x) ^ i) fs)"
    using squash_pow
    by presburger
  have "undo_factorize
     (c, map (λ(f, pow).(index_of ftrs f, pow)) fs) (sign_vec ftrs x) =
    squash
     (c * (xafs. case xa of (f, y)  sign_vec ftrs x ! index_of ftrs f ^ y))"
    unfolding undo_factorize_def apply (auto simp add: o_def)
    by (metis (mono_tags, lifting) case_prod_conv old.prod.exhaust)
  moreover have "... =  squash
     (c * (xafs. case xa of (f, y)  (squash (rpoly f x)) ^ y))"
    using assms(2) sign_vec_index_of map_eq_conv split_cong by (smt (verit, del_insts))
  ultimately show ?thesis using rp
    by (metis (mono_tags, lifting) of_rat_hom.hom_mult squash_idem squash_mult squash_real_of_rat)
qed

lemma length_sign_vec[simp]:
  shows "length (sign_vec ps x) = length ps" unfolding sign_vec_def by auto

lemma factorize_polys_has_factors:
  assumes "factorize_polys ps = (ftrs,data)"
  assumes "p  set ps"
  assumes "factorize_rat_poly_monic p = (c,fs)"
  shows "set (map fst fs)  set ftrs"
  using assms unfolding factorize_polys_def Let_def apply auto
  by (metis UN_iff fst_conv image_eqI snd_conv)

lemma factorize_polys_undo_factorize_polys:
  assumes "factorize_polys ps = (ftrs,data)"
  shows "undo_factorize_polys data (sign_vec ftrs x) = sign_vec ps x"
  unfolding list_eq_iff_nth_eq undo_factorize_polys_def apply auto
proof -
  show leq:"length data = length ps"
    using assms unfolding factorize_polys_def by (auto simp add: Let_def)
  fix i
  assume il:"i < length data"
  obtain c fs where cfs: "factorize_rat_poly_monic (ps ! i) = (c,fs)"
    by (meson surj_pair)
  then have fsts:"set (map fst fs)  set ftrs"
    using assms factorize_polys_has_factors il leq nth_mem by fastforce
  have *:"data ! i = (c,map (λ(f,pow). (index_of ftrs f, pow)) fs)"
    using assms unfolding factorize_polys_def
    using cfs il by (auto simp add: Let_def cfs)
  have "undo_factorize (data ! i) (sign_vec ftrs x) = squash (rpoly (ps ! i) x)" unfolding *
    apply (subst  undo_factorize_correct[of "ps ! i"])
    apply (auto simp add: cfs)
    using fsts by auto
  thus "undo_factorize (data ! i) (sign_vec ftrs x) =  sign_vec ps x ! i" 
    using leq il sign_vec_def by auto
qed

lemma factorize_polys_irreducible_monic:
  assumes "factorize_polys ps = (fs,data)"
  shows "distinct fs" "f. f  set fs  irreducible f  monic f"
  using assms unfolding factorize_polys_def Let_def apply auto
  using factorize_rat_poly_monic_irreducible_monic
  apply (metis prod.collapse) 
  using factorize_rat_poly_monic_irreducible_monic
  by (metis prod.collapse)

lemma factorize_polys_square_free:
  assumes "factorize_polys ps = (fs,data)"
  shows "f. f  set fs  square_free f"
  using assms factorize_polys_irreducible_monic(2) irreducible_imp_square_free by blast

lemma irreducible_monic_coprime:
  assumes f: "monic f" "irreducible (f::rat poly)" 
  assumes g: "monic g" "irreducible (g::rat poly)"
  assumes "f  g"
  shows "coprime f g"
  by (metis (no_types, lifting) assms(5) coprime_0(2) coprime_def' f(1) f(2) g(1) g(2) irreducible_normalized_divisors normalize_dvd_iff normalize_idem normalize_monic)

lemma factorize_polys_coprime:
  assumes "factorize_polys ps = (fs,data)"
  shows "f g. f  set fs  g  set fs  f  g  coprime f g"
  using assms factorize_polys_irreducible_monic(2) irreducible_monic_coprime by auto

lemma coprime_rat_poly_real_poly:
  assumes "coprime p (q::rat poly)"
  shows "coprime (real_of_rat_poly p) ((real_of_rat_poly q)::real poly)"
  by (metis assms gcd_dvd_1 of_rat_hom.map_poly_gcd of_rat_poly_hom.hom_dvd_1)

lemma coprime_rat_poly_iff_coprimereal_poly:
  shows "coprime p (q::rat poly)  coprime (real_of_rat_poly p) ((real_of_rat_poly q)::real poly)"
proof -
  have forward: "coprime p (q::rat poly)  coprime (real_of_rat_poly p) ((real_of_rat_poly q)::real poly)"
    using coprime_rat_poly_real_poly by auto
  have backward: "coprime (real_of_rat_poly p) ((real_of_rat_poly q)::real poly)  coprime p (q::rat poly)"
  proof -
    assume copr_real: "comm_monoid_mult_class.coprime (real_of_rat_poly p) (real_of_rat_poly q)"
    have "degree (gcd p (q::rat poly)) > 0  False" 
    proof -
      assume deg: "degree (gcd p (q::rat poly)) > 0"
      then have "y. y dvd p  y dvd q  degree y > 0"
        by blast
      then obtain y where yprop: "y dvd p  y dvd q  degree y > 0"
        by auto
      then have "(real_of_rat_poly y) dvd (real_of_rat_poly p) 
        (real_of_rat_poly y ) dvd (real_of_rat_poly q)  degree y > 0"
        by simp
      then show "False" 
        using copr_real apply (auto)
        by fastforce 
    qed
    then show "comm_monoid_mult_class.coprime p (q::rat poly)" 
      using comm_monoid_gcd_class.gcd_dvd_1
      by (metis Missing_Polynomial_Factorial.is_unit_field_poly copr_real gcd_zero_iff' neq0_conv of_rat_poly_hom.hom_zero)
  qed
  show ?thesis
    using forward backward by auto
qed

lemma factorize_polys_map_distinct:
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss =  map real_of_rat_poly fs"
  shows "distinct fss"
  using factorize_polys_irreducible_monic[OF assms(1)]
  unfolding assms(2) 
  apply (simp add: distinct_conv_nth)
  by (metis of_rat_eq_iff of_rat_hom.coeff_map_poly_hom poly_eqI)

lemma factorize_polys_map_square_free:
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss =  map real_of_rat_poly fs"
  shows "f. f  set fss  square_free f"
  using factorize_polys_square_free[OF assms(1)]
  using assms(2) field_hom_0'.square_free_map_poly of_rat_hom.field_hom_0'_axioms by auto 

lemma factorize_polys_map_coprime:
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss =  map real_of_rat_poly fs"
  shows "f g. f  set fss  g  set fss  f  g  coprime f g"
  using factorize_polys_coprime[OF assms(1)] coprime_rat_poly_real_poly unfolding assms(2)
  by auto

lemma coprime_prod_list:
  assumes "p. p  set ps  p  0"
  assumes "coprime (prod_list ps) (q::real poly)"  
  shows "p. p  set ps  coprime p q"
proof -
  fix p
  assume "p  set ps"
  then obtain r where r: "prod_list ps = r * p"
    using remove1_retains_prod by blast
  show "coprime p q"
    apply (rule coprime_prod[of r 1])
    using assms r apply auto
    by blast
qed

(* basically copied from square_free_factorizationD' *)
lemma factorize_polys_square_free_prod_list:
  assumes "factorize_polys ps = (fs,data)"
  shows "square_free (prod_list fs)"
proof (rule square_freeI)
  from factorize_polys_coprime[OF assms]
  have coprime: "p q. p  set fs  q  set fs  p  q  coprime p q" .
  from factorize_polys_square_free[OF assms]
  have sq: "p. p  set fs  square_free p" .
  thus "prod_list fs  0" unfolding prod_list_zero_iff
    using square_free_def by blast
  fix q
  assume "degree q > 0" "q * q dvd prod_list fs"
  from irreducibled_factor[OF this(1)] this(2) obtain q where 
    irr: "irreducible q" and dvd: "q * q dvd prod_list fs" unfolding dvd_def by auto
  hence dvd': "q dvd prod_list fs" unfolding dvd_def by auto
  from irreducible_dvd_prod_list[OF irr dvd'] obtain b where 
    mem: "b  set fs" and dvd1: "q dvd b" by auto
  from dvd1 obtain k where b: "b = q * k" unfolding dvd_def by auto
  from split_list[OF mem] b obtain bs1 bs2 where bs: "fs = bs1 @ b # bs2" by auto
  from irr have q0: "q  0" and dq: "degree q > 0" unfolding irreducibled_def by auto
  have "square_free (q * k)" using sq b mem by auto
  from this[unfolded square_free_def, THEN conjunct2, rule_format, OF dq] 
  have qk: "¬ q dvd k" by simp
  from dvd[unfolded bs b] have "q * q dvd q * (k * prod_list (bs1 @ bs2))"
    by (auto simp: ac_simps)
  with q0 have "q dvd k * prod_list (bs1 @ bs2)" by auto
  with irr qk have "q dvd prod_list (bs1 @ bs2)" by auto
  from irreducible_dvd_prod_list[OF irr this] obtain b' where 
    mem': "b'  set (bs1 @ bs2)" and dvd2: "q dvd b'" by fastforce
  from dvd1 dvd2 have "q dvd gcd b b'" by auto
  with dq is_unit_iff_degree[OF q0] have cop: "¬ coprime b b'" by force
  from mem' have "b'  set fs" unfolding bs by auto
  have b': "b' = b" using coprime
    using b'  set fs cop mem by blast
  with mem' bs factorize_polys_irreducible_monic(1)[OF assms] show False by auto
qed

lemma factorize_polys_map_square_free_prod_list:
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss =  map real_of_rat_poly fs"
  shows "square_free (prod_list fss)"
  using  factorize_polys_square_free_prod_list[OF assms(1)] unfolding assms(2)
  by (simp add: of_rat_hom.square_free_map_poly)

lemma factorize_polys_map_coprime_pderiv:
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss = map real_of_rat_poly fs"
  shows "f. f  set fss  coprime f (pderiv (prod_list fss))"
proof -
  fix f
  assume f: "f  set fss"
  from factorize_polys_map_square_free[OF assms]
  have sq: "p. p  set fss  square_free p" .
  have z: "p. p  set fss  p  0" using sq square_free_def by blast
  have c: "coprime (prod_list fss) (pderiv (prod_list fss))"
    apply (simp add: separable_def[symmetric] square_free_iff_separable[symmetric])
    using factorize_polys_map_square_free_prod_list[OF assms] .
  from coprime_prod_list[OF z c f]
  show "coprime f (pderiv (prod_list fss))" by auto
qed

definition pairwise_coprime_list:: "rat poly list  bool"
  where "pairwise_coprime_list qs = 
    (m < length qs.  n < length qs.
     m  n  coprime (qs ! n) (qs ! m))"

(* Restating factorize_polys_map_coprime to match later definitions *)
lemma coprime_factorize:
  fixes qs:: "rat poly list"
  shows "pairwise_coprime_list (fst(factorize_polys qs))"
proof -
  let ?fs = "fst(factorize_polys qs)"
  have "(m < length ?fs.  n < length ?fs.
     m  n  coprime (?fs ! n) (?fs ! m))"
  proof clarsimp 
    fix m n
    assume "m < length (fst (factorize_polys qs))"
    assume "n < length (fst (factorize_polys qs))"
    assume "m  n"
    show " algebraic_semidom_class.coprime (fst (factorize_polys qs) ! n)
            (fst (factorize_polys qs) ! m)"
      by (metis m < length (fst (factorize_polys qs)) m  n n < length (fst (factorize_polys qs)) coprime_iff_coprime distinct_conv_nth factorize_polys_coprime factorize_polys_def factorize_polys_irreducible_monic(1) fstI nth_mem)
  qed
  then show ?thesis unfolding pairwise_coprime_list_def by auto
qed

lemma squarefree_factorization_degree:
  assumes "square_free_factorization p (c,fs)"
  shows "degree p = sum_list (map (λ(f,c). c * degree f) fs)"
proof -
  have "p =
    Polynomial.smult c
     ((a, i)set fs. a ^ i)" using assms unfolding square_free_factorization_def
    by blast
  then have "degree p = degree ((a, i)set fs. a ^ i)"
    using assms square_free_factorizationD(4) by fastforce
  also have "... = degree (prod_list (map (λ(f,c). f ^ c) fs))"
    by (metis assms prod.distinct_set_conv_list square_free_factorizationD(5))
  also have "... = ((a, i)fs. degree (a ^ i))"
    apply (subst degree_prod_list_eq)
    apply (auto simp add: o_def)
    using assms degree_0 square_free_factorizationD(2) apply blast
    by (simp add: prod.case_distrib)
  ultimately show ?thesis
    by (smt (verit, ccfv_SIG) Polynomial.degree_power_eq Suc_eq_plus1 assms degree_0 map_eq_conv split_cong square_free_factorizationD(2))
qed

subsection "Auxiliary Polynomial Lemmas"
definition roots_of_coprime_r:: "real poly list  real set"
  where "roots_of_coprime_r qs = {x. poly (coprime_r qs) x = 0}"

lemma crb_lem_pos: 
  fixes x:: "real"
  fixes p:: "real poly"
  assumes x: "poly p x = 0" 
  assumes p: "p  0" 
  shows "x < crb p"
  using cauchy_root_bound[of p x] apply (auto)
  unfolding crb_def apply (auto)
  using p x 
  by linarith 

lemma crb_lem_neg: 
  fixes x:: "real"
  fixes p:: "real poly"
  assumes x: "poly p x = 0" 
  assumes p: "p  0" 
  shows "x > -crb p"
  using cauchy_root_bound[of p x] apply (auto)
  unfolding crb_def apply (auto)
  using p x by linarith 

(* Show that the product of the polynomial list is 0 at x iff there is a polynomial 
  in the list that is 0 at x *)
lemma prod_zero:
  shows "x . poly (prod_list (qs:: rat poly list)) x = 0  (q  set (qs). poly q x = 0)" 
  apply auto
  using poly_prod_list_zero_iff apply blast
  using poly_prod_list_zero_iff by blast

lemma coprime_r_zero1: "poly (coprime_r qs) (crb (prod_list qs)) = 0"
  by (simp add: coprime_r_def) 

lemma coprime_r_zero2: "poly (coprime_r qs) (-crb (prod_list qs)) = 0"
  by (simp add: coprime_r_def) 

lemma coprime_mult:
  fixes a:: "real poly"
  fixes b:: "real poly"
  fixes c:: "real poly"
  assumes "algebraic_semidom_class.coprime a b"
  assumes "algebraic_semidom_class.coprime a c"
  shows "algebraic_semidom_class.coprime a (b*c)"
  using assms(1) assms(2) by auto

(* Will be needed when we call the BKR roots on coprime_r *)
lemma coprime_r_coprime_prop:
  fixes ps:: "rat poly list"
  assumes "factorize_polys ps = (fs,data)"
  assumes "fss = map real_of_rat_poly fs"
  shows "f. f  set fss  coprime f (coprime_r fss)"
proof clarsimp 
  fix f:: "real poly"
  assume f_in: "f  set fss"
  have nonz_prod: "prod_list fss  0" using factorize_polys_map_square_free apply (auto)
    using assms(1) assms(2) square_free_def by fastforce 
  have nonz_f: "f  0" using f_in factorize_polys_map_square_free apply (auto)
    using assms(1) assms(2) square_free_def by fastforce 
  have copr_pderiv: "algebraic_semidom_class.coprime f (pderiv (prod_list fss))" using factorize_polys_map_coprime_pderiv
    apply (auto)
    using f_in assms(1) assms(2) by auto 
  have z_iff: "x. poly f x = 0  poly (prod_list fss) x = 0"
    using f_in apply (auto)
    using poly_prod_list_zero_iff by blast 
  let ?inf_p = "[:-(crb (prod_list fss)),1:]::real poly"
  have copr_inf: "algebraic_semidom_class.coprime f ([:-(crb (prod_list fss)),1:])"
  proof - 
    have zero_prop: "x. poly ?inf_p x = 0  x = crb (prod_list fss)" 
      by auto
    have  "poly (prod_list fss) (crb (prod_list fss))  0"
    proof -
      have h: "x. poly (prod_list fss) x = 0  x < (crb (prod_list fss))"
        using nonz_prod crb_lem_pos[where p = "prod_list fss"] 
        by auto
      then  show ?thesis by auto
    qed
    then have nonzero: "poly f (crb (prod_list fss))  0"
      using z_iff by auto
    then have "¬(x. poly f x = 0  poly ?inf_p x = 0)"
      by simp
    have is_unit_gcd: "is_unit (gcd ?inf_p f)"
      using prime_elem_imp_gcd_eq  prime_elem_iff_irreducible linear_irreducible_field
      apply (auto) using nonzero
    proof -
      have f1: "x0. - (x0::real) = - 1 * x0"
        by simp
      have "(1::real)  0"
        by auto
      then have "is_unit (gcd (pCons (- 1 * real_of_int (crb (prod_list fss))) 1) f)"
        using f1 by (metis (no_types) is_unit_gcd nonzero one_poly_eq_simps(1) poly_eq_0_iff_dvd prime_elem_imp_coprime prime_elem_linear_field_poly)
      then show "degree (gcd (pCons (- real_of_int (crb (prod_list fss))) 1) f) = 0"
        by simp
    qed 
    then show ?thesis
      using is_unit_gcd
      by (metis gcd.commute gcd_eq_1_imp_coprime is_unit_gcd_iff) 
  qed
  let ?ninf_p = "[:(crb (prod_list fss)),1:]::real poly"
  have copr_neg_inf: "algebraic_semidom_class.coprime f ([:(crb (prod_list fss)),1:])"
  proof - 
    have h: "x. poly f x = 0  poly (prod_list fss) x = 0"
      using f_in apply (auto)
      using poly_prod_list_zero_iff by blast 
    have zero_prop: "x. poly ?ninf_p x = 0  x = -crb (prod_list fss)" 
      by auto
    have  "poly (prod_list fss) (-crb (prod_list fss))  0"
    proof  -
      have h: "x. poly (prod_list fss) x = 0  x > (-crb (prod_list fss))"
        using nonz_prod crb_lem_neg[where p = "prod_list fss"] 
        by auto
      then  show ?thesis by auto
    qed
    then have nonzero: "poly f (-crb (prod_list fss))  0"
      using z_iff by auto
    then have "¬(x. poly f x = 0  poly ?ninf_p x = 0)"
      using zero_prop by auto
    have is_unit_gcd: "is_unit (gcd ?ninf_p f)"
      using prime_elem_imp_gcd_eq  prime_elem_iff_irreducible linear_irreducible_field
      apply (auto) using nonzero
    proof -
      have f1: "(1::real)  0"
        by auto
      have "¬ pCons (real_of_int (crb (prod_list fss))) 1 dvd f"
        using nonzero by auto
      then show "degree (gcd (pCons (real_of_int (crb (prod_list fss))) 1) f) = 0"
        using f1 by (metis (no_types) Missing_Polynomial_Factorial.is_unit_field_poly coprime_imp_gcd_eq_1 is_unit_gcd_iff one_poly_eq_simps(1) prime_elem_imp_coprime prime_elem_linear_field_poly)
    qed
    then show ?thesis 
      using is_unit_gcd
      by (metis gcd.commute gcd_eq_1_imp_coprime is_unit_gcd_iff) 
  qed
  show "algebraic_semidom_class.coprime f (coprime_r fss)" 
    using copr_pderiv coprime_mult unfolding coprime_r_def
    using copr_inf copr_neg_inf by blast 
qed

lemma coprime_r_nonzero:
  fixes ps:: "rat poly list"
  assumes "factorize_polys ps = (fs,data)"
  assumes nonempty_fs: "fs  []"
  assumes fss_is: "fss = map real_of_rat_poly fs"
  shows "(coprime_r fss)  0"
proof - 
  have nonempty_fss: "fss  []" using nonempty_fs fss_is by auto
  have deg_f: "f  set (fs). degree f > 0"
    using factorize_polys_irreducible_monic 
    apply (auto)
    using assms(1) irreducible_degree_field by blast
  then have deg_fss: "f  set (fss). degree f > 0"
    using fss_is by simp 
  then have fss_nonz: "f  set (fss). f  0"
    by auto
  have "fss  []  ((f  set (fss). (degree f > 0  f  0))  degree (prod_list fss) > 0)"
  proof (induct fss)
    case Nil
    then show ?case
      by blast 
  next
    case (Cons a fss)
    show ?case 
    proof clarsimp 
      assume z_lt: "0 < degree a"
      assume anonz: "a  0"
      assume fnonz: "fset fss. 0 < degree f  f  0"
      have h: "degree (a * prod_list fss) = degree a + degree (prod_list fss) "
        using degree_mult_eq[where p = "a", where q = "prod_list fss"] anonz fnonz 
        by auto
      then show "0 < degree (a * prod_list fss)"
        using z_lt Cons.hyps by auto
    qed
  qed
  then have "degree (prod_list fss) > 0"
    using nonempty_fss deg_fss fss_nonz by auto
  then have pderiv_nonzero: "pderiv (prod_list fss)  0"
    by (simp add: pderiv_eq_0_iff)
  have "(([:-(crb (prod_list fss)),1:]) * ([:(crb (prod_list fss)),1:]))  0"
    by auto
  then show ?thesis using pderiv_nonzero
    unfolding coprime_r_def apply (auto)
    by (metis offset_poly_eq_0_lemma right_minus_eq synthetic_div_unique_lemma) 
qed

lemma Rolle_pderiv:
  fixes q:: "real poly"
  fixes x1 x2:: "real"
  shows "(x1 < x2  poly q x1 = 0  poly q x2 = 0)  (w. x1 < w  w < x2  poly (pderiv q) w = 0)"
  using Rolle_deriv apply (auto)
  by (metis DERIV_unique Rolle continuous_at_imp_continuous_on poly_DERIV poly_differentiable poly_isCont) 

lemma coprime_r_roots_prop: 
  fixes qs:: "real poly list"
  assumes pairwise_rel_prime: "q1 q2. (q1  q2  (List.member qs q1)  (List.member qs q2)) coprime q1 q2"
  shows "x1. x2. ((x1 < x2  (q1  set (qs). (poly q1 x1) = 0)  (q2 set(qs). (poly q2 x2) = 0))  (q. x1 < q  q < x2  poly (coprime_r qs) q = 0))"
proof clarsimp
  fix x1:: "real"
  fix x2:: "real"
  fix q1:: "real poly"
  fix q2:: "real poly"
  assume "x1 < x2"
  assume q1_in: "q1  set qs"
  assume q1_0: "poly q1 x1 = 0"
  assume q2_in: "q2  set qs"
  assume q2_0: "poly q2 x2 = 0"
  have prod_z_x1: "poly (prod_list qs) x1 = 0" using q1_in q1_0
    using poly_prod_list_zero_iff by blast  
  have prod_z_x2: "poly (prod_list qs) x2 = 0" using q2_in q2_0 
    using poly_prod_list_zero_iff by blast 
  have "w>x1. w < x2  poly (pderiv (prod_list qs)) w = 0" 
    using Rolle_pderiv[where q = "prod_list qs"] prod_z_x1 prod_z_x2
    using x1 < x2 by blast
  then obtain w where w_def: "w > x1 w < x2  poly (pderiv (prod_list qs)) w = 0"
    by auto
  then have "poly (coprime_r qs) w = 0"
    unfolding coprime_r_def
    by simp 
  then show "q>x1. q < x2  poly (coprime_r qs) q = 0"
    using w_def by blast 
qed

subsection "Setting Up the Procedure: Lemmas"
definition has_no_zeros::"rat list  bool"
  where "has_no_zeros l = (0  set l)"

lemma hnz_prop: "has_no_zeros l  ¬(k < length l. l ! k = 0)"
  unfolding has_no_zeros_def
  by (simp add: in_set_conv_nth)

definition cast_rat_list:: "rat poly list  real poly list"
  where "cast_rat_list qs = map real_of_rat_poly qs"

definition consistent_sign_vectors_r::"real poly list  real set  rat list set"
  where "consistent_sign_vectors_r qs S = (signs_at qs) ` S"

lemma consistent_sign_vectors_consistent_sign_vectors_r:
  shows"consistent_sign_vectors_r (cast_rat_list qs) S = consistent_sign_vectors qs S"
  unfolding consistent_sign_vectors_r_def cast_rat_list_def consistent_sign_vectors_def
    sign_vec_def signs_at_def
  by auto

(* Relies on coprime_rat_poly_real_poly *)
lemma coprime_over_reals_coprime_over_rats:
  fixes qs:: "rat poly list"
  assumes csa_in: "csa  (consistent_sign_vectors qs UNIV)"
  assumes p1p2: "p1p2  p1 < length csa  p2 < length csa  csa ! p1 = 0  csa ! p2 = 0"
  shows "¬ algebraic_semidom_class.coprime (nth qs p1) (nth qs p2) "
proof - 
  have isx: "(x::real). csa = (sign_vec qs x)" 
    using csa_in unfolding consistent_sign_vectors_def by auto
  then obtain x where havex: "csa = (sign_vec qs x)" by auto
  then have expolys: "poly (real_of_rat_poly (nth qs p1)) x = 0  poly (real_of_rat_poly (nth qs p2)) x = 0"
    using havex unfolding sign_vec_def squash_def
    by (smt (verit) class_field.neg_1_not_0 length_map map_map nth_map one_neq_zero p1p2)
  then have "¬ coprime (real_of_rat_poly (nth qs p1)) ((real_of_rat_poly (nth qs p2))::real poly)"
    using coprime_poly_0 by force
  then show ?thesis
    using coprime_rat_poly_real_poly by auto
qed

(* This and the following lemma are designed to show that if you have two sgas that aren't the same, 
   there's a 0 in between! The proof uses IVT. It hones in on the component that's changed sign 
   (either from 1 to {0, -1} or from -1 to {0, 1}) *)
lemma zero_above: 
  fixes qs:: "rat poly list"
  fixes x1:: "real"
  assumes hnz: "has_no_zeros (sign_vec qs x1)"
  shows "( x2 > x1. ((sign_vec qs x1)  (sign_vec qs x2))  
  ((r::real) > x1. (r  x2  ( q  set(qs). rpoly q r = 0))))"
proof clarsimp 
  fix x2
  assume x1_lt: "x1 < x2"
  assume diff_sign_vec: "sign_vec qs x1  sign_vec qs x2"
  then have "q  set qs. squash (rpoly q x1)  squash (rpoly q x2)"
    unfolding sign_vec_def
    by simp 
  then obtain q where q_prop: "q  set qs  squash (rpoly q x1)  squash (rpoly q x2)"
    by auto
  then have q_in: "q  set qs" by auto
  have poss1: "squash (rpoly q x1) = -1  squash (rpoly q x2) = 1  (r>x1. r  x2  (qset qs. rpoly q r = 0))"
    using poly_IVT_pos[of x1 x2] using x1_lt unfolding squash_def apply (auto)
    using q_prop by fastforce 
  have poss2: "squash (rpoly q x1) = 1  squash (rpoly q x2) = -1  (r>x1. r  x2  (qset qs. rpoly q r = 0))"
    using poly_IVT_neg[of x1 x2] using x1_lt unfolding squash_def apply (auto)
    using q_prop by fastforce 
  have poss3: "squash (rpoly q x2) = 0  (r>x1. r  x2  (qset qs. rpoly q r = 0))"
    using x1_lt unfolding squash_def apply (auto)
    using q_prop by blast
  have "(q  set qs  rpoly q x1 = 0)  ¬has_no_zeros (sign_vec qs x1)"
    unfolding has_no_zeros_def sign_vec_def
    by (smt (verit) image_eqI list.set_map o_apply squash_def)
  have not_poss4: "squash (rpoly q x1)  0" 
    using hnz q_in unfolding squash_def
    using q  set qs  rpoly q x1 = 0  ¬ has_no_zeros (sign_vec qs x1) by auto 
  then show "r>x1. r  x2  (qset qs. rpoly q r = 0)"
    using q_prop poss1 poss2 poss3 not_poss4 by (metis (no_types, lifting) squash_def) 
qed

lemma zero_below: 
  fixes qs:: "rat poly list"
  fixes x1:: "real"
  assumes hnz: "has_no_zeros (sign_vec qs x1)"
  shows "x2 < x1. ((sign_vec qs x1)  (sign_vec qs x2))  
  ((r::real) < x1. (r  x2  ( q  set(qs). rpoly q r = 0)))"
proof clarsimp 
  fix x2
  assume x1_gt: "x2 < x1"
  assume diff_sign_vec: "sign_vec qs x1  sign_vec qs x2"
  then have "q  set qs. squash (rpoly q x1)  squash (rpoly q x2)"
    unfolding sign_vec_def
    by simp 
  then obtain q where q_prop: "q  set qs  squash (rpoly q x1)  squash (rpoly q x2)"
    by auto
  then have q_in: "q  set qs" by auto
  have poss1: "squash (rpoly q x1) = -1  squash (rpoly q x2) = 1  (r<x1. (r  x2  ( q  set(qs). rpoly q r = 0)))"
    using poly_IVT_neg[of x2 x1] using x1_gt unfolding squash_def apply (auto)
    using q_prop by fastforce 
  have poss2: "squash (rpoly q x1) = 1  squash (rpoly q x2) = -1  (r<x1. (r  x2  ( q  set(qs). rpoly q r = 0)))"
    using poly_IVT_pos[of x2 x1] using x1_gt unfolding squash_def apply (auto)
    using q_prop by fastforce 
  have poss3: "squash (rpoly q x2) = 0  (r<x1. (r  x2  ( q  set(qs). rpoly q r = 0)))"
    using x1_gt unfolding squash_def apply (auto)
    using q_prop by blast
  have "(q  set qs  rpoly q x1 = 0)  ¬has_no_zeros (sign_vec qs x1)"
    unfolding has_no_zeros_def sign_vec_def
    by (smt (verit) comp_apply image_eqI list.set_map squash_def)
  have not_poss4: "squash (rpoly q x1)  0" 
    using hnz q_in unfolding squash_def
    using q  set qs  rpoly q x1 = 0  ¬ has_no_zeros (sign_vec qs x1) by auto 
  then show "(r<x1. (r  x2  ( q  set(qs). rpoly q r = 0)))"
    using q_prop poss1 poss2 poss3 not_poss4 by (metis (no_types, lifting) squash_def)
qed

lemma sorted_list_lemma:
  fixes l:: "real list"
  fixes a b:: "real"
  fixes n:: "nat"
  assumes "a < b"
  assumes "(n + 1) < length l"
  assumes strict_sort: "sorted_wrt (<) l"
  assumes lt_a: "(l ! n) < a"
  assumes b_lt: "b < (l ! (n+1))"
  shows "¬((x::real). (List.member l x  a  x  x  b))"
proof -
  have sorted_hyp_var: "q1 < length l. q2 < length l. (q1 < q2 
     (l ! q1) < (l ! q2))" 
    using strict_sort by (auto simp: sorted_wrt_iff_nth_less)
  then have sorted_hyp_var2:  "q1 < length l. q2 < length l. ((l ! q1) < (l ! q2))  q1 < q2"
    using linorder_neqE_nat
    by (metis Groups.add_ac(2) add_mono_thms_linordered_field(5) less_irrefl) 
  have "((x::real). (List.member l x  a  x  x  b))  False"
  proof -
    assume "((x::real). (List.member l x  a  x  x  b))"
    then obtain x where x_prop: "List.member l x  a  x  x  b" by auto
    then have l_prop: "List.member l x  (l ! n) < x  x < (l ! (n+1))"
      using lt_a b_lt by auto
    have nth_l: "l ! n < x" using l_prop by auto
    have np1th_l: "x < l ! (n+1)" using l_prop by auto
    have "k. k < length l  nth l k = x" using l_prop
      by (meson in_set_member index_of_lookup(1) index_of_lookup(2))
    then obtain k where k_prop: "k < length l  nth l k = x" by auto
    have n_lt: "n < k"
      using nth_l sorted_hyp_var k_prop add_lessD1 assms(2) linorder_neqE_nat nat_SN.gt_trans
      by (meson sorted_hyp_var2)
    have k_gt: "k < n + 1" 
      using sorted_hyp_var np1th_l k_prop
      using assms(2) sorted_hyp_var2 by blast 
    show False 
      using n_lt k_gt by auto
  qed
  then show ?thesis by auto
qed

(* This lemma shows that any zero of coprime_r is either between two roots or it's smaller than the 
least root (neg inf) or it's greater than the biggest root (pos inf). *)
lemma roots_of_coprime_r_location_property: 
  fixes qs:: "rat poly list"
  fixes sga:: "rat list"
  fixes zer_list
  assumes pairwise_rel_prime: "pairwise_coprime_list qs"
  assumes all_squarefree: "q. q  set qs  rsquarefree q"
  assumes x1_prop: "sga = sign_vec qs x1"
  assumes hnz: "has_no_zeros sga"
  assumes zer_list_prop: "zer_list = sorted_list_of_set {(x::real). (q  set(qs). (rpoly q x = 0))}"
  shows "zer_list  []  ((x1 < (zer_list ! 0))  (x1 > (zer_list ! (length zer_list - 1)) 
    ( n < (length zer_list - 1). x1 > (zer_list ! n)  x1 < (zer_list ! (n+1)))))"
proof - 
  let ?zer_list = "sorted_list_of_set {(x::real). (q  set(qs). (rpoly q x = 0))} :: real list"
  show ?thesis 
  proof -
    have "((q. (List.member qs q)  q  0)  has_no_zeros (sign_vec qs x1))  ¬ List.member ?zer_list x1"
    proof (induct qs)
      case Nil
      then show ?case  apply (auto)
        by (simp add: member_rec(2)) 
    next
      case (Cons a qs)
      then show ?case 
      proof clarsimp 
        assume imp: "((q. List.member qs q  q  0) 
     has_no_zeros (sign_vec qs x1) 
     ¬ List.member (sorted_list_of_set {x. qset qs. rpoly q x = 0})
         x1)"
        assume nonz: "q. List.member (a # qs) q  q  0"
        assume hnz: " has_no_zeros (sign_vec (a # qs) x1)"
        assume mem_list: "List.member
     (sorted_list_of_set {x. rpoly a x = 0  (qset qs. rpoly q x = 0)})
     x1"
        have "has_no_zeros (sign_vec (a # qs) x1)  has_no_zeros (sign_vec qs x1)"
        proof - 
          assume hnz: "has_no_zeros (sign_vec (a # qs) x1)"
          have same_vec: "sign_vec (a # qs) x1 = ((if rpoly a x1 > 0 then 1 else if rpoly a x1 = 0 then 0 else -1) # sign_vec (qs) x1)"
            unfolding sign_vec_def squash_def by auto
          have "has_no_zeros ((if rpoly a x1 > 0 then 1 else if rpoly a x1 = 0 then 0 else -1) # sign_vec (qs) x1)
              has_no_zeros (sign_vec (qs) x1)"
            by (simp add: has_no_zeros_def)
          then show "has_no_zeros (sign_vec qs x1)" using hnz same_vec by auto
        qed
        then have nmem: "¬ List.member (sorted_list_of_set {x. qset qs. rpoly q x = 0}) x1"
          using hnz nonz imp apply (auto)
          by (simp add: member_rec(1))
        have "q set qs. q  0"
          using nonz using in_set_member apply (auto) by fastforce 
        then have "q set qs. finite {x. rpoly q x = 0}"
          by (simp add: poly_roots_finite)
        then have fin_set: "finite {x. qset qs. rpoly q x = 0}"
          by auto
        have not_in: "x1  {x. qset qs. rpoly q x = 0}" using fin_set nmem set_sorted_list_of_set
            all_squarefree
          apply (auto)
          by (simp add: List.member_def finite {x. qset qs. rpoly q x = 0})
        have x1_in: "x1  {x. rpoly a x = 0  (qset qs. rpoly q x = 0)}"
          using mem_list sorted_list_of_set 
        proof -
          have f1: "r R. ((r::real)  R  ¬ List.member (sorted_list_of_set R) r)  infinite R"
            by (metis in_set_member set_sorted_list_of_set)
          have "finite {r. rpoly a (r::real) = 0}"
            by (metis (full_types) List.finite_set member_rec(1) nonz real_roots_of_rat_poly(1))
          then show ?thesis
            using f1 finite {x. qset qs. rpoly q x = 0} mem_list by fastforce
        qed
        have "rpoly a x1  0" using hnz 
          unfolding has_no_zeros_def sign_vec_def squash_def by auto
        then show "False" using not_in x1_in 
          by auto
      qed
    qed
    then have non_mem: "¬ List.member ?zer_list x1"
      using all_squarefree unfolding rsquarefree_def hnz apply (auto)
      using hnz x1_prop
      by (simp add: in_set_member) 
    have "?zer_list  []  ((x1  (?zer_list ! 0))  (x1  (?zer_list ! (length ?zer_list - 1))))
 ( n < (length ?zer_list - 1). x1 > (?zer_list ! n)  x1 < (?zer_list ! (n+1)))"
    proof - 
      assume nonempty: "?zer_list  []"
      assume x1_asm: "(x1  (?zer_list ! 0))  (x1  (?zer_list ! (length ?zer_list - 1)))"
      have nm1: "x1  ?zer_list ! 0" using non_mem
        using sorted_list_of_set {x. qset qs. rpoly q x = 0}  [] in_set_member
        by (metis (no_types, lifting) in_set_conv_nth length_greater_0_conv)
      have nm2: "x1  ?zer_list ! (length ?zer_list -1)" using non_mem
        by (metis (no_types, lifting) One_nat_def sorted_list_of_set {x. qset qs. rpoly q x = 0}  [] diff_Suc_less in_set_member length_greater_0_conv nth_mem) 
      then have x_asm_var: "x1 > (?zer_list ! 0)  x1 < ?zer_list ! (length ?zer_list -1)"
        using x1_asm nm1 nm2 by auto
      have "(n. (n < (length ?zer_list - 1)   x1  (?zer_list ! n)  x1  (?zer_list ! (n+1))))  False"
      proof - 
        assume assump: "(n. (n < (length ?zer_list - 1)   x1  (?zer_list ! n)  x1  (?zer_list ! (n+1))))"
        have zer_case: "x1  ?zer_list ! 0" using x_asm_var by auto
        have all_n: " n. (n < (length ?zer_list - 1)  x1  ?zer_list ! n) "
        proof -
          fix n
          show n_lt: "(n < (length ?zer_list - 1)  x1  ?zer_list ! n)"
          proof (induct n)
            case 0
            then show ?case using zer_case
              by blast 
          next
            case (Suc n)
            then show ?case 
              using assump
              using Suc_eq_plus1 Suc_lessD by presburger 
          qed
        qed
        have "(length ?zer_list - 2)