Theory Sturm_Multiple_Roots

(*
    Title:      Extension of Sturm's theorem for multiple roots
    Author:     Wenda Li <wl302@cam.ac.uk / liwenda1990@hotmail.com>
*)

section ‹Extension of Sturm's theorem for multiple roots›

theory Sturm_Multiple_Roots 
  imports
    BF_Misc
begin

text ‹The classic Sturm's theorem is used to count real roots WITHOUT multiplicity of a polynomial within 
  an interval. Surprisingly, we can also extend Sturm's theorem to count real roots WITH 
  multiplicity by modifying the signed remainder sequence, which seems to be overlooked by many
  textbooks. 

  Our formal proof is inspired by Theorem 10.5.6 in 
    Rahman, Q.I., Schmeisser, G.: Analytic Theory of Polynomials. Oxford University Press (2002).
›

subsection ‹More results for @{term smods}›

lemma last_smods_gcd:
  fixes p q ::"real poly"
  defines "pp ≡ last (smods p q)" 
  assumes "p≠0"
  shows "pp = smult (lead_coeff pp) (gcd p q)"
  using ‹p≠0› unfolding pp_def
proof (induct "smods p q" arbitrary:p q rule:length_induct)
  case 1
  have ?case when "q=0"
    using that smult_normalize_field_eq ‹p≠0› by auto
  moreover have ?case when "q≠0"
  proof -
    define r where "r= - (p mod q)"
    have smods_cons:"smods p q = p # smods q r"
      unfolding r_def using ‹p≠0› by simp
    have "last (smods q r) = smult (lead_coeff (last (smods q r))) (gcd q r)"
      apply (rule 1(1)[rule_format,of "smods q r" q r])
      using smods_cons ‹q≠0› by auto
    moreover have "gcd p q = gcd q r"
      unfolding r_def by (simp add: gcd.commute that)
    ultimately show ?thesis unfolding smods_cons using ‹q≠0›
      by simp
  qed
  ultimately show ?case by argo
qed

lemma last_smods_nzero:
  assumes "p≠0"
  shows "last (smods p q) ≠0"
  by (metis assms last_in_set no_0_in_smods smods_nil_eq)

subsection ‹Alternative signed remainder sequences›

function smods_ext::"real poly ⇒ real poly ⇒ real poly list" where 
  "smods_ext p q = (if p=0 then [] else
                      (if p mod q ≠ 0  
                        then Cons p (smods_ext q (-(p mod q))) 
                        else Cons p (smods_ext q (pderiv q)))
                   )"
  by auto
termination
  apply (relation "measure (λ(p,q).if p=0 then 0 else if q=0 then 1 else 2+degree q)")
  using degree_mod_less by (auto simp add:degree_pderiv pderiv_eq_0_iff)

lemma smods_ext_prefix:
  fixes p q::"real poly"
  defines "pp ≡ last (smods p q)" 
  assumes "p≠0" "q≠0"
  shows "smods_ext p q = smods p q @ tl (smods_ext pp (pderiv pp))"
  unfolding pp_def using assms(2,3)
proof (induct "smods_ext p q" arbitrary:p q rule:length_induct)
  case 1
  have ?case when "p mod q ≠0"
  proof -
    define pp where "pp=last (smods q (- (p mod q)))"
    have smods_cons:"smods p q = p# smods q (- (p mod q))"
      using ‹p≠0› by auto
    then have pp_last:"pp=last (smods p q)" unfolding pp_def
      by (simp add: "1.prems"(2) pp_def)
    have smods_ext_cons:"smods_ext p q = p # smods_ext q (- (p mod q))"
      using that ‹p≠0› by auto
    have "smods_ext q (- (p mod q)) = smods q (- (p mod q)) @ tl (smods_ext pp (pderiv pp))"
      apply (rule 1(1)[rule_format,of "smods_ext q (- (p mod q))" q "- (p mod q)",folded pp_def])
      using smods_ext_cons ‹q≠0› that by auto
    then show ?thesis unfolding pp_last
      apply (subst smods_cons)
      apply (subst smods_ext_cons)
      by auto
  qed
  moreover have ?case when "p mod q =0" "pderiv q = 0"
  proof -
    have "smods p q = [p,q]"
      using ‹p≠0› ‹q≠0› that by auto
    moreover have "smods_ext p q = [p,q]"
      using that ‹p≠0› by auto
    ultimately show ?case using ‹p≠0› ‹q≠0› that(1) by auto
  qed
  moreover have ?case when "p mod q =0" "pderiv q ≠ 0"
  proof -
    have smods_cons:"smods p q = [p,q]"
      using ‹p≠0› ‹q≠0› that by auto
    have smods_ext_cons:"smods_ext p q = p#smods_ext q (pderiv q)"
      using that ‹p≠0› by auto
    show ?case unfolding smods_cons smods_ext_cons
      apply (simp del:smods_ext.simps)
      by (simp add: "1.prems"(2))
  qed
  ultimately show ?case by argo
qed

lemma no_0_in_smods_ext: "0∉set (smods_ext p q)"
  apply (induct "smods_ext p q" arbitrary:p q)
   apply simp
  by (metis list.distinct(1) list.inject set_ConsD smods_ext.simps)

subsection ‹Sign variations on the alternative signed remainder sequences›

definition changes_itv_smods_ext:: "real ⇒ real ⇒real poly ⇒ real poly ⇒  int" where
  "changes_itv_smods_ext a b p q= (let ps= smods_ext p q in changes_poly_at ps a 
        - changes_poly_at ps b)"

definition changes_gt_smods_ext:: "real ⇒real poly ⇒ real poly ⇒  int" where
  "changes_gt_smods_ext a p q= (let ps= smods_ext p q in changes_poly_at ps a 
        - changes_poly_pos_inf ps)"

definition changes_le_smods_ext:: "real ⇒real poly ⇒ real poly ⇒  int" where
  "changes_le_smods_ext b p q= (let ps= smods_ext p q in changes_poly_neg_inf ps 
        - changes_poly_at ps b)"

definition changes_R_smods_ext:: "real poly ⇒ real poly ⇒  int" where
  "changes_R_smods_ext p q= (let ps= smods_ext p q in changes_poly_neg_inf ps 
        - changes_poly_pos_inf ps)"

subsection ‹Extension of Sturm's theorem for multiple roots›

theorem sturm_ext_interval:
  assumes "a<b" "poly p a≠0" "poly p b≠0"
  shows "proots_count p {x. a<x ∧ x<b} = changes_itv_smods_ext a b p (pderiv p)"
  using assms(2,3)
proof (induct "smods_ext p (pderiv p)" arbitrary:p rule:length_induct)
  case 1
  have "p≠0" using ‹poly p a ≠ 0› by auto 
  have ?case when "pderiv p=0"
  proof -
    obtain c where "p=[:c:]" "c≠0"
      using ‹p≠0› ‹pderiv p = 0› pderiv_iszero by force
    then have "proots_count p {x. a < x ∧ x < b} = 0"
      unfolding proots_count_def by auto
    moreover have "changes_itv_smods_ext a b p (pderiv p) = 0"
      unfolding changes_itv_smods_ext_def using ‹p=[:c:]› ‹c≠0› by auto
    ultimately show ?thesis by auto
  qed
  moreover have ?case when "pderiv p≠0"
  proof -
    define pp where "pp = last (smods p (pderiv p))"
    define lp where "lp = lead_coeff pp"
    define S where "S={x. a < x ∧ x< b}"

    have prefix:"smods_ext p (pderiv p) = smods p (pderiv p) @ tl (smods_ext pp (pderiv pp))"
      using smods_ext_prefix[OF ‹p≠0› ‹pderiv p≠0›,folded pp_def] .
    have pp_gcd:"pp = smult lp (gcd p (pderiv p))"
      using last_smods_gcd[OF ‹p≠0›,of "pderiv p",folded pp_def lp_def] .
    have "pp≠0" "lp≠0" unfolding pp_def lp_def
      subgoal by (rule last_smods_nzero[OF ‹p≠0›])
      subgoal using ‹last (smods p (pderiv p)) ≠ 0› by auto
      done
    have "poly pp a≠0" "poly pp b ≠ 0"
      unfolding pp_gcd using ‹poly p a≠0› ‹poly p b≠0› ‹lp≠0› 
      by (simp_all add:poly_gcd_0_iff)

    have "proots_count pp S = changes_itv_smods_ext a b pp (pderiv pp)" unfolding S_def
    proof (rule 1(1)[rule_format,of "smods_ext pp (pderiv pp)" pp])
      show "length (smods_ext pp (pderiv pp)) < length (smods_ext p (pderiv p))"
        unfolding prefix by (simp add: ‹p ≠ 0› that)
    qed (use ‹poly pp a≠0› ‹poly pp b≠0› in simp_all)
    moreover have "proots_count p S = card (proots_within p S) + proots_count pp S"
    proof -
      have "(∑r∈proots_within p S. order r p) = (∑r∈ proots_within p S. order r pp + 1)"
      proof (rule sum.cong)
        fix x assume "x ∈ proots_within p S"
        have "order x pp = order x (gcd p (pderiv p))"
          unfolding pp_gcd using ‹lp≠0› by (simp add:order_smult)
        also have "... = min (order x p) (order x (pderiv p))"
          apply (subst order_gcd)
          using ‹p≠0› ‹pderiv p≠0› by simp_all
        also have "... = order x (pderiv p)"
          apply (subst order_pderiv)
          using ‹pderiv p≠0› ‹p ≠ 0› ‹x ∈ proots_within p S› order_root by auto
        finally have "order x pp = order x (pderiv p)" .
        moreover have "order x p = order x (pderiv p) + 1"
          apply (subst order_pderiv)
          using ‹pderiv p≠0› ‹p ≠ 0› ‹x ∈ proots_within p S› order_root by auto
        ultimately show "order x p = order x pp + 1" by auto
      qed simp
      also have "... = card (proots_within p S) + (∑r∈ proots_within p S. order r pp)"
        apply (subst sum.distrib)
        by auto
      also have "... = card (proots_within p S) + (∑r∈ proots_within pp S. order r pp)"
      proof -
        have "(∑r∈proots_within p S. order r pp) = (∑r∈proots_within pp S. order r pp)"
          apply (rule sum.mono_neutral_right)
          subgoal using ‹p≠0› by auto
          subgoal unfolding pp_gcd using ‹lp≠0› by (auto simp:poly_gcd_0_iff)
          subgoal unfolding pp_gcd using ‹lp≠0› 
            apply (auto simp:poly_gcd_0_iff order_smult)
            apply (subst order_gcd)
            by (auto simp add: order_root)
          done
        then show ?thesis by simp
      qed
      finally show ?thesis unfolding proots_count_def .
    qed
    moreover have "card (proots_within p S) = changes_itv_smods a b p (pderiv p)" 
      using sturm_interval[OF ‹a<b› ‹poly p a≠0› ‹poly p b≠0›,symmetric] 
      unfolding S_def proots_within_def 
      by (auto intro!:arg_cong[where f=card])
    moreover have "changes_itv_smods_ext a b p (pderiv p) 
            = changes_itv_smods a b p (pderiv p) + changes_itv_smods_ext a b pp (pderiv pp)"
    proof -
      define xs ys where "xs=smods p (pderiv p)" and "ys=smods_ext pp (pderiv pp)"
      have xys: "xs≠[]" "ys≠[]" "last xs=hd ys" "poly (last xs) a≠0" "poly (last xs) b≠0"
        subgoal unfolding xs_def using ‹p≠0› by auto
        subgoal unfolding ys_def using ‹pp≠0› by auto
        subgoal using ‹pp≠0› unfolding xs_def ys_def 
          apply (fold pp_def)
          by auto
        subgoal using ‹poly pp a≠0› unfolding pp_def xs_def .
        subgoal using ‹poly pp b≠0› unfolding pp_def xs_def .
        done
      have "changes_poly_at (xs @ tl ys) a = changes_poly_at xs a + changes_poly_at ys a"
      proof -
        have "changes_poly_at (xs @ tl ys) a  = changes_poly_at (xs @ ys) a"
          unfolding changes_poly_at_def
          apply (simp add:map_tl)
          apply (subst changes_drop_dup[symmetric])
          using that xys by (auto simp add: hd_map last_map)
        also have "... = changes_poly_at xs a + changes_poly_at ys a"
          unfolding changes_poly_at_def
          apply (subst changes_append[symmetric])
          using xys by (auto simp add: hd_map last_map)
        finally show ?thesis .
      qed
      moreover have "changes_poly_at (xs @ tl ys) b = changes_poly_at xs b + changes_poly_at ys b"
      proof -
        have "changes_poly_at (xs @ tl ys) b  = changes_poly_at (xs @ ys) b"
          unfolding changes_poly_at_def
          apply (simp add:map_tl)
          apply (subst changes_drop_dup[symmetric])
          using that xys by (auto simp add: hd_map last_map)
        also have "... = changes_poly_at xs b + changes_poly_at ys b"
          unfolding changes_poly_at_def
          apply (subst changes_append[symmetric])
          using xys by (auto simp add: hd_map last_map)
        finally show ?thesis .
      qed
      ultimately show ?thesis unfolding changes_itv_smods_ext_def changes_itv_smods_def
        apply (fold xs_def ys_def,unfold prefix[folded xs_def ys_def] Let_def)
        by auto
    qed
    ultimately show "proots_count p S = changes_itv_smods_ext a b p (pderiv p)"
      by auto
  qed
  ultimately show ?case by argo
qed

theorem sturm_ext_above:
  assumes "poly p a≠0" 
  shows "proots_count p {x. a<x} = changes_gt_smods_ext a p (pderiv p)"
proof -
  define ps where "ps≡smods_ext p (pderiv p)"
  have "p≠0" and "p∈set ps" using ‹poly p a≠0› ps_def by auto
  obtain ub where ub:"∀p∈set ps. ∀x. poly p x=0 ⟶ x<ub"
    and ub_sgn:"∀x≥ub. ∀p∈set ps. sgn (poly p x) = sgn_pos_inf p"
    and "ub>a"
    using root_list_ub[OF no_0_in_smods_ext,of p "pderiv p",folded ps_def]
    by auto
  have "proots_count p {x. a<x} = proots_count p {x. a<x ∧ x<ub}"
    unfolding proots_count_def
    apply (rule sum.cong)
    by (use ub ‹p∈set ps› in auto)
  moreover have "changes_gt_smods_ext a p (pderiv p) = changes_itv_smods_ext a ub p (pderiv p)"
  proof -
    have "map (sgn ∘ (λp. poly p ub)) ps = map sgn_pos_inf ps"
      using ub_sgn[THEN spec,of ub,simplified] 
      by (metis (mono_tags, lifting) comp_def list.map_cong0)
    hence "changes_poly_at ps ub=changes_poly_pos_inf ps"
      unfolding changes_poly_pos_inf_def changes_poly_at_def
      by (subst changes_map_sgn_eq,metis map_map)
    thus ?thesis unfolding changes_gt_smods_ext_def changes_itv_smods_ext_def ps_def
      by metis
  qed
  moreover have "poly p ub≠0" using ub ‹p∈set ps› by auto
  ultimately show ?thesis using sturm_ext_interval[OF ‹ub>a› assms] by auto
qed

theorem sturm_ext_below:
  assumes "poly p b≠0" 
  shows "proots_count p {x. x<b} = changes_le_smods_ext b p (pderiv p)"
proof -
  define ps where "ps≡smods_ext p (pderiv p)"
  have "p≠0" and "p∈set ps" using ‹poly p b≠0› ps_def by auto
  obtain lb where lb:"∀p∈set ps. ∀x. poly p x=0 ⟶ x>lb"
    and lb_sgn:"∀x≤lb. ∀p∈set ps. sgn (poly p x) = sgn_neg_inf p"
    and "lb<b"
    using root_list_lb[OF no_0_in_smods_ext,of p "pderiv p",folded ps_def] 
    by auto
  have "proots_count p {x. x<b} = proots_count p {x. lb<x ∧ x<b}"
    unfolding proots_count_def by (rule sum.cong,insert lb ‹p∈set ps›,auto)
  moreover have "changes_le_smods_ext b p (pderiv p) = changes_itv_smods_ext lb b p (pderiv p)"
  proof -
    have "map (sgn ∘ (λp. poly p lb)) ps = map sgn_neg_inf ps"
      using lb_sgn[THEN spec,of lb,simplified] 
      by (metis (mono_tags, lifting) comp_def list.map_cong0)
    hence "changes_poly_at ps lb=changes_poly_neg_inf ps"
      unfolding changes_poly_neg_inf_def changes_poly_at_def
      by (subst changes_map_sgn_eq,metis map_map)
    thus ?thesis unfolding changes_le_smods_ext_def changes_itv_smods_ext_def ps_def
      by metis
  qed
  moreover have "poly p lb≠0" using lb ‹p∈set ps› by auto
  ultimately show ?thesis using sturm_ext_interval[OF ‹lb<b› _ assms] by auto
qed

theorem sturm_ext_R: 
  assumes "p≠0"
  shows "proots_count p UNIV = changes_R_smods_ext p (pderiv p)"
proof - 
  define ps where "ps≡smods_ext p (pderiv p)"
  have "p∈set ps" using ps_def ‹p≠0› by auto
  obtain lb where lb:"∀p∈set ps. ∀x. poly p x=0 ⟶ x>lb"
    and lb_sgn:"∀x≤lb. ∀p∈set ps. sgn (poly p x) = sgn_neg_inf p"
    and "lb<0"
    using root_list_lb[OF no_0_in_smods_ext,of p "pderiv p",folded ps_def] 
    by auto
  obtain ub where ub:"∀p∈set ps. ∀x. poly p x=0 ⟶ x<ub"
    and ub_sgn:"∀x≥ub. ∀p∈set ps. sgn (poly p x) = sgn_pos_inf p"
    and "ub>0"
    using root_list_ub[OF no_0_in_smods_ext,of p "pderiv p",folded ps_def] 
    by auto
  have "proots_count p UNIV = proots_count p {x. lb<x ∧ x<ub}"
    unfolding proots_count_def by (rule sum.cong,insert lb ub ‹p∈set ps›,auto)
  moreover have "changes_R_smods_ext p (pderiv p) = changes_itv_smods_ext lb ub p (pderiv p)"
  proof -
    have "map (sgn ∘ (λp. poly p lb)) ps = map sgn_neg_inf ps"
      and "map (sgn ∘ (λp. poly p ub)) ps = map sgn_pos_inf ps"
      using lb_sgn[THEN spec,of lb,simplified] ub_sgn[THEN spec,of ub,simplified] 
      by (metis (mono_tags, lifting) comp_def list.map_cong0)+
    hence "changes_poly_at ps lb=changes_poly_neg_inf ps
          ∧ changes_poly_at ps ub=changes_poly_pos_inf ps"
      unfolding changes_poly_neg_inf_def changes_poly_at_def changes_poly_pos_inf_def
      by (subst (1 3)  changes_map_sgn_eq,metis map_map)
    thus ?thesis unfolding changes_R_smods_ext_def changes_itv_smods_ext_def ps_def
      by metis
  qed
  moreover have "poly p lb≠0" and "poly p ub≠0" using lb ub ‹p∈set ps› by auto
  moreover have "lb<ub" using ‹lb<0› ‹0<ub› by auto
  ultimately show ?thesis using sturm_ext_interval by auto
qed

end