Theory Sturm_Rat

section ‹Separation of Roots: Sturm›

text ‹We adapt the existing theory on Sturm's theorem to work on rational numbers instead of real 
  numbers. The reason is that we want to implement real numbers as real algebraic numbers with the
  help of Sturm's theorem to separate the roots. To this end, we just copy the definitions of
  of the algorithms w.r.t. Sturm and let them be executed on rational numbers. We then
  prove that corresponds to a homomorphism and therefore can transfer the existing soundness results.›

theory Sturm_Rat
imports 
  Sturm_Sequences.Sturm_Theorem
  Algebraic_Numbers_Prelim
  Berlekamp_Zassenhaus.Square_Free_Int_To_Square_Free_GFp
begin

hide_const (open) UnivPoly.coeff

(* TODO: Move *)
lemma root_primitive_part [simp]:
  fixes p :: "'a :: {semiring_gcd, semiring_no_zero_divisors} poly"
  shows  "poly (primitive_part p) x = 0 ⟷ poly p x = 0"
proof(cases "p = 0")
  case True
  then show ?thesis by auto
next
  case False
  have "poly p x = content p * poly (primitive_part p) x"
    by (metis content_times_primitive_part poly_smult)
  also have "… = 0 ⟷ poly (primitive_part p) x = 0" by (simp add: False)
  finally show ?thesis by auto
qed

(*TODO: Move*)
lemma irreducible_primitive_part:
  assumes "irreducible p" and "degree p > 0"
  shows "primitive_part p = p"
  using irreducible_content[OF assms(1), unfolded primitive_iff_content_eq_1] assms(2)
  by (auto simp: primitive_part_def abs_poly_def)


subsection ‹Interface for Separating Roots›

text ‹For a given rational polynomial, we need to know how many real roots are in a given closed interval,
  and how many real roots are in an interval $(-\infty,r]$.›

datatype root_info = Root_Info (l_r: "rat ⇒ rat ⇒ nat") (number_root: "rat ⇒ nat")
hide_const (open) l_r
hide_const (open) number_root

definition count_roots_interval_sf :: "real poly ⇒ (real ⇒ real ⇒ nat) × (real ⇒ nat)" where
  "count_roots_interval_sf p = (let ps = sturm_squarefree p
    in ((λ a b. sign_changes ps a - sign_changes ps b + (if poly p a = 0 then 1 else 0)),
       (λ a. sign_changes_neg_inf ps - sign_changes ps a)))"

definition count_roots_interval :: "real poly ⇒ (real ⇒ real ⇒ nat) × (real ⇒ nat)" where
  "count_roots_interval p = (let ps = sturm p
    in ((λ a b. sign_changes ps a - sign_changes ps b + (if poly p a = 0 then 1 else 0)),
       (λ a. sign_changes_neg_inf ps - sign_changes ps a)))"

lemma count_roots_interval_iff: "square_free p ⟹ count_roots_interval p = count_roots_interval_sf p"
  unfolding count_roots_interval_def count_roots_interval_sf_def sturm_squarefree_def
    square_free_iff_separable separable_def by (cases "p = 0", auto)

lemma count_roots_interval_sf: assumes p: "p ≠ 0" 
  and cr: "count_roots_interval_sf p = (cr,nr)"
  shows "a ≤ b ⟹ cr a b = (card {x. a ≤ x ∧ x ≤ b ∧ poly p x = 0})"
    "nr a = card {x. x ≤ a ∧ poly p x = 0}"
proof -
  have id: "a ≤ b ⟹ { x. a ≤ x ∧ x ≤ b ∧ poly p x = 0} = 
    { x. a < x ∧ x ≤ b ∧ poly p x = 0} ∪ (if poly p a = 0 then {a} else {})" 
    (is "_ ⟹ _ = ?R ∪ ?S") using not_less by force
  have RS: "finite ?R" "finite ?S" "?R ∩ ?S = {}"  using p by (auto simp: poly_roots_finite)   
  show "a ≤ b ⟹ cr a b = (card {x. a ≤ x ∧ x ≤ b ∧ poly p x = 0})"
    "nr a = card {x. x ≤ a ∧ poly p x = 0}" using cr unfolding arg_cong[OF id, of card] card_Un_disjoint[OF RS] 
    count_roots_interval_sf_def count_roots_between_correct[symmetric]
    count_roots_below_correct[symmetric] count_roots_below_def
    count_roots_between_def Let_def using p by auto
qed

lemma count_roots_interval: assumes cr: "count_roots_interval p = (cr,nr)"
  and sf: "square_free p"
  shows "a ≤ b ⟹ cr a b = (card {x. a ≤ x ∧ x ≤ b ∧ poly p x = 0})"
    "nr a = card {x. x ≤ a ∧ poly p x = 0}"
  using count_roots_interval_sf[OF _ cr[unfolded count_roots_interval_iff[OF sf]]] 
    sf[unfolded square_free_def] by blast+

definition root_cond :: "int poly × rat × rat ⇒ real ⇒ bool" where
  "root_cond plr x = (case plr of (p,l,r) ⇒ of_rat l ≤ x ∧ x ≤ of_rat r ∧ ipoly p x = 0)"

definition root_info_cond :: "root_info ⇒ int poly ⇒ bool" where
  "root_info_cond ri p ≡ (∀ a b. a ≤ b ⟶ root_info.l_r ri a b = card {x. root_cond (p,a,b) x})
    ∧ (∀ a. root_info.number_root ri a = card {x. x ≤ real_of_rat a ∧ ipoly p x = 0})"

lemma root_info_condD: "root_info_cond ri p ⟹ a ≤ b ⟹ root_info.l_r ri a b = card {x. root_cond (p,a,b) x}"
  "root_info_cond ri p ⟹ root_info.number_root ri a = card {x. x ≤ real_of_rat a ∧ ipoly p x = 0}"
  unfolding root_info_cond_def by auto

  
definition count_roots_interval_sf_rat :: "int poly ⇒ root_info" where
  "count_roots_interval_sf_rat  p = (let pp = real_of_int_poly p;
    (cr,nr) = count_roots_interval_sf pp
  in Root_Info (λ a b. cr (of_rat a) (of_rat b)) (λ a. nr (of_rat a)))"

definition count_roots_interval_rat :: "int poly ⇒ root_info" where
  [code del]: "count_roots_interval_rat  p = (let pp = real_of_int_poly p;
    (cr,nr) = count_roots_interval pp
  in Root_Info (λ a b. cr (of_rat a) (of_rat b)) (λ a. nr (of_rat a)))"

definition count_roots_rat :: "int poly ⇒ nat" where
  [code del]: "count_roots_rat  p = (count_roots (real_of_int_poly p))"

lemma count_roots_interval_sf_rat: assumes p: "p ≠ 0" 
  shows "root_info_cond (count_roots_interval_sf_rat p) p"
proof -  
  let ?p = "real_of_int_poly p"
  let ?r = real_of_rat
  let ?ri = "count_roots_interval_sf_rat p"
  from p have p: "?p ≠ 0" by auto
  obtain cr nr where cr: "count_roots_interval_sf ?p = (cr,nr)" by force
  have "?ri = Root_Info (λa b. cr (?r a) (?r b)) (λa. nr (?r a))"
    unfolding count_roots_interval_sf_rat_def Let_def cr by auto
  hence id: "root_info.l_r ?ri = (λa b. cr (?r a) (?r b))" "root_info.number_root ?ri = (λa. nr (?r a))"
    by auto
  note cr = count_roots_interval_sf[OF p cr]
  show ?thesis unfolding root_info_cond_def id
  proof (intro conjI impI allI)
    fix a
    show "nr (?r a) = card {x. x ≤ (?r a) ∧ ipoly p x = 0}"
      using cr(2)[of "?r a"] by simp
  next
    fix a b :: rat
    assume ab: "a ≤ b"
    from ab have ab: "?r a ≤ ?r b" by (simp add: of_rat_less_eq)
    from cr(1)[OF this] show "cr (?r a) (?r b) = card (Collect (root_cond (p, a, b)))"
      unfolding root_cond_def[abs_def] split by simp
  qed
qed
  
lemma of_rat_of_int_poly: "map_poly of_rat (of_int_poly p) = of_int_poly p" 
  by (subst map_poly_map_poly, auto simp: o_def)
  
lemma square_free_of_int_poly: assumes "square_free p" 
  shows "square_free (of_int_poly p :: 'a :: {field_gcd, field_char_0} poly)" 
proof - 
  have "square_free (map_poly of_rat (of_int_poly p) :: 'a poly)"
    unfolding of_rat_hom.square_free_map_poly by (rule square_free_int_rat[OF assms])
  thus ?thesis unfolding of_rat_of_int_poly .
qed

lemma count_roots_interval_rat: assumes sf: "square_free p"
  shows "root_info_cond (count_roots_interval_rat p) p"
proof -
  from sf have sf: "square_free (real_of_int_poly p)" by (rule square_free_of_int_poly)
  from sf have p: "p ≠ 0" unfolding square_free_def by auto
  show ?thesis
  using count_roots_interval_sf_rat[OF p]
  unfolding count_roots_interval_rat_def count_roots_interval_sf_rat_def 
    Let_def count_roots_interval_iff[OF sf] .
qed


lemma count_roots_rat: "count_roots_rat p = card {x. ipoly p x = (0 :: real)}"
  unfolding count_roots_rat_def count_roots_correct ..

subsection ‹Implementing Sturm on Rational Polynomials›

function sturm_aux_rat where
"sturm_aux_rat (p :: rat poly) q =
    (if degree q = 0 then [p,q] else p # sturm_aux_rat q (-(p mod q)))"
  by (pat_completeness, simp_all)
termination by (relation "measure (degree ∘ snd)",
                simp_all add: o_def degree_mod_less')

lemma sturm_aux_rat: "sturm_aux (real_of_rat_poly p) (real_of_rat_poly q) = 
  map real_of_rat_poly (sturm_aux_rat p q)"
proof (induct p q rule: sturm_aux_rat.induct)
  case (1 p q)
  interpret map_poly_inj_idom_hom of_rat..
  note deg = of_int_hom.degree_map_poly_hom
  show ?case 
    unfolding sturm_aux.simps[of "real_of_rat_poly p"] sturm_aux_rat.simps[of p]
    using 1 by (cases "degree q = 0"; simp add: hom_distribs)
qed

definition sturm_rat where "sturm_rat p = sturm_aux_rat p (pderiv p)"

lemma sturm_rat: "sturm (real_of_rat_poly p) = map real_of_rat_poly (sturm_rat p)"
  unfolding sturm_rat_def sturm_def
  apply (fold of_rat_hom.map_poly_pderiv)
  unfolding sturm_aux_rat..

definition poly_number_rootat :: "rat poly ⇒ rat" where 
  "poly_number_rootat p ≡ sgn (coeff p (degree p))"

definition poly_neg_number_rootat :: "rat poly ⇒ rat" where 
  "poly_neg_number_rootat p ≡ if even (degree p) then sgn (coeff p (degree p))
                                       else -sgn (coeff p (degree p))"

lemma poly_number_rootat: "poly_inf (real_of_rat_poly p) = real_of_rat (poly_number_rootat p)"
  unfolding poly_inf_def poly_number_rootat_def of_int_hom.degree_map_poly_hom of_rat_hom.coeff_map_poly_hom
    real_of_rat_sgn by simp

lemma poly_neg_number_rootat: "poly_neg_inf (real_of_rat_poly p) = real_of_rat (poly_neg_number_rootat p)"
  unfolding poly_neg_inf_def poly_neg_number_rootat_def of_int_hom.degree_map_poly_hom of_rat_hom.coeff_map_poly_hom
    real_of_rat_sgn by (simp add:hom_distribs)

definition sign_changes_rat where
"sign_changes_rat ps (x::rat) =
    length (remdups_adj (filter (λx. x ≠ 0) (map (λp. sgn (poly p x)) ps))) - 1"

definition sign_changes_number_rootat where
  "sign_changes_number_rootat ps = 
    length (remdups_adj (filter (λx. x ≠ 0) (map poly_number_rootat ps))) - 1"

definition sign_changes_neg_number_rootat where
  "sign_changes_neg_number_rootat ps = 
      length (remdups_adj (filter (λx. x ≠ 0) (map poly_neg_number_rootat ps))) - 1"

lemma real_of_rat_list_neq: "list_neq (map real_of_rat xs) 0 
  = map real_of_rat (list_neq xs 0)"
  by (induct xs, auto)

lemma real_of_rat_remdups_adj: "remdups_adj (map real_of_rat xs) = map real_of_rat (remdups_adj xs)"
  by (induct xs rule: remdups_adj.induct, auto)

lemma sign_changes_rat: "sign_changes (map real_of_rat_poly ps) (real_of_rat x)
  = sign_changes_rat ps x" (is "?l = ?r")
proof - 
  define xs where "xs = list_neq (map (λp. sgn (poly p x)) ps) 0"
  have "?l = length (remdups_adj (list_neq (map real_of_rat (map (λxa.  (sgn (poly xa x))) ps)) 0)) - 1"
    by (simp add: sign_changes_def real_of_rat_sgn o_def)
  also have "… = ?r" unfolding sign_changes_rat_def real_of_rat_list_neq 
    unfolding real_of_rat_remdups_adj by simp
  finally show ?thesis .
qed

lemma sign_changes_neg_number_rootat: "sign_changes_neg_inf (map real_of_rat_poly ps)
  =  sign_changes_neg_number_rootat ps" (is "?l = ?r")
proof - 
  have "?l = length (remdups_adj (list_neq (map real_of_rat (map poly_neg_number_rootat ps)) 0)) - 1"
    by (simp add: sign_changes_neg_inf_def o_def real_of_rat_sgn poly_neg_number_rootat)
  also have "… = ?r" unfolding sign_changes_neg_number_rootat_def real_of_rat_list_neq 
    unfolding real_of_rat_remdups_adj by simp
  finally show ?thesis .
qed

lemma sign_changes_number_rootat: "sign_changes_inf (map real_of_rat_poly ps)
  =  sign_changes_number_rootat ps" (is "?l = ?r")
proof - 
  have "?l = length (remdups_adj (list_neq (map real_of_rat (map poly_number_rootat ps)) 0)) - 1"
    unfolding sign_changes_inf_def
    unfolding map_map o_def real_of_rat_sgn poly_number_rootat ..
  also have "… = ?r" unfolding sign_changes_number_rootat_def real_of_rat_list_neq 
    unfolding real_of_rat_remdups_adj by simp
  finally show ?thesis .
qed

lemma count_roots_interval_rat_code[code]:
  "count_roots_interval_rat p = (let rp = map_poly rat_of_int p; ps = sturm_rat rp
    in Root_Info 
      (λ a b. sign_changes_rat ps a - sign_changes_rat ps b + (if poly rp a = 0 then 1 else 0))
      (λ a. sign_changes_neg_number_rootat ps - sign_changes_rat ps a))"
  unfolding count_roots_interval_rat_def Let_def count_roots_interval_def split of_rat_of_int_poly[symmetric, where 'a = real]
    sturm_rat sign_changes_rat 
    by (simp add: sign_changes_neg_number_rootat)

lemma count_roots_rat_code[code]:
  "count_roots_rat p = (let rp = map_poly rat_of_int p in if p = 0 then 0 else let ps = sturm_rat rp
    in sign_changes_neg_number_rootat ps - sign_changes_number_rootat ps)"
  unfolding count_roots_rat_def Let_def sturm_rat count_roots_code of_rat_of_int_poly[symmetric, where 'a = real]
    sign_changes_neg_number_rootat sign_changes_number_rootat
  by simp

hide_const (open) count_roots_interval_sf_rat

text ‹Finally we provide an even more efficient implementation which
  avoids the "poly p x = 0" test, but it is restricted to irreducible polynomials.›

definition root_info :: "int poly ⇒ root_info" where
  "root_info p = (if degree p = 1 then 
    (let x = Rat.Fract (- coeff p 0) (coeff p 1)
     in Root_Info (λ l r. if l ≤ x ∧ x ≤ r then 1 else 0)  (λ b. if x ≤ b then 1 else 0)) else 
    (let rp = map_poly rat_of_int p; ps = sturm_rat rp in 
   Root_Info (λ a b. sign_changes_rat ps a - sign_changes_rat ps b)
      (λ a. sign_changes_neg_number_rootat ps - sign_changes_rat ps a)))" 

lemma root_info:
  assumes irr: "irreducible p" and deg: "degree p > 0"
  shows "root_info_cond (root_info p) p"
proof (cases "degree p = 1")
  case deg: True
  from degree1_coeffs[OF this] obtain a b where p: "p = [:b,a:]" and "a ≠ 0" by auto
  from deg have "degree (real_of_int_poly p) = 1" by simp
  from roots1[OF this, unfolded roots1_def] p
  have id: "(ipoly p x = 0) = ((x :: real) = - b / a)" for x by auto
  have idd: "{x. real_of_rat aa ≤ x ∧
                 x ≤ real_of_rat ba ∧ x = real_of_int (- b) / real_of_int a} 
   = (if real_of_rat aa ≤ real_of_int (- b) / real_of_int a ∧
                 real_of_int (- b) / real_of_int a ≤ real_of_rat ba then {real_of_int (- b) / real_of_int a} else {})" 
    for aa ba by auto
  have iddd: "{x. x ≤ real_of_rat aa ∧ x = real_of_int (- b) / real_of_int a}
    = (if real_of_int (- b) / real_of_int a ≤ real_of_rat aa then {real_of_int (- b) / real_of_int a} else {})" for aa
    by auto
  have id4: "real_of_int x = real_of_rat (rat_of_int x)" for x by simp
  show ?thesis unfolding root_info_def deg unfolding root_info_cond_def id root_cond_def split
    unfolding p Fract_of_int_quotient Let_def idd iddd 
    unfolding id4 of_rat_divide[symmetric] of_rat_less_eq by auto
next
  case False
  have irr_d: "irreducible⇩d p" by (simp add: deg irr irreducible_connect_rev)
  from irreducible⇩d_int_rat[OF this]
  have "irreducible (of_int_poly p :: rat poly)" by auto
  from irreducible_root_free[OF this]
  have idd: "(poly (of_int_poly p) a = 0) = False" for a :: rat
    unfolding root_free_def using False by auto
  have id: "root_info p = count_roots_interval_rat p"
    unfolding root_info_def if_False count_roots_interval_rat_code Let_def idd using False by auto
  show ?thesis unfolding id
    by (rule count_roots_interval_rat[OF irreducible⇩d_square_free[OF irr_d]])
qed

end