Session Sturm_Sequences

Theory Sturm_Library_Document

section ‹Miscellaneous›
(*<*) theory Sturm_Library_Document imports Main begin end (*>*)

Theory Misc_Polynomial

(* Author: Manuel Eberl <eberlm@in.tum.de> *)
theory Misc_Polynomial
imports "HOL-Computational_Algebra.Polynomial" "HOL-Computational_Algebra.Polynomial_Factorial"
begin

subsection ‹Analysis›

lemma fun_eq_in_ivl:
  assumes "a  b" "x::real. a  x  x  b  eventually (λξ. f ξ = f x) (at x)"
  shows "f a = f b"
proof (rule connected_local_const)
  show "connected {a..b}" "a  {a..b}" "b  {a..b}" using a  b by (auto intro: connected_Icc)
  show "aa{a..b}. eventually (λb. f aa = f b) (at aa within {a..b})"
  proof
    fix x assume "x  {a..b}"
    with assms(2)[rule_format, of x]
    show "eventually (λb. f x = f b) (at x within {a..b})"
      by (auto simp: eventually_at_filter elim: eventually_mono)
  qed
qed

subsection ‹Polynomials›

subsubsection ‹General simplification lemmas›

lemma pderiv_div:
  assumes [simp]: "q dvd p" "q  0"
  shows "pderiv (p div q) = (q * pderiv p - p * pderiv q) div (q * q)"
        "q*q dvd (q * pderiv p - p * pderiv q)"
proof-
  from assms obtain r where "p = q * r" unfolding dvd_def by blast
  hence "q * pderiv p - p * pderiv q = (q * q) * pderiv r"
      by (simp add: algebra_simps pderiv_mult)
  thus "q*q dvd (q * pderiv p - p * pderiv q)" by simp
  note 0 = pderiv_mult[of q "p div q"]
  have 1: "q * (p div q) = p" 
    by (metis assms(1) assms(2) dvd_def nonzero_mult_div_cancel_left)
  have f1: "pderiv (p div q) * (q * q) div (q * q) = pderiv (p div q)"
    by simp
  have f2: "pderiv p = q * pderiv (p div q) + p div q * pderiv q"
    by (metis 0 1) 
  have "p * pderiv q = pderiv q * (q * (p div q))"
    by (metis 1 mult.commute) 
  then have "p * pderiv q = q * (p div q * pderiv q)"
    by fastforce
  then have "q * pderiv p - p * pderiv q = q * (q * pderiv (p div q))"
    using f2 by (metis add_diff_cancel_right' distrib_left)
  then show "pderiv (p div q) = (q * pderiv p - p * pderiv q) div (q * q)"
    using f1 by (metis mult.commute mult.left_commute)
qed


subsubsection ‹Divisibility of polynomials›

text ‹
  Two polynomials that are coprime have no common roots.
›
lemma coprime_imp_no_common_roots:
  "¬ (poly p x = 0  poly q x = 0)" if "coprime p q"
     for x :: "'a :: field"
proof clarify
  assume "poly p x = 0" "poly q x = 0"
  then have "[:-x, 1:] dvd p" "[:-x, 1:] dvd q"
    by (simp_all add: poly_eq_0_iff_dvd)
  with that have "is_unit [:-x, 1:]"
    by (rule coprime_common_divisor)
  then show False
    by (auto simp add: is_unit_pCons_iff)
qed

lemma poly_div:
  assumes "poly q x  0" and "(q::'a :: field poly) dvd p"
  shows "poly (p div q) x = poly p x / poly q x"
proof-
  from assms have [simp]: "q  0" by force
  have "poly q x * poly (p div q) x = poly (q * (p div q)) x" by simp
  also have "q * (p div q) = p"
      using assms by (simp add: div_mult_swap)
  finally show "poly (p div q) x = poly p x / poly q x"
      using assms by (simp add: field_simps)
qed

(* TODO: make this less ugly *)
lemma poly_div_gcd_squarefree_aux:
  assumes "pderiv (p::('a::{field_char_0,field_gcd}) poly)  0"
  defines "d  gcd p (pderiv p)"
  shows "coprime (p div d) (pderiv (p div d))" and
        "x. poly (p div d) x = 0  poly p x = 0"
proof -
  obtain r s where "bezout_coefficients p (pderiv p) = (r, s)"
    by (auto simp add: prod_eq_iff)
  then have "r * p + s * pderiv p = gcd p (pderiv p)"
    by (rule bezout_coefficients)
  then have rs: "d = r * p + s * pderiv p"
    by (simp add: d_def)
  define t where "t = p div d"
  define p' where [simp]: "p' = pderiv p"
  define d' where [simp]: "d' = pderiv d"
  define u where "u = p' div d"
  have A: "p = t * d" and B: "p' = u * d"
    by (simp_all add: t_def u_def d_def algebra_simps)
  from poly_squarefree_decomp[OF assms(1) A B[unfolded p'_def] rs]
    show "x. poly (p div d) x = 0  poly p x = 0" by (auto simp: t_def)

  from rs have C: "s*t*d' = d * (1 - r*t - s*pderiv t)"
      by (simp add: A B algebra_simps pderiv_mult)
  from assms have [simp]: "p  0" "d  0" "t  0"
      by (force, force, subst (asm) A, force)

  have "x. x dvd t; x dvd (pderiv t)  x dvd 1"
  proof -
    fix x assume "x dvd t" "x dvd (pderiv t)"
    then obtain v w where vw:
        "t = x*v" "pderiv t = x*w" unfolding dvd_def by blast
    define x' v' where [simp]: "x' = pderiv x" and [simp]: "v' = pderiv v"
    from vw have "x*v' + v*x' = x*w" by (simp add: pderiv_mult)
    hence "v*x' = x*(w - v')" by (simp add: algebra_simps)
    hence "x dvd v*pderiv x" by simp
    then obtain y where y: "v*x' = x*y" unfolding dvd_def by force
    from t  0 and vw have "x  0" by simp

    have x_pow_n_dvd_d: "n. x^n dvd d"
    proof-
      fix n show "x ^ n dvd d"
      proof (induction n, simp, rename_tac n, case_tac n)
        fix n assume "n = (0::nat)"
        from vw and C have "d = x*(d*r*v + d*s*w + s*v*d')"
            by (simp add: algebra_simps)
        with n = 0 show "x^Suc n dvd d" by (force intro: dvdI)
      next
        fix n n' assume IH: "x^n dvd d" and "n = Suc n'"
        hence [simp]: "Suc n' = n" "x * x^n' = x^n" by simp_all
        define c :: "'a poly" where "c = [:of_nat n:]"
        from pderiv_power_Suc[of x n']
            have [simp]: "pderiv (x^n) = c*x^n' * x'" unfolding c_def
            by simp

        from IH obtain z where d: "d = x^n * z" unfolding dvd_def by blast
        define z' where [simp]: "z' = pderiv z"
        from d d  0 have "x^n  0" "z  0" by force+
        from C d have "x^n*z = z*r*v*x^Suc n + z*s*c*x^n*(v*x') +
                          s*v*z'*x^Suc n + s*z*(v*x')*x^n + s*z*v'*x^Suc n"
            by (simp add: algebra_simps vw pderiv_mult)
        also have "... = x^n*x * (z*r*v + z*s*c*y + s*v*z' + s*z*y + s*z*v')"
            by (simp only: y, simp add: algebra_simps)
        finally have "z = x*(z*r*v+z*s*c*y+s*v*z'+s*z*y+s*z*v')"
             using x^n  0 by force
        hence "x dvd z" by (metis dvd_triv_left)
        with d show "x^Suc n dvd d" by simp
     qed
   qed

   have "degree x = 0"
   proof (cases "degree x", simp)
     case (Suc n)
       hence "x  0" by auto
       with Suc have "degree (x ^ (Suc (degree d))) > degree d"
           by (subst degree_power_eq, simp_all)
       moreover from x_pow_n_dvd_d[of "Suc (degree d)"] and d  0
           have "degree (x^Suc (degree d))  degree d"
                by (simp add: dvd_imp_degree_le)
       ultimately show ?thesis by simp
    qed
    then obtain c where [simp]: "x = [:c:]" by (cases x, simp split: if_split_asm)
    moreover from x  0 have "c  0" by simp
    ultimately show "x dvd 1" using dvdI[of 1 x "[:inverse c:]"]
      by simp
  qed

  then show "coprime t (pderiv t)"
    by (rule coprimeI)
qed

lemma normalize_field:
  "normalize (x :: 'a :: {field,normalization_semidom}) = (if x = 0 then 0 else 1)"
  by (auto simp: is_unit_normalize dvd_field_iff)

lemma normalize_field_eq_1 [simp]:
  "x  0  normalize (x :: 'a :: {field,normalization_semidom}) = 1"
  by (simp add: normalize_field)

lemma unit_factor_field [simp]:
  "unit_factor (x :: 'a :: {field,normalization_semidom}) = x"
  by (cases "x = 0") (auto simp: is_unit_unit_factor dvd_field_iff)

text ‹
  Dividing a polynomial by its gcd with its derivative yields
  a squarefree polynomial with the same roots.
›
lemma poly_div_gcd_squarefree:
  assumes "(p :: ('a::{field_char_0,field_gcd}) poly)  0"
  defines "d  gcd p (pderiv p)"
  shows "coprime (p div d) (pderiv (p div d))" (is ?A) and
        "x. poly (p div d) x = 0  poly p x = 0" (is "x. ?B x")
proof-
  have "?A  (x. ?B x)"
  proof (cases "pderiv p = 0")
    case False
      from poly_div_gcd_squarefree_aux[OF this] show ?thesis
          unfolding d_def by auto
  next
    case True
      then obtain c where [simp]: "p = [:c:]" using pderiv_iszero by blast
      from assms(1) have "c  0" by simp
      from True have "d = smult (inverse c) p"
        by (simp add: d_def normalize_poly_def map_poly_pCons field_simps)
      with p  0 c  0 have "p div d = [:c:]"
        by (simp add: pCons_one)
      with c  0 show ?thesis
        by (simp add: normalize_const_poly is_unit_triv)
  qed
  thus ?A and "x. ?B x" by simp_all
qed



subsubsection ‹Sign changes of a polynomial›

text ‹
  If a polynomial has different signs at two points, it has a root inbetween.
›
lemma poly_different_sign_imp_root:
  assumes "a < b" and "sgn (poly p a)  sgn (poly p (b::real))"
  shows "x. a  x  x  b  poly p x = 0"
proof (cases "poly p a = 0  poly p b = 0")
  case True
    thus ?thesis using assms(1)
        by (elim disjE, rule_tac exI[of _ a], simp,
                        rule_tac exI[of _ b], simp)
next
  case False
    hence [simp]: "poly p a  0" "poly p b  0" by simp_all
    show ?thesis
    proof (cases "poly p a < 0")
      case True
        hence "sgn (poly p a) = -1" by simp
        with assms True have "poly p b > 0"
            by (auto simp: sgn_real_def split: if_split_asm)
        from poly_IVT_pos[OF a < b True this] guess x ..
        thus ?thesis by (intro exI[of _ x], simp)
    next
      case False
        hence "poly p a > 0" by (simp add: not_less less_eq_real_def)
        hence "sgn (poly p a) = 1"  by simp
        with assms False have "poly p b < 0"
            by (auto simp: sgn_real_def not_less
                           less_eq_real_def split: if_split_asm)
        from poly_IVT_neg[OF a < b ‹poly p a > 0 this] guess x ..
        thus ?thesis by (intro exI[of _ x], simp)
    qed
qed

lemma poly_different_sign_imp_root':
  assumes "sgn (poly p a)  sgn (poly p (b::real))"
  shows "x. poly p x = 0"
using assms by (cases "a < b", auto dest!: poly_different_sign_imp_root
                                    simp: less_eq_real_def not_less)


lemma no_roots_inbetween_imp_same_sign:
  assumes "a < b" "x. a  x  x  b  poly p x  (0::real)"
  shows "sgn (poly p a) = sgn (poly p b)"
  using poly_different_sign_imp_root assms by auto



subsubsection ‹Limits of polynomials›

lemma poly_neighbourhood_without_roots:
  assumes "(p :: real poly)  0"
  shows "eventually (λx. poly p x  0) (at x0)"
proof-
  {
    fix ε :: real assume "ε > 0"
    have fin: "finite {x. ¦x-x0¦ < ε  x  x0  poly p x = 0}"
        using poly_roots_finite[OF assms] by simp
    with ε > 0have "δ>0. δε  (x. ¦x-x0¦ < δ  x  x0  poly p x  0)"
    proof (induction "card {x. ¦x-x0¦ < ε  x  x0  poly p x = 0}"
           arbitrary: ε rule: less_induct)
    case (less ε)
    let ?A = "{x. ¦x - x0¦ < ε  x  x0  poly p x = 0}"
    show ?case
      proof (cases "card ?A")
      case 0
        hence "?A = {}" using less by auto
        thus ?thesis using less(2) by (rule_tac exI[of _ ε], auto)
      next
      case (Suc _)
        with less(3) have "{x. ¦x - x0¦ < ε  x  x0  poly p x = 0}  {}" by force
        then obtain x where x_props: "¦x - x0¦ < ε" "x  x0" "poly p x = 0" by blast
        define ε' where "ε' = ¦x - x0¦ / 2"
        have "ε' > 0" "ε' < ε" unfolding ε'_def using x_props by simp_all
        from x_props(1,2) and ε > 0
            have "x  {x'. ¦x' - x0¦ < ε'  x'  x0  poly p x' = 0}" (is "_  ?B")
            by (auto simp: ε'_def)
        moreover from x_props
            have "x  {x. ¦x - x0¦ < ε  x  x0  poly p x = 0}" by blast
        ultimately have "?B  ?A" by auto
        hence "card ?B < card ?A" "finite ?B"
            by (rule psubset_card_mono[OF less(3)],
                blast intro: finite_subset[OF _ less(3)])
        from less(1)[OF this(1) ε' > 0 this(2)]
            show ?thesis using ε' < ε by force
      qed
    qed
  }
  from this[of "1"]
    show ?thesis by (auto simp: eventually_at dist_real_def)
qed


lemma poly_neighbourhood_same_sign:
  assumes "poly p (x0 :: real)  0"
  shows "eventually (λx. sgn (poly p x) = sgn (poly p x0)) (at x0)"
proof -
  have cont: "isCont (λx. sgn (poly p x)) x0"
      by (rule isCont_sgn, rule poly_isCont, rule assms)
  then have "eventually (λx. ¦sgn (poly p x) - sgn (poly p x0)¦ < 1) (at x0)"
      by (auto simp: isCont_def tendsto_iff dist_real_def)
  then show ?thesis
    by (rule eventually_mono) (simp add: sgn_real_def split: if_split_asm)
qed

lemma poly_lhopital:
  assumes "poly p (x::real) = 0" "poly q x = 0" "q  0"
  assumes "(λx. poly (pderiv p) x / poly (pderiv q) x) x y"
  shows "(λx. poly p x / poly q x) x y"
using assms
proof (rule_tac lhopital)
  have "isCont (poly p) x" "isCont (poly q) x" by simp_all
  with assms(1,2) show "poly p x 0" "poly q x 0"
       by (simp_all add: isCont_def)
  from q  0 and ‹poly q x = 0 have "pderiv q  0"
      by (auto dest: pderiv_iszero)
  from poly_neighbourhood_without_roots[OF this]
      show "eventually (λx. poly (pderiv q) x  0) (at x)" .
qed (auto intro: poly_DERIV poly_neighbourhood_without_roots)


lemma poly_roots_bounds:
  assumes "p  0"
  obtains l u
  where "l  (u :: real)"
    and "poly p l  0"
    and "poly p u  0"
    and "{x. x > l  x  u  poly p x = 0 } = {x. poly p x = 0}"
    and "x. x  l  sgn (poly p x) = sgn (poly p l)"
    and "x. x  u  sgn (poly p x) = sgn (poly p u)"
proof
  from assms have "finite {x. poly p x = 0}" (is "finite ?roots")
      using poly_roots_finite by fast
  let ?roots' = "insert 0 ?roots"

  define l where "l = Min ?roots' - 1"
  define u where "u = Max ?roots' + 1"

  from ‹finite ?roots have A: "finite ?roots'"  by auto
  from Min_le[OF this, of 0] and Max_ge[OF this, of 0]
      show "l   u" by (simp add: l_def u_def)
  from Min_le[OF A] have l_props: "x. xl  poly p x  0"
      by (fastforce simp: l_def)
  from Max_ge[OF A] have u_props: "x. xu  poly p x  0"
      by (fastforce simp: u_def)
  from l_props u_props show [simp]: "poly p l  0" "poly p u  0" by auto

  from l_props have "x. poly p x = 0  x > l" by (metis not_le)
  moreover from u_props have "x. poly p x = 0  x  u" by (metis linear)
  ultimately show "{x. x > l  x  u  poly p x = 0} = ?roots" by auto

  {
    fix x assume A: "x < l" "sgn (poly p x)  sgn (poly p l)"
    with poly_IVT_pos[OF A(1), of p] poly_IVT_neg[OF A(1), of p] A(2)
        have False by (auto split: if_split_asm
                         simp: sgn_real_def l_props not_less less_eq_real_def)
  }
  thus "x. x  l  sgn (poly p x) = sgn (poly p l)"
      by (case_tac "x = l", auto simp: less_eq_real_def)

  {
    fix x assume A: "x > u" "sgn (poly p x)  sgn (poly p u)"
    with u_props poly_IVT_neg[OF A(1), of p] poly_IVT_pos[OF A(1), of p] A(2)
        have False by (auto split: if_split_asm
                         simp: sgn_real_def l_props not_less less_eq_real_def)
  }
  thus "x. x  u  sgn (poly p x) = sgn (poly p u)"
      by (case_tac "x = u", auto simp: less_eq_real_def)
qed



definition poly_inf :: "('a::real_normed_vector) poly  'a" where
  "poly_inf p  sgn (coeff p (degree p))"
definition poly_neg_inf :: "('a::real_normed_vector) poly  'a" where
  "poly_neg_inf p  if even (degree p) then sgn (coeff p (degree p))
                                       else -sgn (coeff p (degree p))"
lemma poly_inf_0_iff[simp]:
    "poly_inf p = 0  p = 0" "poly_neg_inf p = 0  p = 0"
    by (auto simp: poly_inf_def poly_neg_inf_def sgn_zero_iff)

lemma poly_inf_mult[simp]:
  fixes p :: "('a::real_normed_field) poly"
  shows "poly_inf (p*q) = poly_inf p * poly_inf q"
        "poly_neg_inf (p*q) = poly_neg_inf p * poly_neg_inf q"
unfolding poly_inf_def poly_neg_inf_def
by ((cases "p = 0  q = 0",auto simp: sgn_zero_iff
         degree_mult_eq[of p q] coeff_mult_degree_sum Real_Vector_Spaces.sgn_mult)[])+


lemma poly_neq_0_at_infinity:
  assumes "(p :: real poly)  0"
  shows "eventually (λx. poly p x  0) at_infinity"
proof-
  from poly_roots_bounds[OF assms] guess l u .
  note lu_props = this
  define b where "b = max (-l) u"
  show ?thesis
  proof (subst eventually_at_infinity, rule exI[of _ b], clarsimp)
    fix x assume A: "¦x¦  b" and B: "poly p x = 0"
    show False
    proof (cases "x  0")
      case True
        with A have "x  u" unfolding b_def by simp
        with lu_props(3, 6) show False by (metis sgn_zero_iff B)
    next
      case False
        with A have "x  l" unfolding b_def by simp
        with lu_props(2, 5) show False by (metis sgn_zero_iff B)
    qed
  qed
qed




lemma poly_limit_aux:
  fixes p :: "real poly"
  defines "n  degree p"
  shows "((λx. poly p x / x ^ n)  coeff p n) at_infinity"
proof (subst filterlim_cong, rule refl, rule refl)
  show "eventually (λx. poly p x / x^n = (in. coeff p i / x^(n-i)))
            at_infinity"
  proof (rule eventually_mono)
    show "eventually (λx::real. x  0) at_infinity"
        by (simp add: eventually_at_infinity, rule exI[of _ 1], auto)
    fix x :: real assume [simp]: "x  0"
    show "poly p x / x ^ n = (in. coeff p i / x ^ (n - i))"
        by (simp add: n_def sum_divide_distrib power_diff poly_altdef)
  qed

  let ?a = "λi. if i = n then coeff p n else 0"
  have "i{..n}. ((λx. coeff p i / x ^ (n - i))  ?a i) at_infinity"
  proof
    fix i assume "i  {..n}"
    hence "i  n" by simp
    show "((λx. coeff p i / x ^ (n - i))  ?a i) at_infinity"
    proof (cases "i = n")
      case True
        thus ?thesis by (intro tendstoI, subst eventually_at_infinity,
                         intro exI[of _ 1], simp add: dist_real_def)
    next
      case False
        hence "n - i > 0" using i  n by simp
        from tendsto_inverse_0 and divide_real_def[of 1]
            have "((λx. 1 / x :: real)  0) at_infinity" by simp
        from tendsto_power[OF this, of "n - i"]
            have "((λx::real. 1 / x ^ (n - i))  0) at_infinity"
                using n - i > 0 by (simp add: power_0_left power_one_over)
        from tendsto_mult_right_zero[OF this, of "coeff p i"]
            have "((λx. coeff p i / x ^ (n - i))  0) at_infinity"
                by (simp add: field_simps)
        thus ?thesis using False by simp
    qed
  qed
  hence "((λx. in. coeff p i / x^(n-i))  (in. ?a i)) at_infinity"
    by (force intro!: tendsto_sum)
  also have "(in. ?a i) = coeff p n" by (subst sum.delta, simp_all)
  finally show "((λx. in. coeff p i / x^(n-i))  coeff p n) at_infinity" .
qed



lemma poly_at_top_at_top:
  fixes p :: "real poly"
  assumes "degree p  1" "coeff p (degree p) > 0"
  shows "LIM x at_top. poly p x :> at_top"
proof-
  let ?n = "degree p"
  define f g where "f x = poly p x / x^?n" and "g x = x ^ ?n" for x :: real

  from poly_limit_aux have "(f  coeff p (degree p)) at_top"
      using tendsto_mono at_top_le_at_infinity unfolding f_def by blast
  moreover from assms
      have "LIM x at_top. g x :> at_top"
        by (auto simp add: g_def intro!: filterlim_pow_at_top filterlim_ident)
  ultimately have "LIM x at_top. f x * g x :> at_top"
      using filterlim_tendsto_pos_mult_at_top assms by simp
  also have "eventually (λx. f x * g x = poly p x) at_top"
      unfolding f_def g_def
      by (subst eventually_at_top_linorder, rule exI[of _ 1],
          simp add: poly_altdef field_simps sum_distrib_left power_diff)
  note filterlim_cong[OF refl refl this]
  finally show ?thesis .
qed

lemma poly_at_bot_at_top:
  fixes p :: "real poly"
  assumes "degree p  1" "coeff p (degree p) < 0"
  shows "LIM x at_top. poly p x :> at_bot"
proof-
  from poly_at_top_at_top[of "-p"] and assms
      have "LIM x at_top. -poly p x :> at_top" by simp
  thus ?thesis by (simp add: filterlim_uminus_at_bot)
qed

lemma poly_lim_inf:
  "eventually (λx::real. sgn (poly p x) = poly_inf p) at_top"
proof (cases "degree p  1")
  case False
    hence "degree p = 0" by simp
    then obtain c where "p = [:c:]" by (cases p, auto split: if_split_asm)
    thus ?thesis
        by (simp add: eventually_at_top_linorder poly_inf_def)
next
  case True
    note deg = this
    let ?lc = "coeff p (degree p)"
    from True have "?lc  0" by force
    show ?thesis
    proof (cases "?lc > 0")
      case True
        from poly_at_top_at_top[OF deg this]
          obtain x0 where "x. x  x0  poly p x  1"
              by (fastforce simp: filterlim_at_top
                      eventually_at_top_linorder less_eq_real_def)
        hence "x. x  x0  sgn (poly p x) = 1" by force
        thus ?thesis by (simp only: eventually_at_top_linorder poly_inf_def,
                             intro exI[of _ x0], simp add: True)
    next
      case False
        hence "?lc < 0" using ?lc  0 by linarith
        from poly_at_bot_at_top[OF deg this]
          obtain x0 where "x. x  x0  poly p x  -1"
              by (fastforce simp: filterlim_at_bot
                      eventually_at_top_linorder less_eq_real_def)
        hence "x. x  x0  sgn (poly p x) = -1" by force
        thus ?thesis by (simp only: eventually_at_top_linorder poly_inf_def,
                             intro exI[of _ x0], simp add: ?lc < 0)
    qed
qed

lemma poly_at_top_or_bot_at_bot:
  fixes p :: "real poly"
  assumes "degree p  1" "coeff p (degree p) > 0"
  shows "LIM x at_bot. poly p x :> (if even (degree p) then at_top else at_bot)"
proof-
  let ?n = "degree p"
  define f g where "f x = poly p x / x ^ ?n" and "g x = x ^ ?n" for x :: real

  from poly_limit_aux have "(f  coeff p (degree p)) at_bot"
      using tendsto_mono at_bot_le_at_infinity by (force simp: f_def [abs_def])
  moreover from assms
      have "LIM x at_bot. g x :> (if even (degree p) then at_top else at_bot)"
        by (auto simp add: g_def split: if_split_asm intro: filterlim_pow_at_bot_even filterlim_pow_at_bot_odd filterlim_ident)
  ultimately have "LIM x at_bot. f x * g x :>
                      (if even ?n then at_top else at_bot)"
      by (auto simp: assms intro: filterlim_tendsto_pos_mult_at_top
                                  filterlim_tendsto_pos_mult_at_bot)
  also have "eventually (λx. f x * g x = poly p x) at_bot"
      unfolding f_def g_def
      by (subst eventually_at_bot_linorder, rule exI[of _ "-1"],
          simp add: poly_altdef field_simps sum_distrib_left power_diff)
  note filterlim_cong[OF refl refl this]
  finally show ?thesis .
qed


lemma poly_at_bot_or_top_at_bot:
  fixes p :: "real poly"
  assumes "degree p  1" "coeff p (degree p) < 0"
  shows "LIM x at_bot. poly p x :> (if even (degree p) then at_bot else at_top)"
proof-
  from poly_at_top_or_bot_at_bot[of "-p"] and assms
      have "LIM x at_bot. -poly p x :>
                (if even (degree p) then at_top else at_bot)" by simp
  thus ?thesis by (auto simp: filterlim_uminus_at_bot)
qed

lemma poly_lim_neg_inf:
  "eventually (λx::real. sgn (poly p x) = poly_neg_inf p) at_bot"
proof (cases "degree p  1")
  case False
    hence "degree p = 0" by simp
    then obtain c where "p = [:c:]" by (cases p, auto split: if_split_asm)
    thus ?thesis
        by (simp add: eventually_at_bot_linorder poly_neg_inf_def)
next
  case True
    note deg = this
    let ?lc = "coeff p (degree p)"
    from True have "?lc  0" by force
    show ?thesis
    proof (cases "?lc > 0")
      case True
        note lc_pos = this
        show ?thesis
        proof (cases "even (degree p)")
          case True
            from poly_at_top_or_bot_at_bot[OF deg lc_pos] and True
              obtain x0 where "x. x  x0  poly p x  1"
                by (fastforce simp add: filterlim_at_top filterlim_at_bot
                        eventually_at_bot_linorder less_eq_real_def)
                hence "x. x  x0  sgn (poly p x) = 1" by force
              thus ?thesis
                by (simp add: True eventually_at_bot_linorder poly_neg_inf_def,
                    intro exI[of _ x0], simp add: lc_pos)
       next
          case False
            from poly_at_top_or_bot_at_bot[OF deg lc_pos] and False
              obtain x0 where "x. x  x0  poly p x  -1"
                by (fastforce simp add: filterlim_at_bot
                        eventually_at_bot_linorder less_eq_real_def)
                hence "x. x  x0  sgn (poly p x) = -1" by force
              thus ?thesis
                by (simp add: False eventually_at_bot_linorder poly_neg_inf_def,
                              intro exI[of _ x0], simp add: lc_pos)
      qed
    next
      case False
        hence lc_neg: "?lc < 0" using ?lc  0 by linarith
        show ?thesis
        proof (cases "even (degree p)")
          case True
            with poly_at_bot_or_top_at_bot[OF deg lc_neg]
              obtain x0 where "x. x  x0  poly p x  -1"
                  by (fastforce simp: filterlim_at_bot
                          eventually_at_bot_linorder less_eq_real_def)
              hence "x. x  x0  sgn (poly p x) = -1" by force
              thus ?thesis
                by (simp only: True eventually_at_bot_linorder poly_neg_inf_def,
                               intro exI[of _ x0], simp add: lc_neg)
        next
          case False
            with poly_at_bot_or_top_at_bot[OF deg lc_neg]
              obtain x0 where "x. x  x0  poly p x  1"
                  by (fastforce simp: filterlim_at_top
                          eventually_at_bot_linorder less_eq_real_def)
              hence "x. x  x0  sgn (poly p x) = 1" by force
              thus ?thesis
                by (simp only: False eventually_at_bot_linorder poly_neg_inf_def,
                               intro exI[of _ x0], simp add: lc_neg)
        qed
    qed
qed





subsubsection ‹Signs of polynomials for sufficiently large values›

lemma polys_inf_sign_thresholds:
  assumes "finite (ps :: real poly set)"
  obtains l u
  where "l  u"
    and "p. p  ps; p  0 
              {x. l < x  x  u  poly p x = 0} = {x. poly p x = 0}"
    and "p x. p  ps; x  u  sgn (poly p x) = poly_inf p"
    and "p x. p  ps; x  l  sgn (poly p x) = poly_neg_inf p"
proof goal_cases
  case prems: 1
  have "l u. l  u  (p x. p  ps  x  u  sgn (poly p x) = poly_inf p) 
              (p x. p  ps  x  l  sgn (poly p x) = poly_neg_inf p)"
      (is "l u. ?P ps l u")
  proof (induction rule: finite_subset_induct[OF assms(1), where A = UNIV])
    case 1
      show ?case by simp
  next
    case 2
      show ?case by (intro exI[of _ 42], simp)
  next
    case prems: (3 p ps)
      from prems(4) obtain l u where lu_props: "?P ps l u" by blast
      from poly_lim_inf obtain u'
          where u'_props: "xu'. sgn (poly p x) = poly_inf p"
          by (force simp add: eventually_at_top_linorder)
      from poly_lim_neg_inf obtain l'
          where l'_props: "xl'. sgn (poly p x) = poly_neg_inf p"
          by (force simp add: eventually_at_bot_linorder)
      show ?case
          by (rule exI[of _ "min l l'"], rule exI[of _ "max u u'"],
              insert lu_props l'_props u'_props, auto)
  qed
  then obtain l u where lu_props: "l  u"
        "p x. p  ps  u  x  sgn (poly p x) = poly_inf p"
        "p x. p  ps  x  l  sgn (poly p x) = poly_neg_inf p" by blast
  moreover {
    fix p x assume A: "p  ps" "p  0" "poly p x = 0"
    from A have "l < x" "x < u"
        by (auto simp: not_le[symmetric] dest: lu_props(2,3))
  }
  note A = this
  have "p. p  ps  p  0 
                 {x. l < x  x  u  poly p x = 0} = {x. poly p x = 0}"
      by (auto dest: A)

  from prems[OF lu_props(1) this lu_props(2,3)] show thesis .
qed


subsubsection ‹Positivity of polynomials›

lemma poly_pos:
  "(x::real. poly p x > 0)  poly_inf p = 1  (x. poly p x  0)"
proof (intro iffI conjI)
  assume A: "x::real. poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. poly p x  0" by simp

  from poly_lim_inf obtain x where "sgn (poly p x) = poly_inf p"
      by (auto simp: eventually_at_top_linorder)
  with A show "poly_inf p = 1"
      by (simp add: sgn_real_def split: if_split_asm)
next
  assume "poly_inf p = 1  (x. poly p x  0)"
  hence A: "poly_inf p = 1" and B: "(x. poly p x  0)" by simp_all
  from poly_lim_inf obtain x where C: "sgn (poly p x) = poly_inf p"
      by (auto simp: eventually_at_top_linorder)
  show "x. poly p x > 0"
  proof (rule ccontr)
    assume "¬(x. poly p x > 0)"
    then obtain x' where "poly p x'  0" by (auto simp: not_less)
    with A and C have "sgn (poly p x')  sgn (poly p x)"
        by (auto simp: sgn_real_def split: if_split_asm)
    from poly_different_sign_imp_root'[OF this] and B
        show False by blast
  qed
qed

lemma poly_pos_greater:
  "(x::real. x > a  poly p x > 0) 
    poly_inf p = 1  (x. x > a  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x::real. x > a  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. x > a  poly p x  0" by auto

  from poly_lim_inf obtain x0 where
      "xx0. sgn (poly p x) = poly_inf p"
      by (auto simp: eventually_at_top_linorder)
  hence "poly_inf p = sgn (poly p (max x0 (a + 1)))" by simp
  also from A have "... = 1" by force
  finally show "poly_inf p = 1" .
next
  assume "poly_inf p = 1  (x. x > a  poly p x  0)"
  hence A: "poly_inf p = 1" and
        B: "(x. x > a  poly p x  0)" by simp_all
  from poly_lim_inf obtain x0 where
      C: "xx0. sgn (poly p x) = poly_inf p"
      by (auto simp: eventually_at_top_linorder)
  hence "sgn (poly p (max x0 (a+1))) = poly_inf p" by simp
  with A have D: "sgn (poly p (max x0 (a+1))) = 1" by simp
  show "x. x > a  poly p x > 0"
  proof (rule ccontr)
    assume "¬(x. x > a  poly p x > 0)"
    then obtain x' where "x' > a" "poly p x'  0" by (auto simp: not_less)
    with A and D have E: "sgn (poly p x')  sgn (poly p (max x0(a+1)))"
        by (auto simp: sgn_real_def split: if_split_asm)
    show False
        apply (cases x' "max x0 (a+1)" rule: linorder_cases)
        using B E x' > a
            apply (force dest!: poly_different_sign_imp_root[of _ _ p])+
        done
  qed
qed

lemma poly_pos_geq:
  "(x::real. x  a  poly p x > 0) 
    poly_inf p = 1  (x. x  a  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x::real. x  a  poly p x > 0"
  hence "x::real. x > a  poly p x > 0" by simp
  also note poly_pos_greater
  finally have "poly_inf p = 1" "(x>a. poly p x  0)" by simp_all
  moreover from A have "poly p a > 0" by simp
  ultimately show "poly_inf p = 1" "xa. poly p x  0"
      by (auto simp: less_eq_real_def)
next
  assume "poly_inf p = 1  (x. x  a  poly p x  0)"
  hence A: "poly_inf p = 1" and
        B: "poly p a  0" and C: "x>a. poly p x  0" by simp_all
  from A and C and poly_pos_greater have "x>a. poly p x > 0" by simp
  moreover with B C poly_IVT_pos[of a "a+1" p] have "poly p a > 0" by force
  ultimately show "xa. poly p x > 0" by (auto simp: less_eq_real_def)
qed

lemma poly_pos_less:
  "(x::real. x < a  poly p x > 0) 
    poly_neg_inf p = 1  (x. x < a  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x::real. x < a  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. x < a  poly p x  0" by auto

  from poly_lim_neg_inf obtain x0 where
      "xx0. sgn (poly p x) = poly_neg_inf p"
      by (auto simp: eventually_at_bot_linorder)
  hence "poly_neg_inf p = sgn (poly p (min x0 (a - 1)))" by simp
  also from A have "... = 1" by force
  finally show "poly_neg_inf p = 1" .
next
  assume "poly_neg_inf p = 1  (x. x < a  poly p x  0)"
  hence A: "poly_neg_inf p = 1" and
        B: "(x. x < a  poly p x  0)" by simp_all
  from poly_lim_neg_inf obtain x0 where
      C: "xx0. sgn (poly p x) = poly_neg_inf p"
      by (auto simp: eventually_at_bot_linorder)
  hence "sgn (poly p (min x0 (a - 1))) = poly_neg_inf p" by simp
  with A have D: "sgn (poly p (min x0 (a - 1))) = 1" by simp
  show "x. x < a  poly p x > 0"
  proof (rule ccontr)
    assume "¬(x. x < a  poly p x > 0)"
    then obtain x' where "x' < a" "poly p x'  0" by (auto simp: not_less)
    with A and D have E: "sgn (poly p x')  sgn (poly p (min x0 (a - 1)))"
        by (auto simp: sgn_real_def split: if_split_asm)
    show False
        apply (cases x' "min x0 (a - 1)" rule: linorder_cases)
        using B E x' < a
            apply (auto dest!: poly_different_sign_imp_root[of _ _ p])+
        done
  qed
qed

lemma poly_pos_leq:
  "(x::real. x  a  poly p x > 0) 
    poly_neg_inf p = 1  (x. x  a  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x::real. x  a  poly p x > 0"
  hence "x::real. x < a  poly p x > 0" by simp
  also note poly_pos_less
  finally have "poly_neg_inf p = 1" "(x<a. poly p x  0)" by simp_all
  moreover from A have "poly p a > 0" by simp
  ultimately show "poly_neg_inf p = 1" "xa. poly p x  0"
      by (auto simp: less_eq_real_def)
next
  assume "poly_neg_inf p = 1  (x. x  a  poly p x  0)"
  hence A: "poly_neg_inf p = 1" and
        B: "poly p a  0" and C: "x<a. poly p x  0" by simp_all
  from A and C and poly_pos_less have "x<a. poly p x > 0" by simp
  moreover with B C poly_IVT_neg[of "a - 1" a p] have "poly p a > 0" by force
  ultimately show "xa. poly p x > 0" by (auto simp: less_eq_real_def)
qed

lemma poly_pos_between_less_less:
  "(x::real. a < x  x < b  poly p x > 0) 
    (a  b  poly p ((a+b)/2) > 0)  (x. a < x  x < b  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x. a < x  x < b  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. a < x  x < b  poly p x  0" by auto
  from A show "a  b  poly p ((a+b)/2) > 0" by (cases "a < b", auto)
next
  assume "(b  a  0 < poly p ((a+b)/2))  (x. a<x  x<b  poly p x  0)"
  hence A: "b  a  0 < poly p ((a+b)/2)" and
        B: "x. a<x  x<b  poly p x  0" by simp_all
  show "x. a < x  x < b  poly p x > 0"
  proof (cases "a  b", simp, clarify, rule_tac ccontr,
         simp only: not_le not_less)
    fix x assume "a < b" "a < x" "x < b" "poly p x  0"
    with B have "poly p x < 0" by (simp add: less_eq_real_def)
    moreover from A and a < b have "poly p ((a+b)/2) > 0" by simp
    ultimately have "sgn (poly p x)  sgn (poly p ((a+b)/2))" by simp
    thus False using B
       apply (cases x "(a+b)/2" rule: linorder_cases)
       apply (drule poly_different_sign_imp_root[of _ _ p], assumption,
              insert a < b a < x x < b , force) []
       apply simp
       apply (drule poly_different_sign_imp_root[of _ _ p], simp,
              insert a < b a < x x < b , force)
       done
  qed
qed

lemma poly_pos_between_less_leq:
  "(x::real. a < x  x  b  poly p x > 0) 
    (a  b  poly p b > 0)  (x. a < x  x  b  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x. a < x  x  b  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. a < x  x  b  poly p x  0" by auto
  from A show "a  b  poly p b > 0" by (cases "a < b", auto)
next
  assume "(b  a  0 < poly p b)  (x. a<x  xb  poly p x  0)"
  hence A: "b  a  0 < poly p b" and B: "x. a<x  xb  poly p x  0"
      by simp_all
  show "x. a < x  x  b  poly p x > 0"
  proof (cases "a  b", simp, clarify, rule_tac ccontr,
         simp only: not_le not_less)
    fix x assume "a < b" "a < x" "x  b" "poly p x  0"
    with B have "poly p x < 0" by (simp add: less_eq_real_def)
    moreover from A and a < b have "poly p b > 0" by simp
    ultimately have "x < b" using x  b by (auto simp: less_eq_real_def)
    from ‹poly p x < 0 and ‹poly p b > 0
        have "sgn (poly p x)  sgn (poly p b)" by simp
    from poly_different_sign_imp_root[OF x < b this] and B and x > a
        show False by auto
  qed
qed

lemma poly_pos_between_leq_less:
  "(x::real. a  x  x < b  poly p x > 0) 
    (a  b  poly p a > 0)  (x. a  x  x < b  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x. a  x  x < b  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. a  x  x < b  poly p x  0" by auto
  from A show "a  b  poly p a > 0" by (cases "a < b", auto)
next
  assume "(b  a  0 < poly p a)  (x. ax  x<b  poly p x  0)"
  hence A: "b  a  0 < poly p a" and B: "x. ax  x<b  poly p x  0"
      by simp_all
  show "x. a  x  x < b  poly p x > 0"
  proof (cases "a  b", simp, clarify, rule_tac ccontr,
         simp only: not_le not_less)
    fix x assume "a < b" "a  x" "x < b" "poly p x  0"
    with B have "poly p x < 0" by (simp add: less_eq_real_def)
    moreover from A and a < b have "poly p a > 0" by simp
    ultimately have "x > a" using x  a by (auto simp: less_eq_real_def)
    from ‹poly p x < 0 and ‹poly p a > 0
        have "sgn (poly p a)  sgn (poly p x)" by simp
    from poly_different_sign_imp_root[OF x > a this] and B and x < b
        show False by auto
  qed
qed

lemma poly_pos_between_leq_leq:
  "(x::real. a  x  x  b  poly p x > 0) 
    (a > b  poly p a > 0)  (x. a  x  x  b  poly p x  0)"
proof (intro iffI conjI)
  assume A: "x. a  x  x  b  poly p x > 0"
  have "x. poly p (x::real) > 0  poly p x  0" by simp
  with A show "x::real. a  x  x  b  poly p x  0" by auto
  from A show "a > b  poly p a > 0" by (cases "a  b", auto)
next
  assume "(b < a  0 < poly p a)  (x. ax  xb  poly p x  0)"
  hence A: "b < a  0 < poly p a" and B: "x. ax  xb  poly p x  0"
      by simp_all
  show "x. a  x  x  b  poly p x > 0"
  proof (cases "a > b", simp, clarify, rule_tac ccontr,
         simp only: not_le not_less)
    fix x assume "a  b" "a  x" "x  b" "poly p x  0"
    with B have "poly p x < 0" by (simp add: less_eq_real_def)
    moreover from A and a  b have "poly p a > 0" by simp
    ultimately have "x > a" using x  a by (auto simp: less_eq_real_def)
    from ‹poly p x < 0 and ‹poly p a > 0
        have "sgn (poly p a)  sgn (poly p x)" by simp
    from poly_different_sign_imp_root[OF x > a this] and B and x  b
        show False by auto
  qed
qed

end

Theory Sturm_Library

section ‹Miscellaneous›
theory Sturm_Library
imports Misc_Polynomial
begin
end

Theory Sturm_Theorem

section ‹Proof of Sturm's Theorem›
(* Author: Manuel Eberl <eberlm@in.tum.de> *)
theory Sturm_Theorem
  imports "HOL-Computational_Algebra.Polynomial"
    "Lib/Sturm_Library" "HOL-Computational_Algebra.Field_as_Ring"
begin

subsection ‹Sign changes of polynomial sequences›

text ‹
  For a given sequence of polynomials, this function computes the number of sign changes
  of the sequence of polynomials evaluated at a given position $x$. A sign change is a
  change from a negative value to a positive one or vice versa; zeros in the sequence are
  ignored.
›

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

text ‹
  The number of sign changes of a sequence distributes over a list in the sense that
  the number of sign changes of a sequence $p_1, \ldots, p_i, \ldots, p_n$ at $x$ is the same
  as the sum of the sign changes of the sequence $p_1, \ldots, p_i$ and $p_i, \ldots, p_n$
  as long as $p_i(x)\neq 0$.
›

lemma sign_changes_distrib:
  "poly p x  0 
      sign_changes (ps1 @ [p] @ ps2) x =
      sign_changes (ps1 @ [p]) x + sign_changes ([p] @ ps2) x"
  by (simp add: sign_changes_def sgn_zero_iff, subst remdups_adj_append, simp)

text ‹
  The following two congruences state that the number of sign changes is the same
  if all the involved signs are the same.
›

lemma sign_changes_cong:
  assumes "length ps = length ps'"
  assumes "i < length ps. sgn (poly (ps!i) x) = sgn (poly (ps'!i) y)"
  shows "sign_changes ps x = sign_changes ps' y"
proof-
 from assms(2) have A: "map (λp. sgn (poly p x)) ps = map (λp. sgn (poly p y)) ps'"
  proof (induction rule: list_induct2[OF assms(1)])
    case 1
      then show ?case by simp
  next
    case (2 p ps p' ps')
      from 2(3)
      have "i<length ps. sgn (poly (ps ! i) x) =
                         sgn (poly (ps' ! i) y)" by auto
      from 2(2)[OF this] 2(3) show ?case by auto
  qed
  show ?thesis unfolding sign_changes_def by (simp add: A)
qed

lemma sign_changes_cong':
  assumes "p  set ps. sgn (poly p x) = sgn (poly p y)"
  shows "sign_changes ps x = sign_changes ps y"
using assms by (intro sign_changes_cong, simp_all)

text ‹
  For a sequence of polynomials of length 3, if the first and the third
  polynomial have opposite and nonzero sign at some $x$, the number of
  sign changes is always 1, irrespective of the sign of the second
  polynomial.
›

lemma sign_changes_sturm_triple:
  assumes "poly p x  0" and "sgn (poly r x) = - sgn (poly p x)"
  shows "sign_changes [p,q,r] x = 1"
unfolding sign_changes_def by (insert assms, auto simp: sgn_real_def)

text ‹
  Finally, we define two additional functions that count the sign changes ``at infinity''.
›

definition sign_changes_inf where
"sign_changes_inf ps =
    length (remdups_adj (filter (λx. x  0) (map poly_inf ps))) - 1"

definition sign_changes_neg_inf where
"sign_changes_neg_inf ps =
    length (remdups_adj (filter (λx. x  0) (map poly_neg_inf ps))) - 1"



subsection ‹Definition of Sturm sequences locale›

text ‹
  We first define the notion of a ``Quasi-Sturm sequence'', which is a weakening of
  a Sturm sequence that captures the properties that are fulfilled by a nonempty
  suffix of a Sturm sequence:
  \begin{itemize}
    \item The sequence is nonempty.
    \item The last polynomial does not change its sign.
    \item If the middle one of three adjacent polynomials has a root at $x$, the other
          two have opposite and nonzero signs at $x$.
  \end{itemize}
›

locale quasi_sturm_seq =
  fixes ps :: "(real poly) list"
  assumes last_ps_sgn_const[simp]:
      "x y. sgn (poly (last ps) x) = sgn (poly (last ps) y)"
  assumes ps_not_Nil[simp]: "ps  []"
  assumes signs: "i x. i < length ps - 2; poly (ps ! (i+1)) x = 0
                      (poly (ps ! (i+2)) x) * (poly (ps ! i) x) < 0"


text ‹
  Now we define a Sturm sequence $p_1,\ldots,p_n$ of a polynomial $p$ in the following way:
  \begin{itemize}
    \item The sequence contains at least two elements.
    \item $p$ is the first polynomial, i.\,e. $p_1 = p$.
    \item At any root $x$ of $p$, $p_2$ and $p$ have opposite sign left of $x$ and
          the same sign right of $x$ in some neighbourhood around $x$.
    \item The first two polynomials in the sequence have no common roots.
    \item If the middle one of three adjacent polynomials has a root at $x$, the other
          two have opposite and nonzero signs at $x$.
  \end{itemize}
›

locale sturm_seq = quasi_sturm_seq +
  fixes p :: "real poly"
  assumes hd_ps_p[simp]: "hd ps = p"
  assumes length_ps_ge_2[simp]: "length ps  2"
  assumes deriv: "x0. poly p x0 = 0 
      eventually (λx. sgn (poly (p * ps!1) x) =
                      (if x > x0 then 1 else -1)) (at x0)"
  assumes p_squarefree: "x. ¬(poly p x = 0  poly (ps!1) x = 0)"
begin

  text ‹
    Any Sturm sequence is obviously a Quasi-Sturm sequence.
›
  lemma quasi_sturm_seq: "quasi_sturm_seq ps" ..

(*<*)
  lemma ps_first_two:
    obtains q ps' where "ps = p # q # ps'"
    using hd_ps_p length_ps_ge_2
      by (cases ps, simp, clarsimp, rename_tac ps', case_tac ps', auto)

  lemma ps_first: "ps ! 0 = p" by (rule ps_first_two, simp)

  lemma [simp]: "p  set ps" using hd_in_set[OF ps_not_Nil] by simp
(*>*)
end

(*<*)
lemma [simp]: "¬quasi_sturm_seq []" by (simp add: quasi_sturm_seq_def)
(*>*)

text ‹
  Any suffix of a Quasi-Sturm sequence is again a Quasi-Sturm sequence.
›

lemma quasi_sturm_seq_Cons:
  assumes "quasi_sturm_seq (p#ps)" and "ps  []"
  shows "quasi_sturm_seq ps"
proof (unfold_locales)
  show "ps  []" by fact
next
  from assms(1) interpret quasi_sturm_seq "p#ps" .
  fix x y
  from last_ps_sgn_const and ps  []
      show "sgn (poly (last ps) x) = sgn (poly (last ps) y)" by simp_all
next
  from assms(1) interpret quasi_sturm_seq "p#ps" .
  fix i x
  assume "i < length ps - 2" and "poly (ps ! (i+1)) x = 0"
  with signs[of "i+1"]
      show "poly (ps ! (i+2)) x * poly (ps ! i) x < 0" by simp
qed



subsection ‹Auxiliary lemmas about roots and sign changes›

lemma sturm_adjacent_root_aux:
  assumes "i < length (ps :: real poly list) - 1"
  assumes "poly (ps ! i) x = 0" and "poly (ps ! (i + 1)) x = 0"
  assumes "i x. i < length ps - 2; poly (ps ! (i+1)) x = 0
                    sgn (poly (ps ! (i+2)) x) = - sgn (poly (ps ! i) x)"
  shows "ji+1. poly (ps ! j) x = 0"
using assms
proof (induction i)
  case 0 thus ?case by (clarsimp, rename_tac j, case_tac j, simp_all)
next
  case (Suc i)
    from Suc.prems(1,2)
        have "sgn (poly (ps ! (i + 2)) x) = - sgn (poly (ps ! i) x)"
        by (intro assms(4)) simp_all
    with Suc.prems(3) have "poly (ps ! i) x = 0" by (simp add: sgn_zero_iff)
    with Suc.prems have "ji+1. poly (ps ! j) x = 0"
        by (intro Suc.IH, simp_all)
    with Suc.prems(3) show ?case
      by (clarsimp, rename_tac j, case_tac "j = Suc (Suc i)", simp_all)
qed


text ‹
  This function splits the sign list of a Sturm sequence at a
  position @{term x} that is not a root of @{term p} into a
  list of sublists such that the number of sign changes within
  every sublist is constant in the neighbourhood of @{term x},
  thus proving that the total number is also constant.
›
fun split_sign_changes where
"split_sign_changes [p] (x :: real) = [[p]]" |
"split_sign_changes [p,q] x = [[p,q]]" |
"split_sign_changes (p#q#r#ps) x =
    (if poly p x  0  poly q x = 0 then
       [p,q,r] # split_sign_changes (r#ps) x
     else
       [p,q] # split_sign_changes (q#r#ps) x)"

lemma (in quasi_sturm_seq) split_sign_changes_subset[dest]:
  "ps'  set (split_sign_changes ps x)  set ps'  set ps"
apply (insert ps_not_Nil)
apply (induction ps x rule: split_sign_changes.induct)
apply (simp, simp, rename_tac p q r ps x,
       case_tac "poly p x  0  poly q x = 0", auto)
done

text ‹
  A custom induction rule for @{term split_sign_changes} that
  uses the fact that all the intermediate parameters in calls
  of @{term split_sign_changes} are quasi-Sturm sequences.
›
lemma (in quasi_sturm_seq) split_sign_changes_induct:
  "p x. P [p] x; p q x. quasi_sturm_seq [p,q]  P [p,q] x;
    p q r ps x. quasi_sturm_seq (p#q#r#ps) 
       poly p x  0  poly q x = 0  P (r#ps) x;
        poly q x  0  P (q#r#ps) x;
        poly p x = 0  P (q#r#ps) x
            P (p#q#r#ps) x  P ps x"
proof goal_cases
  case prems: 1
  have "quasi_sturm_seq ps" ..
  with prems show ?thesis
  proof (induction ps x rule: split_sign_changes.induct)
    case (3 p q r ps x)
      show ?case
      proof (rule 3(5)[OF 3(6)])
        assume A: "poly p x  0" "poly q x = 0"
        from 3(6) have "quasi_sturm_seq (r#ps)"
            by (force dest: quasi_sturm_seq_Cons)
        with 3 A show "P (r # ps) x" by blast
      next
        assume A: "poly q x  0"
        from 3(6) have "quasi_sturm_seq (q#r#ps)"
            by (force dest: quasi_sturm_seq_Cons)
        with 3 A show "P (q # r # ps) x" by blast
      next
        assume A: "poly p x = 0"
        from 3(6) have "quasi_sturm_seq (q#r#ps)"
            by (force dest: quasi_sturm_seq_Cons)
        with 3 A show "P (q # r # ps) x" by blast
      qed
  qed simp_all
qed

text ‹
  The total number of sign changes in the split list is the same
  as the number of sign changes in the original list.
›
lemma (in quasi_sturm_seq) split_sign_changes_correct:
  assumes "poly (hd ps) x0  0"
  defines "sign_changes'  λps x.
               ps'split_sign_changes ps x. sign_changes ps' x"
  shows "sign_changes' ps x0 = sign_changes ps x0"
using assms(1)
proof (induction x0 rule: split_sign_changes_induct)
case (3 p q r ps x0)
  hence "poly p x0  0" by simp
  note IH = 3(2,3,4)
  show ?case
  proof (cases "poly q x0 = 0")
    case True
      from 3 interpret quasi_sturm_seq "p#q#r#ps" by simp
      from signs[of 0] and True have
           sgn_r_x0: "poly r x0 * poly p x0 < 0" by simp
      with 3 have "poly r x0  0" by force
      from sign_changes_distrib[OF this, of "[p,q]" ps]
        have "sign_changes (p#q#r#ps) x0 =
                  sign_changes ([p, q, r]) x0 + sign_changes (r # ps) x0" by simp
      also have "sign_changes (r#ps) x0 = sign_changes' (r#ps) x0"
          using ‹poly q x0 = 0 ‹poly p x0  0 3(5)‹poly r x0  0
          by (intro IH(1)[symmetric], simp_all)
      finally show ?thesis unfolding sign_changes'_def
          using True ‹poly p x0  0 by simp
  next
    case False
      from sign_changes_distrib[OF this, of "[p]" "r#ps"]
          have "sign_changes (p#q#r#ps) x0 =
                  sign_changes ([p,q]) x0 + sign_changes (q#r#ps) x0" by simp
      also have "sign_changes (q#r#ps) x0 = sign_changes' (q#r#ps) x0"
          using ‹poly q x0  0 ‹poly p x0  0 3(5)
          by (intro IH(2)[symmetric], simp_all)
      finally show ?thesis unfolding sign_changes'_def
          using False by simp
    qed
qed (simp_all add: sign_changes_def sign_changes'_def)


text ‹
  We now prove that if $p(x)\neq 0$, the number of sign changes of a Sturm sequence of $p$
  at $x$ is constant in a neighbourhood of $x$.
›

lemma (in quasi_sturm_seq) split_sign_changes_correct_nbh:
  assumes "poly (hd ps) x0  0"
  defines "sign_changes'  λx0 ps x.
               ps'split_sign_changes ps x0. sign_changes ps' x"
  shows "eventually (λx. sign_changes' x0 ps x = sign_changes ps x) (at x0)"
proof (rule eventually_mono)
  show "eventually (λx. p{p  set ps. poly p x0  0}. sgn (poly p x) = sgn (poly p x0)) (at x0)"
      by (rule eventually_ball_finite, auto intro: poly_neighbourhood_same_sign)
next
  fix x
  show "(p{p  set ps. poly p x0  0}. sgn (poly p x) = sgn (poly p x0)) 
        sign_changes' x0 ps x = sign_changes ps x"
  proof -
    fix x assume nbh: "p{p  set ps. poly p x0  0}. sgn (poly p x) = sgn (poly p x0)"
    thus "sign_changes' x0 ps x = sign_changes ps x" using assms(1)
    proof (induction x0 rule: split_sign_changes_induct)
    case (3 p q r ps x0)
      hence "poly p x0  0" by simp
      note IH = 3(2,3,4)
      show ?case
      proof (cases "poly q x0 = 0")
        case True
          from 3 interpret quasi_sturm_seq "p#q#r#ps" by simp
          from signs[of 0] and True have
               sgn_r_x0: "poly r x0 * poly p x0 < 0" by simp
          with 3 have "poly r x0  0" by force
          with nbh 3(5) have "poly r x  0" by (auto simp: sgn_zero_iff)
          from sign_changes_distrib[OF this, of "[p,q]" ps]
            have "sign_changes (p#q#r#ps) x =
                      sign_changes ([p, q, r]) x + sign_changes (r # ps) x" by simp
          also have "sign_changes (r#ps) x = sign_changes' x0 (r#ps) x"
              using ‹poly q x0 = 0 nbh ‹poly p x0  0 3(5)‹poly r x0  0
              by (intro IH(1)[symmetric], simp_all)
          finally show ?thesis unfolding sign_changes'_def
              using True ‹poly p x0  0by simp
      next
        case False
          with nbh 3(5) have "poly q x  0" by (auto simp: sgn_zero_iff)
          from sign_changes_distrib[OF this, of "[p]" "r#ps"]
              have "sign_changes (p#q#r#ps) x =
                      sign_changes ([p,q]) x + sign_changes (q#r#ps) x" by simp
          also have "sign_changes (q#r#ps) x = sign_changes' x0 (q#r#ps) x"
              using ‹poly q x0  0 nbh ‹poly p x0  0 3(5)
              by (intro IH(2)[symmetric], simp_all)
          finally show ?thesis unfolding sign_changes'_def
              using False by simp
        qed
    qed (simp_all add: sign_changes_def sign_changes'_def)
  qed
qed



lemma (in quasi_sturm_seq) hd_nonzero_imp_sign_changes_const_aux:
  assumes "poly (hd ps) x0  0" and "ps'  set (split_sign_changes ps x0)"
  shows "eventually (λx. sign_changes ps' x = sign_changes ps' x0) (at x0)"
using assms
proof (induction x0 rule: split_sign_changes_induct)
  case (1 p x)
    thus ?case by (simp add: sign_changes_def)
next
  case (2 p q x0)
    hence [simp]: "ps' = [p,q]" by simp
    from 2 have "poly p x0  0" by simp
    from 2(1) interpret quasi_sturm_seq "[p,q]" .
    from poly_neighbourhood_same_sign[OF ‹poly p x0  0]
        have "eventually (λx. sgn (poly p x) = sgn (poly p x0)) (at x0)" .
    moreover from last_ps_sgn_const
        have sgn_q: "x. sgn (poly q x) = sgn (poly q x0)" by simp
    ultimately have A:  "eventually (λx. pset[p,q]. sgn (poly p x) =
                           sgn (poly p x0)) (at x0)" by simp
    thus ?case by (force intro: eventually_mono[OF A]
                                sign_changes_cong')
next
  case (3 p q r ps'' x0)
    hence p_not_0: "poly p x0  0" by simp
    note sturm = 3(1)
    note IH = 3(2,3)
    note ps''_props = 3(6)
    show ?case
    proof (cases "poly q x0 = 0")
      case True
        note q_0 = this
        from sturm interpret quasi_sturm_seq "p#q#r#ps''" .
        from signs[of 0] and q_0
            have signs': "poly r x0 * poly p x0 < 0" by simp
        with p_not_0 have r_not_0: "poly r x0  0" by force
        show ?thesis
        proof (cases "ps'  set (split_sign_changes (r # ps'') x0)")
          case True
            show ?thesis by (rule IH(1), fact, fact, simp add: r_not_0, fact)
        next
          case False
            with ps''_props p_not_0 q_0 have ps'_props: "ps' = [p,q,r]" by simp
            from signs[of 0] and q_0
                have sgn_r: "poly r x0 * poly p x0 < 0" by simp
            from p_not_0 sgn_r
              have A: "eventually (λx. sgn (poly p x) = sgn (poly p x0) 
                                     sgn (poly r x) = sgn (poly r x0)) (at x0)"
                  by (intro eventually_conj poly_neighbourhood_same_sign,
                      simp_all add: r_not_0)
            show ?thesis
            proof (rule eventually_mono[OF A], clarify,
                   subst ps'_props, subst sign_changes_sturm_triple)
              fix x assume A: "sgn (poly p x) = sgn (poly p x0)"
                       and B: "sgn (poly r x) = sgn (poly r x0)"
              have prod_neg: "a (b::real). a>0; b>0; a*b<0  False"
                             "a (b::real). a<0; b<0; a*b<0  False"
                  by (drule mult_pos_pos, simp, simp,
                      drule mult_neg_neg, simp, simp)
              from A and ‹poly p x0  0 show "poly p x  0"
                  by (force simp: sgn_zero_iff)

              with sgn_r p_not_0 r_not_0 A B
                  have "poly r x * poly p x < 0" "poly r x  0"
                  by (metis sgn_less sgn_mult, metis sgn_0_0)
              with sgn_r show sgn_r': "sgn (poly r x) = - sgn (poly p x)"
                  apply (simp add: sgn_real_def not_le not_less
                             split: if_split_asm, intro conjI impI)
                  using prod_neg[of "poly r x" "poly p x"] apply force+
                  done

              show "1 = sign_changes ps' x0"
                  by (subst ps'_props, subst sign_changes_sturm_triple,
                      fact, metis A B sgn_r', simp)
            qed
        qed
    next
      case False
        note q_not_0 = this
        show ?thesis
        proof (cases "ps'  set (split_sign_changes (q # r # ps'') x0)")
          case True
            show ?thesis by (rule IH(2), fact, simp add: q_not_0, fact)
        next
          case False
            with ps''_props and q_not_0 have "ps' = [p, q]" by simp
            hence [simp]: "pset ps'. poly p x0  0"
                using q_not_0 p_not_0 by simp
            show ?thesis
            proof (rule eventually_mono)
              fix x assume "pset ps'. sgn (poly p x) = sgn (poly p x0)"
              thus "sign_changes ps' x = sign_changes ps' x0"
                  by (rule sign_changes_cong')
            next
              show "eventually (λx. pset ps'.
                        sgn (poly p x) = sgn (poly p x0)) (at x0)"
                  by (force intro: eventually_ball_finite
                                   poly_neighbourhood_same_sign)
            qed
    qed
  qed
qed


lemma (in quasi_sturm_seq) hd_nonzero_imp_sign_changes_const:
  assumes "poly (hd ps) x0  0"
  shows "eventually (λx. sign_changes ps x = sign_changes ps x0) (at x0)"
proof-
  let ?pss = "split_sign_changes ps x0"
  let ?f = "λpss x. ps'pss. sign_changes ps' x"
  {
    fix pss assume "ps'. ps'set pss 
        eventually (λx. sign_changes ps' x = sign_changes ps' x0) (at x0)"
    hence "eventually (λx. ?f pss x = ?f pss x0) (at x0)"
    proof (induction pss)
      case (Cons ps' pss)
      then show ?case
        apply (rule eventually_mono[OF eventually_conj])
        apply (auto simp add: Cons.prems)
        done
    qed simp
  }
  note A = this[of ?pss]
  have B: "eventually (λx. ?f ?pss x = ?f ?pss x0) (at x0)"
      by (rule A, rule hd_nonzero_imp_sign_changes_const_aux[OF assms], simp)
  note C = split_sign_changes_correct_nbh[OF assms]
  note D = split_sign_changes_correct[OF assms]
  note E = eventually_conj[OF B C]
  show ?thesis by (rule eventually_mono[OF E], auto simp: D)
qed

(*<*)
hide_fact quasi_sturm_seq.split_sign_changes_correct_nbh
hide_fact quasi_sturm_seq.hd_nonzero_imp_sign_changes_const_aux
(*>*)

lemma (in sturm_seq) p_nonzero_imp_sign_changes_const:
  "poly p x0  0 
       eventually (λx. sign_changes ps x = sign_changes ps x0) (at x0)"
  using hd_nonzero_imp_sign_changes_const by simp


text ‹
  If $x$ is a root of $p$ and $p$ is not the zero polynomial, the
  number of sign changes of a Sturm chain of $p$ decreases by 1 at $x$.
›
lemma (in sturm_seq) p_zero:
  assumes "poly p x0 = 0" "p  0"
  shows "eventually (λx. sign_changes ps x =
      sign_changes ps x0 + (if x<x0 then 1 else 0)) (at x0)"
proof-
  from ps_first_two obtain q ps' where [simp]: "ps = p#q#ps'" .
  hence "ps!1 = q" by simp
  have "eventually (λx. x  x0) (at x0)"
      by (simp add: eventually_at, rule exI[of _ 1], simp)
  moreover from p_squarefree and assms(1) have "poly q x0  0" by simp
  {
      have A: "quasi_sturm_seq ps" ..
      with quasi_sturm_seq_Cons[of p "q#ps'"]
          interpret quasi_sturm_seq "q#ps'" by simp
      from ‹poly q x0  0 have "eventually (λx. sign_changes (q#ps') x =
                                     sign_changes (q#ps') x0) (at x0)"
      using hd_nonzero_imp_sign_changes_const[where x0=x0] by simp
  }
  moreover note poly_neighbourhood_without_roots[OF assms(2)] deriv[OF assms(1)]
  ultimately
      have A: "eventually (λx. x  x0  poly p x  0 
                   sgn (poly (p*ps!1) x) = (if x > x0 then 1 else -1) 
                   sign_changes (q#ps') x = sign_changes (q#ps') x0) (at x0)"
           by (simp only: ps!1 = q, intro eventually_conj)
  show ?thesis
  proof (rule eventually_mono[OF A], clarify, goal_cases)
    case prems: (1 x)
    from zero_less_mult_pos have zero_less_mult_pos':
        "a b. (0::real) < a*b; 0 < b  0 < a"
        by (subgoal_tac "a*b = b*a", auto)
    from prems have "poly q x  0" and q_sgn: "sgn (poly q x) =
              (if x < x0 then -sgn (poly p x) else sgn (poly p x))"
        by (auto simp add: sgn_real_def elim: linorder_neqE_linordered_idom
                 dest: mult_neg_neg zero_less_mult_pos
                 zero_less_mult_pos' split: if_split_asm)
     from sign_changes_distrib[OF ‹poly q x  0, of "[p]" ps']
        have "sign_changes ps x = sign_changes [p,q] x + sign_changes (q#ps') x"
            by simp
    also from q_sgn and ‹poly p x  0
        have "sign_changes [p,q] x = (if x<x0 then 1 else 0)"
        by (simp add: sign_changes_def sgn_zero_iff split: if_split_asm)
    also note prems(4)
    also from assms(1) have "sign_changes (q#ps') x0 = sign_changes ps x0"
        by (simp add: sign_changes_def)
    finally show ?case by simp
  qed
qed

text ‹
  With these two results, we can now show that if $p$ is nonzero, the number
  of roots in an interval of the form $(a;b]$ is the difference of the sign changes
  of a Sturm sequence of $p$ at $a$ and $b$.\\
  First, however, we prove the following auxiliary lemma that shows that
  if a function $f: \RR\to\NN$ is locally constant at any $x\in(a;b]$, it is constant
  across the entire interval $(a;b]$:
›

lemma count_roots_between_aux:
  assumes "a  b"
  assumes "x::real. a < x  x  b  eventually (λξ. f ξ = (f x::nat)) (at x)"
  shows "x. a < x  x  b  f x = f b"
proof (clarify)
  fix x assume "x > a" "x  b"
  with assms have "x'. x  x'  x'  b 
                       eventually (λξ. f ξ = f x') (at x')" by auto
  from fun_eq_in_ivl[OF x  b this] show "f x = f b" .
qed

text ‹
  Now we can prove the actual root-counting theorem:
›
theorem (in sturm_seq) count_roots_between:
  assumes [simp]: "p  0" "a  b"
  shows "sign_changes ps a - sign_changes ps b =
             card {x. x > a  x  b  poly p x = 0}"
proof-
  have "sign_changes ps a - int (sign_changes ps b) =
             card {x. x > a  x  b  poly p x = 0}" using a  b
  proof (induction "card {x. x > a  x  b  poly p x = 0}" arbitrary: a b
             rule: less_induct)
    case (less a b)
      show ?case
      proof (cases "x. a < x  x  b  poly p x = 0")
        case False
          hence no_roots: "{x. a < x  x  b  poly p x = 0} = {}"
              (is "?roots=_") by auto
          hence card_roots: "card ?roots = (0::int)" by (subst no_roots, simp)
          show ?thesis
          proof (simp only: card_roots eq_iff_diff_eq_0[symmetric] of_nat_eq_iff,
                 cases "poly p a = 0")
            case False
              with no_roots show "sign_changes ps a = sign_changes ps b"
                  by (force intro: fun_eq_in_ivl a  b
                                   p_nonzero_imp_sign_changes_const)
          next
            case True
              have A: "x. a < x  x  b  sign_changes ps x = sign_changes ps b"
                  apply (rule count_roots_between_aux, fact, clarify)
                  apply (rule p_nonzero_imp_sign_changes_const)
                  apply (insert False, simp)
                  done
              have "eventually (λx. x > a 
                        sign_changes ps x = sign_changes ps a) (at a)"
                  apply (rule eventually_mono [OF p_zero[OF ‹poly p a = 0 p  0]])
                  apply force
                  done
              then obtain δ where δ_props:
                  "δ > 0" "x. x > a  x < a+δ 
                                   sign_changes ps x = sign_changes ps a"
                  by (auto simp: eventually_at dist_real_def)

              show "sign_changes ps a = sign_changes ps b"
              proof (cases "a = b")
                case False
                  define x where "x = min (a+δ/2) b"
                  with False have "a < x" "x < a+δ" "x  b"
                     using δ > 0 a  b by simp_all
                  from δ_props a < x x < a+δ
                      have "sign_changes ps a = sign_changes ps x" by simp
                  also from A a < x x  b have "... = sign_changes ps b"
                      by blast
                  finally show ?thesis .
              qed simp
          qed

      next
        case True
          from poly_roots_finite[OF assms(1)]
            have fin: "finite {x. x > a  x  b  poly p x = 0}"
            by (force intro: finite_subset)
          from True have "{x. x > a  x  b  poly p x = 0}  {}" by blast
          with fin have card_greater_0:
              "card {x. x > a  x  b  poly p x = 0} > 0" by fastforce

          define x2 where "x2 = Min {x. x > a  x  b  poly p x = 0}"
          from Min_in[OF fin] and True
              have x2_props: "x2 > a" "x2  b" "poly p x2 = 0"
              unfolding x2_def by blast+
          from Min_le[OF fin] x2_props
              have x2_le: "x'. x' > a; x'  b; poly p x' = 0  x2  x'"
              unfolding x2_def by simp

          have left: "{x. a < x  x  x2  poly p x = 0} = {x2}"
              using x2_props x2_le by force
          hence [simp]: "card {x. a < x  x  x2  poly p x = 0} = 1" by simp

          from p_zero[OF ‹poly p x2 = 0 p  0,
              unfolded eventually_at dist_real_def] guess ε ..
          hence ε_props: "ε > 0"
              "x. x  x2  ¦x - x2¦ < ε 
                   sign_changes ps x = sign_changes ps x2 +
                       (if x < x2 then 1 else 0)" by auto
          define x1 where "x1 = max (x2 - ε / 2) a"
          have "¦x1 - x2¦ < ε" using ε > 0 x2_props by (simp add: x1_def)
          hence "sign_changes ps x1 =
              (if x1 < x2 then sign_changes ps x2 + 1 else sign_changes ps x2)"
              using ε_props(2) by (cases "x1 = x2", auto)
          hence "sign_changes ps x1 - sign_changes ps x2 = 1"
              unfolding x1_def using x2_props ε > 0 by simp

          also have "x2  {x. a < x  x  x1  poly p x = 0}"
              unfolding x1_def using ε > 0 by force
          with left have "{x. a < x  x  x1  poly p x = 0} = {}" by force
          with less(1)[of a x1] have "sign_changes ps x1 = sign_changes ps a"
              unfolding x1_def ε > 0 by (force simp: card_greater_0)

          finally have signs_left:
              "sign_changes ps a - int (sign_changes ps x2) = 1" by simp

          have "{x. x > a  x  b  poly p x = 0} =
                {x. a < x  x  x2  poly p x = 0} 
                {x. x2 < x  x  b  poly p x = 0}" using x2_props by auto
          also note left
          finally have A: "card {x. x2 < x  x  b  poly p x = 0} + 1 =
              card {x. a < x  x  b  poly p x = 0}" using fin by simp
          hence "card {x. x2 < x  x  b  poly p x = 0} <
                 card {x. a < x  x  b  poly p x = 0}" by simp
          from less(1)[OF this x2_props(2)] and A
              have signs_right: "sign_changes ps x2 - int (sign_changes ps b) + 1 =
                  card {x. a < x  x  b  poly p x = 0}" by simp

          from signs_left and signs_right show ?thesis by simp
        qed
  qed
  thus ?thesis by simp
qed

text ‹
  By applying this result to a sufficiently large upper bound, we can effectively count
  the number of roots ``between $a$ and infinity'', i.\,e. the roots greater than $a$:
›
lemma (in sturm_seq) count_roots_above:
  assumes "p  0"
  shows "sign_changes ps a - sign_changes_inf ps =
             card {x. x > a  poly p x = 0}"
proof-
  have "p  set ps" using hd_in_set[OF ps_not_Nil] by simp
  have "finite (set ps)" by simp
  from polys_inf_sign_thresholds[OF this] guess l u .
  note lu_props = this
  let ?u = "max a u"
  {fix x assume "poly p x = 0" hence "x  ?u"
   using lu_props(3)[OF p  set ps, of x] p  0
       by (cases "u  x", auto simp: sgn_zero_iff)
  } note [simp] = this

  from lu_props
    have "map (λp. sgn (poly p ?u)) ps = map poly_inf ps" by simp
  hence "sign_changes ps a - sign_changes_inf ps =
             sign_changes ps a - sign_changes ps ?u"
      by (simp_all only: sign_changes_def sign_changes_inf_def)
  also from count_roots_between[OF assms] lu_props
      have "... =  card {x. a < x  x  ?u  poly p x = 0}" by simp
  also have "{x. a < x  x  ?u  poly p x = 0} = {x. a < x  poly p x = 0}"
      using lu_props by auto
  finally show ?thesis .
qed

text ‹
  The same works analogously for the number of roots below $a$ and the
  total number of roots.
›

lemma (in sturm_seq) count_roots_below:
  assumes "p  0"
  shows "sign_changes_neg_inf ps - sign_changes ps a =
             card {x. x  a  poly p x = 0}"
proof-
  have "p  set ps" using hd_in_set[OF ps_not_Nil] by simp
  have "finite (set ps)" by simp
  from polys_inf_sign_thresholds[OF this] guess l u .
  note lu_props = this
  let ?l = "min a l"
  {fix x assume "poly p x = 0" hence "x > ?l"
   using lu_props(4)[OF p  set ps, of x] p  0
       by (cases "l < x", auto simp: sgn_zero_iff)
  } note [simp] = this

  from lu_props
    have "map (λp. sgn (poly p ?l)) ps = map poly_neg_inf ps" by simp
  hence "sign_changes_neg_inf ps - sign_changes ps a =
             sign_changes ps ?l - sign_changes ps a"
      by (simp_all only: sign_changes_def sign_changes_neg_inf_def)
  also from count_roots_between[OF assms] lu_props
      have "... =  card {x. ?l < x  x  a  poly p x = 0}" by simp
  also have "{x. ?l < x  x  a  poly p x = 0} = {x. a  x  poly p x = 0}"
      using lu_props by auto
  finally show ?thesis .
qed

lemma (in sturm_seq) count_roots:
  assumes "p  0"
  shows "sign_changes_neg_inf ps - sign_changes_inf ps =
             card {x. poly p x = 0}"
proof-
  have "finite (set ps)" by simp
  from polys_inf_sign_thresholds[OF this] guess l u .
  note lu_props = this

  from lu_props
    have "map (λp. sgn (poly p l)) ps = map poly_neg_inf ps"
         "map (λp. sgn (poly p u)) ps = map poly_inf ps" by simp_all
  hence "sign_changes_neg_inf ps - sign_changes_inf ps =
             sign_changes ps l - sign_changes ps u"
      by (simp_all only: sign_changes_def sign_changes_inf_def
                         sign_changes_neg_inf_def)
  also from count_roots_between[OF assms] lu_props
      have "... =  card {x. l < x  x  u  poly p x = 0}" by simp
  also have "{x. l < x  x  u  poly p x = 0} = {x. poly p x = 0}"
      using lu_props assms by simp
  finally show ?thesis .
qed



subsection ‹Constructing Sturm sequences›

subsection ‹The canonical Sturm sequence›

text ‹
  In this subsection, we will present the canonical Sturm sequence construction for
  a polynomial $p$ without multiple roots that is very similar to the Euclidean
  algorithm:
  $$p_i = \begin{cases}
    p & \text{for}\ i = 1\\
    p' & \text{for}\ i = 2\\
    -p_{i-2}\ \text{mod}\ p_{i-1} & \text{otherwise}
  \end{cases}$$
  We break off the sequence at the first constant polynomial.
›

(*<*)
lemma degree_mod_less': "degree q  0  degree (p mod q) < degree q"
  by (metis degree_0 degree_mod_less not_gr0)
(*>*)

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

(*<*)
declare sturm_aux.simps[simp del]
(*>*)

definition sturm where "sturm p = sturm_aux p (pderiv p)"

text ‹Next, we show some simple facts about this construction:›

lemma sturm_0[simp]: "sturm 0 = [0,0]"
    by (unfold sturm_def, subst sturm_aux.simps, simp)

lemma [simp]: "sturm_aux p q = []  False"
    by (induction p q rule: sturm_aux.induct, subst sturm_aux.simps, auto)

lemma sturm_neq_Nil[simp]: "sturm p  []" unfolding sturm_def by simp

lemma [simp]: "hd (sturm p) = p"
  unfolding sturm_def by (subst sturm_aux.simps, simp)

lemma [simp]: "p  set (sturm p)"
  using hd_in_set[OF sturm_neq_Nil] by simp

lemma [simp]: "length (sturm p)  2"
proof-
  {fix q have "length (sturm_aux p q)  2"
           by (induction p q rule: sturm_aux.induct, subst sturm_aux.simps, auto)
  }
  thus ?thesis unfolding sturm_def .
qed

lemma [simp]: "degree (last (sturm p)) = 0"
proof-
  {fix q have "degree (last (sturm_aux p q)) = 0"
           by (induction p q rule: sturm_aux.induct, subst sturm_aux.simps, simp)
  }
  thus ?thesis unfolding sturm_def .
qed

lemma [simp]: "sturm_aux p q ! 0 = p"
    by (subst sturm_aux.simps, simp)
lemma [simp]: "sturm_aux p q ! Suc 0 = q"
    by (subst sturm_aux.simps, simp)

lemma [simp]: "sturm p ! 0 = p"
    unfolding sturm_def by simp
lemma [simp]: "sturm p ! Suc 0 = pderiv p"
    unfolding sturm_def by simp


lemma sturm_indices:
  assumes "i < length (sturm p) - 2"
  shows "sturm p!(i+2) = -(sturm p!i mod sturm p!(i+1))"
proof-
 {fix ps q
  have "ps = sturm_aux p q; i < length ps - 2
             ps!(i+2) = -(ps!i mod ps!(i+1))"
  proof (induction p q arbitrary: ps i rule: sturm_aux.induct)
    case (1 p q)
      show ?case
      proof (cases "i = 0")
        case False
          then obtain i' where [simp]: "i = Suc i'" by (cases i, simp_all)
          hence "length ps  4" using 1 by simp
          with 1(2) have deg: "degree q  0"
              by (subst (asm) sturm_aux.simps, simp split: if_split_asm)
          with 1(2) obtain ps' where [simp]: "ps = p # ps'"
              by (subst (asm) sturm_aux.simps, simp)
          with 1(2) deg have ps': "ps' = sturm_aux q (-(p mod q))"
              by (subst (asm) sturm_aux.simps, simp)
          from ‹length ps  4 and ps = p # ps' 1(3) False
              have "i - 1 < length ps' - 2" by simp
          from 1(1)[OF deg ps' this]
              show ?thesis by simp
      next
        case True
          with 1(3) have "length ps  3" by simp
          with 1(2) have "degree q  0"
              by (subst (asm) sturm_aux.simps, simp split: if_split_asm)
          with 1(2) have [simp]: "sturm_aux p q ! Suc (Suc 0) = -(p mod q)"
              by (subst sturm_aux.simps, simp)
          from True have "ps!i = p" "ps!(i+1) = q" "ps!(i+2) = -(p mod q)"
              by (simp_all add: 1(2))
          thus ?thesis by simp
      qed
    qed}
  from this[OF sturm_def assms] show ?thesis .
qed

text ‹
  If the Sturm sequence construction is applied to polynomials $p$ and $q$,
  the greatest common divisor of $p$ and $q$ a divisor of every element in the
  sequence. This is obvious from the similarity to Euclid's algorithm for
  computing the GCD.
›

lemma sturm_aux_gcd: "r  set (sturm_aux p q)  gcd p q dvd r"
proof (induction p q rule: sturm_aux.induct)
  case (1 p q)
    show ?case
    proof (cases "r = p")
      case False
        with 1(2) have r: "r  set (sturm_aux q (-(p mod q)))"
          by (subst (asm) sturm_aux.simps, simp split: if_split_asm,
              subst sturm_aux.simps, simp)
        show ?thesis
        proof (cases "degree q = 0")
          case False
            hence "q  0" by force
            with 1(1) [OF False r] show ?thesis
              by (simp add: gcd_mod_right ac_simps)
        next
          case True
            with 1(2) and r  p have "r = q"
                by (subst (asm) sturm_aux.simps, simp)
            thus ?thesis by simp
        qed
    qed simp
qed

lemma sturm_gcd: "r  set (sturm p)  gcd p (pderiv p) dvd r"
    unfolding sturm_def by (rule sturm_aux_gcd)

text ‹
  If two adjacent polynomials in the result of the canonical Sturm chain construction
  both have a root at some $x$, this $x$ is a root of all polynomials in the sequence.
›

lemma sturm_adjacent_root_propagate_left:
  assumes "i < length (sturm (p :: real poly)) - 1"
  assumes "poly (sturm p ! i) x = 0"
      and "poly (sturm p ! (i + 1)) x = 0"
  shows "ji+1. poly (sturm p ! j) x = 0"
using assms(2)
proof (intro sturm_adjacent_root_aux[OF assms(1,2,3)], goal_cases)
  case prems: (1 i x)
    let ?p = "sturm p ! i"
    let ?q = "sturm p ! (i + 1)"
    let ?r = "sturm p ! (i + 2)"
    from sturm_indices[OF prems(2)] have "?p = ?p div ?q * ?q - ?r"
        by (simp add: div_mult_mod_eq)
    hence "poly ?p x = poly (?p div ?q * ?q - ?r) x" by simp
    hence "poly ?p x = -poly ?r x" using prems(3) by simp
    thus ?case by (simp add: sgn_minus)
qed

text ‹
  Consequently, if this is the case in the canonical Sturm chain of $p$,
  $p$ must have multiple roots.
›
lemma sturm_adjacent_root_not_squarefree:
  assumes "i < length (sturm (p :: real poly)) - 1"
          "poly (sturm p ! i) x = 0" "poly (sturm p ! (i + 1)) x = 0"
  shows "¬rsquarefree p"
proof-
  from sturm_adjacent_root_propagate_left[OF assms]
      have "poly p x = 0" "poly (pderiv p) x = 0" by auto
  thus ?thesis by (auto simp: rsquarefree_roots)
qed


text ‹
  Since the second element of the sequence is chosen to be the derivative of $p$,
  $p_1$ and $p_2$ fulfil the property demanded by the definition of a Sturm sequence
  that they locally have opposite sign left of a root $x$ of $p$ and the same sign
  to the right of $x$.
›

lemma sturm_firsttwo_signs_aux:
  assumes "(p :: real poly)  0" "q  0"
  assumes q_pderiv:
      "eventually (λx. sgn (poly q x) = sgn (poly (pderiv p) x)) (at x0)"
  assumes p_0: "poly p (x0::real) = 0"
  shows "eventually (λx. sgn (poly (p*q) x) = (if x > x0 then 1 else -1)) (at x0)"
proof-
  have A: "eventually (λx. poly p x  0  poly q x  0 
               sgn (poly q x) = sgn (poly (pderiv p) x)) (at x0)"
      using p  0  q  0
      by (intro poly_neighbourhood_same_sign q_pderiv
                poly_neighbourhood_without_roots eventually_conj)
  then obtain ε where ε_props: "ε > 0" "x. x  x0  ¦x - x0¦ < ε 
      poly p x  0  poly q x  0  sgn (poly (pderiv p) x) = sgn (poly q x)"
      by (auto simp: eventually_at dist_real_def)
  have sqr_pos: "x::real. x  0  sgn x * sgn x = 1"
      by (auto simp: sgn_real_def)

  show ?thesis
  proof (simp only: eventually_at dist_real_def, rule exI[of _ ε],
         intro conjI, fact ε > 0, clarify)
    fix x assume "x  x0" "¦x - x0¦ < ε"
    with ε_props have [simp]: "poly p x  0" "poly q x  0"
        "sgn (poly (pderiv p) x) = sgn (poly q x)" by auto
    show "sgn (poly (p*q) x) = (if x > x0 then 1 else -1)"
    proof (cases "x  x0")
      case True
        with x  x0 have "x > x0" by simp
        from poly_MVT[OF this, of p] guess ξ ..
        note ξ_props = this
        with ¦x - x0¦ < ε ‹poly p x0 = 0 x > x0 ε_props
            have "¦ξ - x0¦ < ε" "sgn (poly p x) = sgn (x - x0) * sgn (poly q ξ)"
            by (auto simp add: q_pderiv sgn_mult)
        moreover from ξ_props ε_props ¦x - x0¦ < ε
            have "t. ξ  t  t  x  poly q t  0" by auto
        hence "sgn (poly q ξ) = sgn (poly q x)" using ξ_props ε_props
            by (intro no_roots_inbetween_imp_same_sign, simp_all)
        ultimately show ?thesis using True x  x0 ε_props ξ_props
            by (auto simp: sgn_mult sqr_pos)
    next
      case False
        hence "x < x0" by simp
        hence sgn: "sgn (x - x0) = -1" by simp
        from poly_MVT[OF x < x0, of p] guess ξ ..
        note ξ_props = this
        with ¦x - x0¦ < ε ‹poly p x0 = 0 x < x0 ε_props
            have "¦ξ - x0¦ < ε" "poly p x = (x - x0) * poly (pderiv p) ξ"
                 "poly p ξ  0" by (auto simp: field_simps)
        hence "sgn (poly p x) = sgn (x - x0) * sgn (poly q ξ)"
            using ε_props ξ_props by (auto simp: q_pderiv sgn_mult)
        moreover from ξ_props ε_props ¦x - x0¦ < ε
            have "t. x  t  t  ξ  poly q t  0" by auto
        hence "sgn (poly q ξ) = sgn (poly q x)" using ξ_props ε_props
            by (rule_tac sym, intro no_roots_inbetween_imp_same_sign, simp_all)
        ultimately show ?thesis using False x  x0
            by (auto simp: sgn_mult sqr_pos)
    qed
  qed
qed

lemma sturm_firsttwo_signs:
  fixes ps :: "real poly list"
  assumes squarefree: "rsquarefree p"
  assumes p_0: "poly p (x0::real) = 0"
  shows "eventually (λx. sgn (poly (p * sturm p ! 1) x) =
             (if x > x0 then 1 else -1)) (at x0)"
proof-
  from assms have [simp]: "p  0" by (auto simp add: rsquarefree_roots)
  with squarefree p_0 have [simp]: "pderiv p  0"
      by (auto simp  add:rsquarefree_roots)
  from assms show ?thesis
      by (intro sturm_firsttwo_signs_aux,
          simp_all add: rsquarefree_roots)
qed


text ‹
  The construction also obviously fulfils the property about three
  adjacent polynomials in the sequence.
›

lemma sturm_signs:
  assumes squarefree: "rsquarefree p"
  assumes i_in_range: "i < length (sturm (p :: real poly)) - 2"
  assumes q_0: "poly (sturm p ! (i+1)) x = 0" (is "poly ?q x = 0")
  shows "poly (sturm p ! (i+2)) x * poly (sturm p ! i) x < 0"
            (is "poly ?p x * poly ?r x < 0")
proof-
  from sturm_indices[OF i_in_range]
      have "sturm p ! (i+2) = - (sturm p ! i mod sturm p ! (i+1))"
           (is "?r = - (?p mod ?q)") .
  hence "-?r = ?p mod ?q" by simp
  with div_mult_mod_eq[of ?p ?q] have "?p div ?q * ?q - ?r = ?p" by simp
  hence "poly (?p div ?q) x * poly ?q x - poly ?r x = poly ?p x"
      by (metis poly_diff poly_mult)
  with q_0 have r_x: "poly ?r x = -poly ?p x" by simp
  moreover have sqr_pos: "x::real. x  0  x * x > 0" apply (case_tac "x  0")
      by (simp_all add: mult_neg_neg)
  from sturm_adjacent_root_not_squarefree[of i p] assms r_x
      have "poly ?p x * poly ?p x > 0" by (force intro: sqr_pos)
  ultimately show "poly ?r x * poly ?p x < 0" by simp
qed


text ‹
  Finally, if $p$ contains no multiple roots, @{term "sturm p"}, i.e.
  the canonical Sturm sequence for $p$, is a Sturm sequence
  and can be used to determine the number of roots of $p$.
›
lemma sturm_seq_sturm[simp]:
   assumes "rsquarefree p"
   shows "sturm_seq (sturm p) p"
proof
  show "sturm p  []" by simp
  show "hd (sturm p) = p" by simp
  show "length (sturm p)  2" by simp
  from assms show "x. ¬(poly p x = 0  poly (sturm p ! 1) x = 0)"
      by (simp add: rsquarefree_roots)
next
  fix x :: real and y :: real
  have "degree (last (sturm p)) = 0" by simp
  then obtain c where "last (sturm p) = [:c:]"
      by (cases "last (sturm p)", simp split: if_split_asm)
  thus "x y. sgn (poly (last (sturm p)) x) =
            sgn (poly (last (sturm p)) y)" by simp
next
  from sturm_firsttwo_signs[OF assms]
    show "x0. poly p x0 = 0 
         eventually (λx. sgn (poly (p*sturm p ! 1) x) =
                         (if x > x0 then 1 else -1)) (at x0)" by simp
next
  from sturm_signs[OF assms]
    show "i x. i < length (sturm p) - 2; poly (sturm p ! (i + 1)) x = 0
           poly (sturm p ! (i + 2)) x * poly (sturm p ! i) x < 0" by simp
qed


subsubsection ‹Canonical squarefree Sturm sequence›

text ‹
  The previous construction does not work for polynomials with multiple roots,
  but we can simply ``divide away'' multiple roots by dividing $p$ by the
  GCD of $p$ and $p'$. The resulting polynomial has the same roots as $p$,
  but with multiplicity 1, allowing us to again use the canonical construction.
›
definition sturm_squarefree where
  "sturm_squarefree p = sturm (p div (gcd p (pderiv p)))"

lemma sturm_squarefree_not_Nil[simp]: "sturm_squarefree p  []"
  by (simp add: sturm_squarefree_def)


lemma sturm_seq_sturm_squarefree:
  assumes [simp]: "p  0"
  defines [simp]: "p'  p div gcd p (pderiv p)"
  shows "sturm_seq (sturm_squarefree p) p'"
proof
  have "rsquarefree p'"
  proof (subst rsquarefree_roots, clarify)
    fix x assume "poly p' x = 0" "poly (pderiv p') x = 0"
    hence "[:-x,1:] dvd gcd p' (pderiv p')" by (simp add: poly_eq_0_iff_dvd)
    also from poly_div_gcd_squarefree(1)[OF assms(1)]
        have "gcd p' (pderiv p') = 1" by simp
    finally show False by (simp add: poly_eq_0_iff_dvd[symmetric])
  qed

  from sturm_seq_sturm[OF ‹rsquarefree p']
      interpret sturm_seq: sturm_seq "sturm_squarefree p" p'
      by (simp add: sturm_squarefree_def)

  show "x y. sgn (poly (last (sturm_squarefree p)) x) =
      sgn (poly (last (sturm_squarefree p)) y)" by simp
  show "sturm_squarefree p  []" by simp
  show "hd (sturm_squarefree p) = p'" by (simp add: sturm_squarefree_def)
  show "length (sturm_squarefree p)  2" by simp

  have [simp]: "sturm_squarefree p ! 0 = p'"
               "sturm_squarefree p ! Suc 0 = pderiv p'"
      by (simp_all add: sturm_squarefree_def)

  from ‹rsquarefree p'
      show "x. ¬ (poly p' x = 0  poly (sturm_squarefree p ! 1) x = 0)"
      by (simp add: rsquarefree_roots)

  from sturm_seq.signs show "i x. i < length (sturm_squarefree p) - 2;
                                 poly (sturm_squarefree p ! (i + 1)) x = 0
                                  poly (sturm_squarefree p ! (i + 2)) x *
                                         poly (sturm_squarefree p ! i) x < 0" .

  from sturm_seq.deriv show "x0. poly p' x0 = 0 
         eventually (λx. sgn (poly (p' * sturm_squarefree p ! 1) x) =
                         (if x > x0 then 1 else -1)) (at x0)" .
qed


subsubsection ‹Optimisation for multiple roots›

text ‹
  We can also define the following non-canonical Sturm sequence that
  is obtained by taking the canonical Sturm sequence of $p$
  (possibly with multiple roots) and then dividing the entire
  sequence by the GCD of $p$ and its derivative.
›
definition sturm_squarefree' where
"sturm_squarefree' p = (let d = gcd p (pderiv p)
                         in map (λp'. p' div d) (sturm p))"

text ‹
  This construction also has all the desired properties:
›

lemma sturm_squarefree'_adjacent_root_propagate_left:
  assumes "p  0"
  assumes "i < length (sturm_squarefree' (p :: real poly)) - 1"
  assumes "poly (sturm_squarefree' p ! i) x = 0"
      and "poly (sturm_squarefree' p ! (i + 1)) x = 0"
  shows "ji+1. poly (sturm_squarefree' p ! j) x = 0"
proof (intro sturm_adjacent_root_aux[OF assms(2,3,4)], goal_cases)
  case prems: (1 i x)
    define q where "q = sturm p ! i"
    define r where "r = sturm p ! (Suc i)"
    define s where "s = sturm p ! (Suc (Suc i))"
    define d where "d = gcd p (pderiv p)"
    define q' r' s' where "q' = q div d" and "r' = r div d" and "s' = s div d"
    from p  0 have "d  0" unfolding d_def by simp
    from prems(1) have i_in_range: "i < length (sturm p) - 2"
        unfolding sturm_squarefree'_def Let_def by simp
    have [simp]: "d dvd q" "d dvd r" "d dvd s" unfolding q_def r_def s_def d_def
        using i_in_range by (auto intro: sturm_gcd)
    hence qrs_simps: "q = q' * d" "r = r' * d" "s = s' * d"
        unfolding q'_def r'_def s'_def by (simp_all)
    with prems(2) i_in_range have r'_0: "poly r' x = 0"
        unfolding r'_def r_def d_def sturm_squarefree'_def Let_def by simp
    hence r_0: "poly r x = 0" by (simp add: r = r' * d)
    from sturm_indices[OF i_in_range] have "q = q div r * r - s"
        unfolding q_def r_def s_def by (simp add: div_mult_mod_eq)
    hence "q' = (q div r * r - s) div d" by (simp add: q'_def)
    also have "... = (q div r * r) div d - s'"
      by (simp add: s'_def poly_div_diff_left)
    also have "... = q div r * r' - s'"
        using dvd_div_mult[OF d dvd r, of "q div r"]
        by (simp add: algebra_simps r'_def)
    also have "q div r = q' div r'" by (simp add: qrs_simps d  0)
    finally have "poly q' x = poly (q' div r' * r' - s') x" by simp
    also from r'_0 have "... = -poly s' x" by simp
    finally have "poly s' x = -poly q' x" by simp
    thus ?case using i_in_range
        unfolding q'_def s'_def q_def s_def sturm_squarefree'_def Let_def
        by (simp add: d_def sgn_minus)
qed

lemma sturm_squarefree'_adjacent_roots:
  assumes "p  0"
           "i < length (sturm_squarefree' (p :: real poly)) - 1"
          "poly (sturm_squarefree' p ! i) x = 0"
          "poly (sturm_squarefree' p ! (i + 1)) x = 0"
  shows False
proof-
  define d where "d = gcd p (pderiv p)"
  from sturm_squarefree'_adjacent_root_propagate_left[OF assms]
      have "poly (sturm_squarefree' p ! 0) x = 0"
           "poly (sturm_squarefree' p ! 1) x = 0" by auto
  hence "poly (p div d) x = 0" "poly (pderiv p div d) x = 0"
      using assms(2)
      unfolding sturm_squarefree'_def Let_def d_def by auto
  moreover from div_gcd_coprime assms(1)
      have "coprime (p div d) (pderiv p div d)" unfolding d_def by auto
  ultimately show False using coprime_imp_no_common_roots by auto
qed

lemma sturm_squarefree'_signs:
  assumes "p  0"
  assumes i_in_range: "i < length (sturm_squarefree' (p :: real poly)) - 2"
  assumes q_0: "poly (sturm_squarefree' p ! (i+1)) x = 0" (is "poly ?q x = 0")
  shows "poly (sturm_squarefree' p ! (i+2)) x *
         poly (sturm_squarefree' p ! i) x < 0"
            (is "poly ?r x * poly ?p x < 0")
proof-
  define d where "d = gcd p (pderiv p)"
  with p  0 have [simp]: "d  0" by simp
  from poly_div_gcd_squarefree(1)[OF p  0]
       coprime_imp_no_common_roots
      have rsquarefree: "rsquarefree (p div d)"
      by (auto simp: rsquarefree_roots d_def)

  from i_in_range have i_in_range': "i < length (sturm p) - 2"
      unfolding sturm_squarefree'_def by simp
  hence "d dvd (sturm p ! i)" (is "d dvd ?p'")
        "d dvd (sturm p ! (Suc i))" (is "d dvd ?q'")
        "d dvd (sturm p ! (Suc (Suc i)))" (is "d dvd ?r'")
      unfolding d_def by (auto intro: sturm_gcd)
  hence pqr_simps: "?p' = ?p * d" "?q' = ?q * d" "?r' = ?r * d"
    unfolding sturm_squarefree'_def Let_def d_def using i_in_range'
    by (auto simp: dvd_div_mult_self)
  with q_0 have q'_0: "poly ?q' x = 0" by simp
  from sturm_indices[OF i_in_range']
      have "sturm p ! (i+2) = - (sturm p ! i mod sturm p ! (i+1))" .
  hence "-?r' = ?p' mod ?q'" by simp
  with div_mult_mod_eq[of ?p' ?q'] have "?p' div ?q' * ?q' - ?r' = ?p'" by simp
  hence "d*(?p div ?q * ?q - ?r) = d* ?p" by (simp add: pqr_simps algebra_simps)
  hence "?p div ?q * ?q - ?r = ?p" by simp
  hence "poly (?p div ?q) x * poly ?q x - poly ?r x = poly ?p x"
      by (metis poly_diff poly_mult)
  with q_0 have r_x: "poly ?r x = -poly ?p x" by simp

  from sturm_squarefree'_adjacent_roots[OF p  0] i_in_range q_0
      have "poly ?p x  0" by force
  moreover have sqr_pos: "x::real. x  0  x * x > 0" apply (case_tac "x  0")
      by (simp_all add: mult_neg_neg)
  ultimately show ?thesis using r_x by simp
qed


text ‹
  This approach indeed also yields a valid squarefree Sturm sequence
  for the polynomial $p/\text{gcd}(p,p')$.
›
lemma sturm_seq_sturm_squarefree':
  assumes "(p :: real poly)  0"
  defines "d  gcd p (pderiv p)"
  shows "sturm_seq (sturm_squarefree' p) (p div d)"
      (is "sturm_seq ?ps' ?p'")
proof
  show "?ps'  []" "hd ?ps' = ?p'" "2  length ?ps'"
      by (simp_all add: sturm_squarefree'_def d_def hd_map)

  from assms have "d  0" by simp
  {
    have "d dvd last (sturm p)" unfolding d_def
        by (rule sturm_gcd, simp)
    hence *: "last (sturm p) = last ?ps' * d"
        by (simp add: sturm_squarefree'_def last_map d_def dvd_div_mult_self)
    then have "last ?ps' dvd last (sturm p)" by simp
    with * dvd_imp_degree_le[OF this] have "degree (last ?ps')  degree (last (sturm p))"
        using d  0 by (cases "last ?ps' = 0") auto
    hence "degree (last ?ps') = 0" by simp
    then obtain c where "last ?ps' = [:c:]"
        by (cases "last ?ps'", simp split: if_split_asm)
    thus "x y. sgn (poly (last ?ps') x) = sgn (poly (last ?ps') y)" by simp
  }

  have squarefree: "rsquarefree ?p'" using p  0
    by (subst rsquarefree_roots, unfold d_def,
        intro allI coprime_imp_no_common_roots poly_div_gcd_squarefree)
  have [simp]: "sturm_squarefree' p ! Suc 0 = pderiv p div d"
      unfolding sturm_squarefree'_def Let_def sturm_def d_def
          by (subst sturm_aux.simps, simp)
  have coprime: "coprime ?p' (pderiv p div d)"
      unfolding d_def using div_gcd_coprime p  0 by blast
  thus squarefree':
      "x. ¬ (poly (p div d) x = 0  poly (sturm_squarefree' p ! 1) x = 0)"
      using coprime_imp_no_common_roots by simp

  from sturm_squarefree'_signs[OF p  0]
      show "i x. i < length ?ps' - 2; poly (?ps' ! (i + 1)) x = 0
                 poly (?ps' ! (i + 2)) x * poly (?ps' ! i) x < 0" .

  have [simp]: "?p'  0" using squarefree by (simp add: rsquarefree_def)
  have A: "?p' = ?ps' ! 0" "pderiv p div d = ?ps' ! 1"
      by (simp_all add: sturm_squarefree'_def Let_def d_def sturm_def,
          subst sturm_aux.simps, simp)
  have [simp]: "?ps' ! 0  0" using squarefree
      by (auto simp: A rsquarefree_def)

  fix x0 :: real
  assume "poly ?p' x0 = 0"
  hence "poly p x0 = 0" using poly_div_gcd_squarefree(2)[OF p  0]
      unfolding d_def by simp
  hence "pderiv p  0" using p  0 by (auto dest: pderiv_iszero)
  with p  0 ‹poly p x0 = 0
      have A: "eventually (λx. sgn (poly (p * pderiv p) x) =
                              (if x0 < x then 1 else -1)) (at x0)"
      by (intro sturm_firsttwo_signs_aux, simp_all)
  note ev = eventually_conj[OF A poly_neighbourhood_without_roots[OF d  0]]

  show "eventually (λx. sgn (poly (p div d * sturm_squarefree' p ! 1) x) =
                        (if x0 < x then 1 else -1)) (at x0)"
  proof (rule eventually_mono[OF ev], goal_cases)
      have [intro]:
          "a (b::real). b  0  a < 0  a / (b * b) < 0"
          "a (b::real). b  0  a > 0  a / (b * b) > 0"
          by ((case_tac "b > 0",
              auto simp: mult_neg_neg field_simps) [])+
    case prems: (1 x)
      hence  [simp]: "poly d x * poly d x > 0"
           by (cases "poly d x > 0", auto simp: mult_neg_neg)
      from poly_div_gcd_squarefree_aux(2)[OF ‹pderiv p  0]
          have "poly (p div d) x = 0  poly p x = 0" by (simp add: d_def)
      moreover have "d dvd p" "d dvd pderiv p" unfolding d_def by simp_all
      ultimately show ?case using prems
          by (auto simp: sgn_real_def poly_div not_less[symmetric]
                         zero_less_divide_iff split: if_split_asm)
  qed
qed


text ‹
  This construction is obviously more expensive to compute than the one that \emph{first}
  divides $p$ by $\text{gcd}(p,p')$ and \emph{then} applies the canonical construction.
  In this construction, we \emph{first} compute the canonical Sturm sequence of $p$ as if
  it had no multiple roots and \emph{then} divide by the GCD.
  However, it can be seen quite easily that unless $x$ is a multiple root of $p$,
  i.\,e. as long as $\text{gcd}(P,P')\neq 0$, the number of sign changes in a sequence of
  polynomials does not actually change when we divide the polynomials by $\text{gcd}(p,p')$.\\
  There\-fore we can use the ca\-no\-ni\-cal Sturm se\-quence even in the non-square\-free
  case as long as the borders of the interval we are interested in are not multiple roots
  of the polynomial.
›

lemma sign_changes_mult_aux:
  assumes "d  (0::real)"
  shows "length (remdups_adj (filter (λx. x  0) (map ((*) d  f) xs))) =
         length (remdups_adj (filter (λx. x  0) (map f xs)))"
proof-
  from assms have inj: "inj ((*) d)" by (auto intro: injI)
  from assms have [simp]: "filter (λx. ((*) d  f) x  0) = filter (λx. f x  0)"
                          "filter ((λx. x  0)  f) = filter (λx. f x  0)"
      by (simp_all add: o_def)
  have "filter (λx. x  0) (map ((*) d  f) xs) =
        map ((*) d  f) (filter (λx. ((*) d  f) x  0) xs)"
      by (simp add: filter_map o_def)
  thus ?thesis using remdups_adj_map_injective[OF inj] assms
      by (simp add: filter_map map_map[symmetric] del: map_map)
qed

lemma sturm_sturm_squarefree'_same_sign_changes:
  fixes p :: "real poly"
  defines "ps  sturm p" and "ps'  sturm_squarefree' p"
  shows "poly p x  0  poly (pderiv p) x  0 
             sign_changes ps' x = sign_changes ps x"
        "p  0  sign_changes_inf ps' = sign_changes_inf ps"
        "p  0  sign_changes_neg_inf ps' = sign_changes_neg_inf ps"
proof-
  define d where "d = gcd p (pderiv p)"
  define p' where "p' = p div d"
  define s' where "s' = poly_inf d"
  define s'' where "s'' = poly_neg_inf d"

  {
    fix x :: real and q :: "real poly"
    assume "q  set ps"
    hence "d dvd q" unfolding d_def ps_def using sturm_gcd by simp
    hence q_prod: "q = (q div d) * d" unfolding p'_def d_def
        by (simp add: algebra_simps dvd_mult_div_cancel)

    have "poly q x = poly d x * poly (q div d) x"  by (subst q_prod, simp)
    hence s1: "sgn (poly q x) = sgn (poly d x) * sgn (poly (q div d) x)"
        by (subst q_prod, simp add: sgn_mult)
    from poly_inf_mult have s2: "poly_inf q = s' * poly_inf (q div d)"
        unfolding s'_def by (subst q_prod, simp)
    from poly_inf_mult have s3: "poly_neg_inf q = s'' * poly_neg_inf (q div d)"
        unfolding s''_def by (subst q_prod, simp)
    note s1 s2 s3
  }
  note signs = this

  {
    fix f :: "real poly  real" and s :: real
    assume f: "q. q  set ps  f q = s * f (q div d)" and s: "s  0"
    hence "inverse s  0" by simp
    {fix q assume "q  set ps"
     hence "f (q div d) = inverse s * f q"
         by (subst f[of q], simp_all add: s)
    } note f' = this
    have "length (remdups_adj [xmap f (map (λq. q div d) ps). x  0]) - 1 =
           length (remdups_adj [xmap (λq. f (q div d)) ps . x  0]) - 1"
        by (simp only: sign_changes_def o_def map_map)
    also have "map (λq. q div d) ps = ps'"
        by (simp add: ps_def ps'_def sturm_squarefree'_def Let_def d_def)
    also from f' have "map (λq. f (q div d)) ps =
                      map (λx. ((*)(inverse s)  f) x) ps" by (simp add: o_def)
    also note sign_changes_mult_aux[OF ‹inverse s  0, of f ps]
    finally have
        "length (remdups_adj [xmap f ps' . x  0]) - 1 =
         length (remdups_adj [xmap f ps . x  0]) - 1" by simp
  }
  note length_remdups_adj = this

  {
    fix x assume A: "poly p x  0  poly (pderiv p) x  0"
    have "d dvd p" "d dvd pderiv p" unfolding d_def by simp_all
    with A have "sgn (poly d x)  0"
        by (auto simp add: sgn_zero_iff elim: dvdE)
    thus "sign_changes ps' x = sign_changes ps x" using signs(1)
        unfolding sign_changes_def
        by (intro length_remdups_adj[of "λq. sgn (poly q x)"], simp_all)
  }

  assume "p  0"
  hence "d  0" unfolding d_def by simp
  hence "s'  0" "s''  0" unfolding s'_def s''_def by simp_all
  from length_remdups_adj[of poly_inf s', OF signs(2) s'  0]
      show "sign_changes_inf ps' = sign_changes_inf ps"
      unfolding sign_changes_inf_def .
  from length_remdups_adj[of poly_neg_inf s'', OF signs(3) s''  0]
      show "sign_changes_neg_inf ps' = sign_changes_neg_inf ps"
      unfolding sign_changes_neg_inf_def .
qed



subsection ‹Root-counting functions›

text ‹
  With all these results, we can now define functions that count roots
  in bounded and unbounded intervals:
›

definition count_roots_between where
"count_roots_between p a b = (if a  b  p  0 then
  (let ps = sturm_squarefree p
    in sign_changes ps a - sign_changes ps b) else 0)"

definition count_roots where
"count_roots p = (if (p::real poly) = 0 then 0 else
  (let ps = sturm_squarefree p
    in sign_changes_neg_inf ps - sign_changes_inf ps))"

definition count_roots_above where
"count_roots_above p a = (if (p::real poly) = 0 then 0 else
  (let ps = sturm_squarefree p
    in sign_changes ps a - sign_changes_inf ps))"

definition count_roots_below where
"count_roots_below p a = (if (p::real poly) = 0 then 0 else
  (let ps = sturm_squarefree p
    in sign_changes_neg_inf ps - sign_changes ps a))"


lemma count_roots_between_correct:
  "count_roots_between p a b = card {x. a < x  x  b  poly p x = 0}"
proof (cases "p  0  a  b")
  case False
    note False' = this
    hence "card {x. a < x  x  b  poly p x = 0} = 0"
    proof (cases "a < b")
      case True
        with False have [simp]: "p = 0" by simp
        have subset: "{a<..<b}  {x. a < x  x  b  poly p x = 0}" by auto
        from infinite_Ioo[OF True] have "¬finite {a<..<b}" .
        hence "¬finite {x. a < x  x  b  poly p x = 0}"
            using finite_subset[OF subset] by blast
        thus ?thesis by simp
    next
      case False
        with False' show ?thesis by (auto simp: not_less card_eq_0_iff)
    qed
    thus ?thesis unfolding count_roots_between_def Let_def using False by auto
next
  case True
  hence "p  0" "a  b" by simp_all
  define p' where "p' = p div (gcd p (pderiv p))"
  from poly_div_gcd_squarefree(1)[OF p  0] have "p'  0"
      unfolding p'_def by clarsimp

  from sturm_seq_sturm_squarefree[OF p  0]
      interpret sturm_seq "sturm_squarefree p" p'
      unfolding p'_def .
  from poly_roots_finite[OF p'  0]
      have "finite {x. a < x  x  b  poly p' x = 0}" by fast
  have "count_roots_between p a b = card {x. a < x  x  b  poly p' x = 0}"
      unfolding count_roots_between_def Let_def
      using True count_roots_between[OF p'  0 a  b] by simp
  also from poly_div_gcd_squarefree(2)[OF p  0]
      have "{x. a < x  x  b  poly p' x = 0} =
            {x. a < x  x  b  poly p x = 0}" unfolding p'_def by blast
  finally show ?thesis .
qed

lemma count_roots_correct:
  fixes p :: "real poly"
  shows "count_roots p = card {x. poly p x = 0}" (is "_ = card ?S")
proof (cases "p = 0")
  case True
    with finite_subset[of "{0<..<1}" ?S]
    have "¬finite {x. poly p x = 0}" by (auto simp: infinite_Ioo)
    thus ?thesis by (simp add: count_roots_def True)
next
  case False
  define p' where "p' = p div (gcd p (pderiv p))"
  from poly_div_gcd_squarefree(1)[OF p  0] have "p'  0"
      unfolding p'_def by clarsimp

  from sturm_seq_sturm_squarefree[OF p  0]
      interpret sturm_seq "sturm_squarefree p" p'
      unfolding p'_def .
  from count_roots[OF p'  0]
      have "count_roots p = card {x. poly p' x = 0}"
      unfolding count_roots_def Let_def by (simp add: p  0)
  also from poly_div_gcd_squarefree(2)[OF p  0]
      have "{x. poly p' x = 0} = {x. poly p x = 0}" unfolding p'_def by blast
  finally show ?thesis .
qed

lemma count_roots_above_correct:
  fixes p :: "real poly"
  shows "count_roots_above p a = card {x. x > a  poly p x = 0}"
         (is "_ = card ?S")
proof (cases "p = 0")
  case True
  with finite_subset[of "{a<..<a+1}" ?S]
    have "¬finite {x. x > a  poly p x = 0}" by (auto simp: infinite_Ioo subset_eq)
  thus ?thesis by (simp add: count_roots_above_def True)
next
  case False
  define p' where "p' = p div (gcd p (pderiv p))"
  from poly_div_gcd_squarefree(1)[OF p  0] have "p'  0"
      unfolding p'_def by clarsimp

  from sturm_seq_sturm_squarefree[OF p  0]
      interpret sturm_seq "sturm_squarefree p" p'
      unfolding p'_def .
  from count_roots_above[OF p'  0]
      have "count_roots_above p a = card {x. x > a  poly p' x = 0}"
      unfolding count_roots_above_def Let_def by (simp add: p  0)
  also from poly_div_gcd_squarefree(2)[OF p  0]
      have "{x. x > a  poly p' x = 0} = {x. x > a  poly p x = 0}"
      unfolding p'_def by blast
  finally show ?thesis .
qed

lemma count_roots_below_correct:
  fixes p :: "real poly"
  shows "count_roots_below p a = card {x. x  a  poly p x = 0}"
         (is "_ = card ?S")
proof (cases "p = 0")
  case True
    with finite_subset[of "{a - 1<..<a}" ?S]
        have "¬finite {x. x  a  poly p x = 0}" by (auto simp: infinite_Ioo subset_eq)
    thus ?thesis by (simp add: count_roots_below_def True)
next
  case False
  define p' where "p' = p div (gcd p (pderiv p))"
  from poly_div_gcd_squarefree(1)[OF p  0] have "p'  0"
      unfolding p'_def by clarsimp

  from sturm_seq_sturm_squarefree[OF p  0]
      interpret sturm_seq "sturm_squarefree p" p'
      unfolding p'_def .
  from count_roots_below[OF p'  0]
      have "count_roots_below p a = card {x. x  a  poly p' x = 0}"
      unfolding count_roots_below_def Let_def by (simp add: p  0)
  also from poly_div_gcd_squarefree(2)[OF p  0]
      have "{x. x  a  poly p' x = 0} = {x. x  a  poly p x = 0}"
      unfolding p'_def by blast
  finally show ?thesis .
qed

text ‹
  The optimisation explained above can be used to prove more efficient code equations that
  use the more efficient construction in the case that the interval borders are not
  multiple roots:
›

lemma count_roots_between[code]:
  "count_roots_between p a b =
     (let q = pderiv p
       in if a > b  p = 0 then 0
       else if (poly p a  0  poly q a  0)  (poly p b  0  poly q b  0)
            then (let ps = sturm p
                   in sign_changes ps a - sign_changes ps b)
            else (let ps = sturm_squarefree p
                   in sign_changes ps a - sign_changes ps b))"
proof (cases "a > b  p = 0")
  case True
    thus ?thesis by (auto simp add: count_roots_between_def Let_def)
next
  case False
    note False1 = this
    hence "a  b" "p  0" by simp_all
    thus ?thesis
    proof (cases "(poly p a  0  poly (pderiv p) a  0) 
                  (poly p b  0  poly (pderiv p) b  0)")
    case False
      thus ?thesis using False1
          by (auto simp add: Let_def count_roots_between_def)
    next
    case True
      hence A: "poly p a  0  poly (pderiv p) a  0" and
            B: "poly p b  0  poly (pderiv p) b  0" by auto
      define d where "d = gcd p (pderiv p)"
      from p  0 have [simp]: "p div d  0"
          using poly_div_gcd_squarefree(1)[OF p  0] by (auto simp add: d_def)
      from sturm_seq_sturm_squarefree'[OF p  0]
          interpret sturm_seq "sturm_squarefree' p" "p div d"
          unfolding sturm_squarefree'_def Let_def d_def .
      note count_roots_between_correct
      also have "{x. a < x  x  b  poly p x = 0} =
                 {x. a < x  x  b  poly (p div d) x = 0}"
          unfolding d_def using poly_div_gcd_squarefree(2)[OF p  0] by simp
      also note count_roots_between[OF p div d  0 a  b, symmetric]
      also note sturm_sturm_squarefree'_same_sign_changes(1)[OF A]
      also note sturm_sturm_squarefree'_same_sign_changes(1)[OF B]
      finally show ?thesis using True False by (simp add: Let_def)
    qed
qed


lemma count_roots_code[code]:
  "count_roots (p::real poly) =
    (if p = 0 then 0
     else let ps = sturm p
           in sign_changes_neg_inf ps - sign_changes_inf ps)"
proof (cases "p = 0", simp add: count_roots_def)
  case False
    define d where "d = gcd p (pderiv p)"
    from p  0 have [simp]: "p div d  0"
        using poly_div_gcd_squarefree(1)[OF p  0] by (auto simp add: d_def)
    from sturm_seq_sturm_squarefree'[OF p  0]
        interpret sturm_seq "sturm_squarefree' p" "p div d"
        unfolding sturm_squarefree'_def Let_def d_def .

    note count_roots_correct
    also have "{x. poly p x = 0} = {x. poly (p div d) x = 0}"
        unfolding d_def using poly_div_gcd_squarefree(2)[OF p  0] by simp
    also note count_roots[OF p div d  0, symmetric]
    also note sturm_sturm_squarefree'_same_sign_changes(2)[OF p  0]
    also note sturm_sturm_squarefree'_same_sign_changes(3)[OF p  0]
    finally show ?thesis using False unfolding Let_def by simp
qed


lemma count_roots_above_code[code]:
  "count_roots_above p a =
     (let q = pderiv p
       in if p = 0 then 0
       else if poly p a  0  poly q a  0
            then (let ps = sturm p
                   in sign_changes ps a - sign_changes_inf ps)
            else (let ps = sturm_squarefree p
                   in sign_changes ps a - sign_changes_inf ps))"
proof (cases "p = 0")
  case True
    thus ?thesis by (auto simp add: count_roots_above_def Let_def)
next
  case False
    note False1 = this
    hence "p  0" by simp_all
    thus ?thesis
    proof (cases "(poly p a  0  poly (pderiv p) a  0)")
    case False
      thus ?thesis using False1
          by (auto simp add: Let_def count_roots_above_def)
    next
    case True
      hence A: "poly p a  0  poly (pderiv p) a  0" by simp
      define d where "d = gcd p (pderiv p)"
      from p  0 have [simp]: "p div d  0"
          using poly_div_gcd_squarefree(1)[OF p  0] by (auto simp add: d_def)
      from sturm_seq_sturm_squarefree'[OF p  0]
          interpret sturm_seq "sturm_squarefree' p" "p div d"
          unfolding sturm_squarefree'_def Let_def d_def .
      note count_roots_above_correct
      also have "{x. a < x  poly p x = 0} =
                 {x. a < x  poly (p div d) x = 0}"
          unfolding d_def using poly_div_gcd_squarefree(2)[OF p  0] by simp
      also note count_roots_above[OF p div d  0, symmetric]
      also note sturm_sturm_squarefree'_same_sign_changes(1)[OF A]
      also note sturm_sturm_squarefree'_same_sign_changes(2)[OF p  0]
      finally show ?thesis using True False by (simp add: Let_def)
    qed
qed

lemma count_roots_below_code[code]:
  "count_roots_below p a =
     (let q = pderiv p
       in if p = 0 then 0
       else if poly p a  0  poly q a  0
            then (let ps = sturm p
                   in sign_changes_neg_inf ps - sign_changes ps a)
            else (let ps = sturm_squarefree p
                   in sign_changes_neg_inf ps - sign_changes ps a))"
proof (cases "p = 0")
  case True
    thus ?thesis by (auto simp add: count_roots_below_def Let_def)
next
  case False
    note False1 = this
    hence "p  0" by simp_all
    thus ?thesis
    proof (cases "(poly p a  0  poly (pderiv p) a  0)")
    case False
      thus ?thesis using False1
          by (auto simp add: Let_def count_roots_below_def)
    next
    case True
      hence A: "poly p a  0  poly (pderiv p) a  0" by simp
      define d where "d = gcd p (pderiv p)"
      from p  0 have [simp]: "p div d  0"
          using poly_div_gcd_squarefree(1)[OF p  0] by (auto simp add: d_def)
      from sturm_seq_sturm_squarefree'[OF p  0]
          interpret sturm_seq "sturm_squarefree' p" "p div d"
          unfolding sturm_squarefree'_def Let_def d_def .
      note count_roots_below_correct
      also have "{x. x  a  poly p x = 0} =
                 {x. x  a  poly (p div d) x = 0}"
          unfolding d_def using poly_div_gcd_squarefree(2)[OF p  0] by simp
      also note count_roots_below[OF p div d  0, symmetric]
      also note sturm_sturm_squarefree'_same_sign_changes(1)[OF A]
      also note sturm_sturm_squarefree'_same_sign_changes(3)[OF p  0]
      finally show ?thesis using True False by (simp add: Let_def)
    qed
qed

end

Theory Sturm_Method

section ‹The ``sturm'' proof method›

(* Author: Manuel Eberl <eberlm@in.tum.de> *)
theory Sturm_Method
imports Sturm_Theorem
begin

subsection ‹Preliminary lemmas›

text ‹
  In this subsection, we prove lemmas that reduce root counting and
  related statements to simple, computable expressions using the 
  @{term "count_roots"} function family.
›

lemma poly_card_roots_less_leq:
  "card {x. a < x  x  b  poly p x = 0} = count_roots_between p a b"
  by (simp add: count_roots_between_correct)

lemma poly_card_roots_leq_leq:
  "card {x. a  x  x  b  poly p x = 0} = 
       ( count_roots_between p a b + 
      (if (a  b  poly p a = 0  p  0)  (a = b  p = 0) then 1 else 0))"
proof (cases "(a  b  poly p a = 0  p  0)  (a = b  p = 0)")
  case False
    note False' = this
    thus ?thesis
    proof (cases "p = 0")
      case False
        with False' have "poly p a  0  a > b" by auto
        hence "{x. a  x  x  b  poly p x = 0} = 
               {x. a < x  x  b  poly p x = 0}"
        by (auto simp: less_eq_real_def)
        thus ?thesis using poly_card_roots_less_leq False' 
            by (auto split: if_split_asm)
    next
      case True
        have "{x. a  x  x  b} = {a..b}"
             "{x. a < x  x  b} = {a<..b}" by auto
        with True False have "card {x. a < x  x  b} = 0" "card {x. a  x  x  b} = 0"
          by (auto simp add: card_eq_0_iff infinite_Ioc infinite_Icc)
        with True False show ?thesis
            using count_roots_between_correct by (simp add: )
    qed
next
  case True
    note True' = this
    have fin: "finite {x. a  x  x  b  poly p x = 0}" 
    proof (cases "p = 0")
      case True
        with True' have "a = b" by simp
        hence "{x. a  x  x  b  poly p x = 0} = {b}" using True by auto
        thus ?thesis by simp
    next
      case False
        from poly_roots_finite[OF this] show ?thesis by fast
    qed
    with True have "{x. a  x  x  b  poly p x = 0} =
        insert a {x. a < x  x  b  poly p x = 0}" by auto
    hence "card {x. a  x  x  b  poly p x = 0} =
           Suc (card {x. a < x  x  b  poly p x = 0})" using fin by force
    thus ?thesis using True count_roots_between_correct by simp
qed

lemma poly_card_roots_less_less:
  "card {x. a < x  x < b  poly p x = 0} = 
      ( count_roots_between p a b -
              (if poly p b = 0  a < b  p  0 then 1 else 0))"
proof (cases "poly p b = 0  a < b  p  0")
  case False
    note False' = this
    show ?thesis
    proof (cases "p = 0")
      case True
        have [simp]: "{x. a < x  x < b} = {a<..<b}"
                     "{x. a < x  x  b} = {a<..b}" by auto
        with True False have "card {x. a < x  x  b} = 0" "card {x. a < x  x < b} = 0"
          by (auto simp add: card_eq_0_iff infinite_Ioo infinite_Ioc)
        with True False' show ?thesis 
            by (auto simp: count_roots_between_correct)
    next
      case False
        with False' have "{x. a < x  x < b  poly p x = 0} = 
                          {x. a < x  x  b  poly p x = 0}"
          by (auto simp: less_eq_real_def)
      thus ?thesis using poly_card_roots_less_leq False by auto
  qed
next
  case True
    with poly_roots_finite
        have fin: "finite {x. a < x  x < b  poly p x = 0}" by fast
    from True have "{x. a < x  x  b  poly p x = 0} =
        insert b {x. a < x  x < b  poly p x = 0}" by auto
    hence "Suc (card {x. a < x  x < b  poly p x = 0}) =
           card {x. a < x  x  b  poly p x = 0}" using fin by force
    also note count_roots_between_correct[symmetric]
    finally show ?thesis using True by simp
qed

lemma poly_card_roots_leq_less:
  "card {x::real. a  x  x < b  poly p x = 0} =
      ( count_roots_between p a b +
      (if p  0  a < b  poly p a = 0 then 1 else 0) -
      (if p  0  a < b  poly p b = 0 then 1 else 0))"
proof (cases "p = 0  a  b")
  case True
    note True' = this
    show ?thesis
    proof (cases "a  b")
      case False
        hence "{x. a < x  x  b} = {a<..b}"
              "{x. a  x  x < b} = {a..<b}" by auto
        with True False have "card {x. a < x  x  b} = 0" "card {x. a  x  x < b} = 0"
          by (auto simp add: card_eq_0_iff infinite_Ico infinite_Ioc)
        with False True' show ?thesis 
            by (simp add: count_roots_between_correct)
    next
      case True
        with True' have "{x. a  x  x < b  poly p x = 0} = 
                          {x. a < x  x  b  poly p x = 0}"
          by (auto simp: less_eq_real_def)
      thus ?thesis using poly_card_roots_less_leq True by simp
  qed
next
  case False
    let ?A = "{x. a  x  x < b  poly p x = 0}"
    let ?B = "{x. a < x  x  b  poly p x = 0}"
    let ?C = "{x. x = b  poly p x = 0}"
    let ?D = "{x. x = a  poly p a = 0}"
    have CD_if: "?C = (if poly p b = 0 then {b} else {})"
                "?D = (if poly p a = 0 then {a} else {})" by auto
    from False poly_roots_finite 
        have [simp]: "finite ?A" "finite ?B" "finite ?C" "finite ?D"
            by (fast, fast, simp_all)
    
    from False have "?A = (?B  ?D) - ?C" by (auto simp: less_eq_real_def)
    with False have "card ?A = card ?B + (if poly p a = 0 then 1 else 0) -
                       (if poly p b = 0 then 1 else 0)" by (auto simp: CD_if)
    also note count_roots_between_correct[symmetric]
    finally show ?thesis using False by simp
qed

lemma poly_card_roots:
  "card {x::real. poly p x = 0} = count_roots p"
  using count_roots_correct by simp

lemma poly_no_roots:
  "(x. poly p x  0)  ( p  0  count_roots p = 0)"
    by (auto simp: count_roots_correct dest: poly_roots_finite)

lemma poly_pos:
  "(x. poly p x > 0)  ( 
        p  0  poly_inf p = 1  count_roots p = 0)"
  by (simp only: Let_def poly_pos poly_no_roots, blast)


lemma poly_card_roots_greater:
  "card {x::real. x > a  poly p x = 0} = count_roots_above p a"
  using count_roots_above_correct by simp

lemma poly_card_roots_leq:
  "card {x::real. x  a  poly p x = 0} = count_roots_below p a"
  using count_roots_below_correct by simp

lemma poly_card_roots_geq:
  "card {x::real. x  a  poly p x = 0} = (
      count_roots_above p a + (if poly p a = 0  p  0 then 1 else 0))"
proof (cases "poly p a = 0  p  0")
  case False
    hence "card {x. x  a  poly p x = 0} = card {x. x > a  poly p x = 0}"
    proof (cases rule: disjE)
      assume "p = 0"
      have "¬finite {a<..<a+1}"
        by (metis infinite_Ioo less_add_one) 
      moreover have "{a<..<a+1}  {x. x  a  poly p x = 0}"
                    "{a<..<a+1}  {x. x > a  poly p x = 0}" 
          using p = 0 by auto
      ultimately have "¬finite {x. x  a  poly p x = 0}" 
                      "¬finite {x. x > a  poly p x = 0}" 
        by (auto dest!: finite_subset[of "{a<..<a+1}"] simp: infinite_Ioo)
      thus ?thesis by simp
    next
      assume "poly p a  0"
      hence "{x. x  a  poly p x = 0} = {x. x > a  poly p x = 0}" 
          by (auto simp: less_eq_real_def)
      thus ?thesis by simp
    qed auto
    thus ?thesis using False 
        by (auto intro: poly_card_roots_greater)
next
  case True
    hence "finite {x. x > a  poly p x = 0}" using poly_roots_finite by force
    moreover have "{x. x  a  poly p x = 0} = 
                       insert a {x. x > a  poly p x = 0}" using True by auto
    ultimately have "card {x. x  a  poly p x = 0} = 
                         Suc (card {x. x > a  poly p x = 0})"
        using card_insert_disjoint by auto
    thus ?thesis using True by (auto intro: poly_card_roots_greater)
qed

lemma poly_card_roots_less:
  "card {x::real. x < a  poly p x = 0} =
       (count_roots_below p a - (if poly p a = 0  p  0 then 1 else 0))"
proof (cases "poly p a = 0  p  0")
  case False
    hence "card {x. x < a  poly p x = 0} = card {x. x  a  poly p x = 0}"
    proof (cases rule: disjE)
      assume "p = 0"
      have "¬finite {a - 1<..<a}" 
        by (metis infinite_Ioo diff_add_cancel less_add_one) 
      moreover have "{a - 1<..<a}  {x. x  a  poly p x = 0}"
                    "{a - 1<..<a}  {x. x < a  poly p x = 0}" 
          using p = 0 by auto
      ultimately have "¬finite {x. x  a  poly p x = 0}" 
                     "¬finite {x. x < a  poly p x = 0}" 
          by (auto dest: finite_subset[of "{a - 1<..<a}"] simp: infinite_Ioo)
      thus ?thesis by simp
    next
      assume "poly p a  0"
      hence "{x. x < a  poly p x = 0} = {x. x  a  poly p x = 0}" 
          by (auto simp: less_eq_real_def)
      thus ?thesis by simp
    qed auto
    thus ?thesis using False 
        by (auto intro: poly_card_roots_leq)
next
  case True
    hence "finite {x. x < a  poly p x = 0}" using poly_roots_finite by force
    moreover have "{x. x  a  poly p x = 0} = 
                       insert a {x. x < a  poly p x = 0}" using True by auto
    ultimately have "Suc (card {x. x < a  poly p x = 0}) =
                     (card {x. x  a  poly p x = 0})"
        using card_insert_disjoint by auto
    also note count_roots_below_correct[symmetric]
    finally show ?thesis using True by simp
qed


lemma poly_no_roots_less_leq:
  "(x. a < x  x  b  poly p x  0) 
   ((a  b  (p  0  count_roots_between p a b = 0)))"
  by (auto simp: count_roots_between_correct card_eq_0_iff not_le 
           dest: poly_roots_finite)

lemma poly_pos_between_less_leq:
  "(x. a < x  x  b  poly p x > 0) 
   ((a  b  (p  0  poly p b > 0  count_roots_between p a b = 0)))"
  by (simp only: poly_pos_between_less_leq Let_def 
                 poly_no_roots_less_leq, blast)


lemma poly_no_roots_leq_leq:
  "(x. a  x  x  b  poly p x  0) 
   ((a > b  (p  0  poly p a  0  count_roots_between p a b = 0)))"
apply (intro iffI)
apply (force simp add: count_roots_between_correct card_eq_0_iff)
apply (elim conjE disjE, simp, intro allI)
apply (rename_tac x, case_tac "x = a")
apply (auto simp add: count_roots_between_correct card_eq_0_iff
            dest: poly_roots_finite)
done

lemma poly_pos_between_leq_leq:
  "(x. a  x  x  b  poly p x > 0)  
   ((a > b  (p  0  poly p a > 0  
                count_roots_between p a b = 0)))"
by (simp only: poly_pos_between_leq_leq Let_def poly_no_roots_leq_leq, force)



lemma poly_no_roots_less_less:
  "(x. a < x  x < b  poly p x  0)  
   ((a  b  p  0  count_roots_between p a b = 
       (if poly p b = 0 then 1 else 0)))"
proof (standard, goal_cases)
  case A: 1
    show ?case
    proof (cases "a  b")
      case True
      with A show ?thesis by simp
    next
      case False
      with A have [simp]: "p  0" using dense[of a b] by auto
      have B: "{x. a < x  x  b  poly p x = 0} =
                {x. a < x  x < b  poly p x = 0} 
                (if poly p b = 0 then {b} else {})" using A False by auto
      have "count_roots_between p a b =
                 card {x. a < x  x < b  poly p x = 0} +
                (if poly p b = 0 then 1 else 0)"
         by (subst count_roots_between_correct, subst B, subst card_Un_disjoint, 
             rule finite_subset[OF _ poly_roots_finite], blast, simp_all)
      also from A have "{x. a < x  x < b  poly p x = 0} = {}" by simp
      finally show ?thesis by auto
    qed
next
  case prems: 2
    hence "card {x. a < x  x < b  poly p x = 0} = 0"
        by (subst poly_card_roots_less_less, auto simp: count_roots_between_def)
    thus ?case using prems
        by (cases "p = 0", simp, subst (asm) card_eq_0_iff, 
            auto dest: poly_roots_finite)
qed

lemma poly_pos_between_less_less:
  "(x. a < x  x < b  poly p x > 0) 
   ((a  b  (p  0  poly p ((a+b)/2) > 0  
       count_roots_between p a b = (if poly p b = 0 then 1 else 0))))"
  by (simp only: poly_pos_between_less_less Let_def 
                 poly_no_roots_less_less, blast)

lemma poly_no_roots_leq_less:
  "(x. a  x  x < b  poly p x  0) 
   ((a  b  p  0  poly p a  0  count_roots_between p a b = 
       (if a < b  poly p b = 0 then 1 else 0)))"
proof (standard, goal_cases)
  case prems: 1
    hence "x. a < x  x < b  poly p x  0" by simp
    thus ?case using prems by (subst (asm) poly_no_roots_less_less, auto)
next
  case prems: 2
    hence "(b  a  p  0  count_roots_between p a b = 
               (if poly p b = 0 then 1 else 0))" by auto
    thus ?case using prems unfolding Let_def
        by (subst (asm) poly_no_roots_less_less[symmetric, unfolded Let_def], 
        auto split: if_split_asm simp: less_eq_real_def) 
qed

lemma poly_pos_between_leq_less:
  "(x. a  x  x < b  poly p x > 0)  
   ((a  b  (p  0  poly p a > 0  count_roots_between p a b = 
        (if a < b  poly p b = 0 then 1 else 0))))"
 by (simp only: poly_pos_between_leq_less Let_def 
                poly_no_roots_leq_less, force)


lemma poly_no_roots_greater:
  "(x. x > a  poly p x  0) 
       ((p  0  count_roots_above p a = 0))"
proof-
  have "x. ¬ a < x  False" by (metis gt_ex)
  thus ?thesis by (auto simp: count_roots_above_correct card_eq_0_iff
                        intro: poly_roots_finite )
qed

lemma poly_pos_greater:
  "(x. x > a  poly p x > 0)  (
       p  0  poly_inf p = 1  count_roots_above p a = 0)"
  unfolding Let_def
  by (subst poly_pos_greater, subst poly_no_roots_greater, force)

lemma poly_no_roots_leq:
  "(x. x  a  poly p x  0)  
       ( (p  0  count_roots_below p a = 0))"
    by (auto simp: Let_def count_roots_below_correct card_eq_0_iff
             intro: poly_roots_finite)

lemma poly_pos_leq:
  "(x. x  a  poly p x > 0)  
   ( p  0  poly_neg_inf p = 1  count_roots_below p a = 0)"
  by (simp only: poly_pos_leq Let_def poly_no_roots_leq, blast)



lemma poly_no_roots_geq:
  "(x. x  a  poly p x  0) 
       ( (p  0  poly p a  0  count_roots_above p a = 0))"
proof (standard, goal_cases)
  case prems: 1
  hence "x>a. poly p x  0" by simp
  thus ?case using prems by (subst (asm) poly_no_roots_greater, auto)
next
  case prems: 2
  hence "(p  0  count_roots_above p a = 0)" by simp
  thus ?case using prems 
      by (subst (asm) poly_no_roots_greater[symmetric, unfolded Let_def], 
          auto simp: less_eq_real_def)
qed

lemma poly_pos_geq:
  "(x. x  a  poly p x > 0) 
   (p  0  poly_inf p = 1  poly p a  0  count_roots_above p a = 0)"
  by (simp only: poly_pos_geq Let_def poly_no_roots_geq, blast)

lemma poly_no_roots_less:
  "(x. x < a  poly p x  0) 
       ((p  0  count_roots_below p a = (if poly p a = 0 then 1 else 0)))"
proof (standard, goal_cases)
  case prems: 1
  hence "{x. x  a  poly p x = 0} = (if poly p a = 0 then {a} else {})"
      by (auto simp: less_eq_real_def)
  moreover have "x. ¬ x < a  False" by (metis lt_ex)
  ultimately show ?case using prems by (auto simp: count_roots_below_correct)
next
  case prems: 2
  have A: "{x. x  a  poly p x = 0} = {x. x < a  poly p x = 0}  
            (if poly p a = 0 then {a} else {})" by (auto simp: less_eq_real_def)
  have "count_roots_below p a = card {x. x < a  poly p x = 0} +
            (if poly p a = 0 then 1 else 0)" using prems
      by (subst count_roots_below_correct, subst A, subst card_Un_disjoint,
          auto intro: poly_roots_finite)
  with prems have "card {x. x < a  poly p x = 0} = 0" by simp
  thus ?case using prems
      by (subst (asm) card_eq_0_iff, auto intro: poly_roots_finite)
qed

lemma poly_pos_less:
  "(x. x < a  poly p x > 0) 
   (p  0  poly_neg_inf p = 1  count_roots_below p a = 
       (if poly p a = 0 then 1 else 0))"
  by (simp only: poly_pos_less Let_def poly_no_roots_less, blast)


lemmas sturm_card_substs = poly_card_roots poly_card_roots_less_leq 
    poly_card_roots_leq_less poly_card_roots_less_less poly_card_roots_leq_leq
    poly_card_roots_less poly_card_roots_leq poly_card_roots_greater
    poly_card_roots_geq 

lemmas sturm_prop_substs = poly_no_roots poly_no_roots_less_leq 
    poly_no_roots_leq_leq poly_no_roots_less_less poly_no_roots_leq_less
    poly_no_roots_leq poly_no_roots_less poly_no_roots_geq 
    poly_no_roots_greater
    poly_pos poly_pos_greater poly_pos_geq poly_pos_less poly_pos_leq
    poly_pos_between_leq_less poly_pos_between_less_leq
    poly_pos_between_leq_leq poly_pos_between_less_less


subsection ‹Reification›

text ‹
  This subsection defines a number of equations to automatically convert 
  statements about roots of polynomials into a canonical form so that they 
  can be proven using the above substitutions.
›

definition "PR_TAG x  x"

lemma sturm_id_PR_prio0:
  "{x::real. P x} = {x::real. (PR_TAG P) x}"
  "(x::real. f x < g x) = (x::real. PR_TAG (λx. f x < g x) x)"
  "(x::real. P x) = (x::real. ¬(PR_TAG (λx. ¬P x)) x)"
  by (simp_all add: PR_TAG_def)

lemma sturm_id_PR_prio1:
  "{x::real. x < a  P x} = {x::real. x < a  (PR_TAG P) x}"
  "{x::real. x  a  P x} = {x::real. x  a  (PR_TAG P) x}"
  "{x::real. x  b  P x} = {x::real. x  b  (PR_TAG P) x}"
  "{x::real. x > b  P x} = {x::real. x > b  (PR_TAG P) x}"
  "(x::real < a. f x < g x) = (x::real < a. PR_TAG (λx. f x < g x) x)"
  "(x::real  a. f x < g x) = (x::real  a. PR_TAG (λx. f x < g x) x)"
  "(x::real > a. f x < g x) = (x::real > a. PR_TAG (λx. f x < g x) x)"
  "(x::real  a. f x < g x) = (x::real  a. PR_TAG (λx. f x < g x) x)"
  "(x::real < a. P x) = (x::real < a. ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real > a. P x) = (x::real > a. ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real  a. P x) = (x::real  a. ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real  a. P x) = (x::real  a. ¬(PR_TAG (λx. ¬P x)) x)"
  by (simp_all add: PR_TAG_def)

lemma sturm_id_PR_prio2:
  "{x::real. x > a  x  b  P x} = 
       {x::real. x > a  x  b  PR_TAG P x}"
  "{x::real. x  a  x  b  P x} = 
       {x::real. x  a  x  b  PR_TAG P x}"
  "{x::real. x  a  x < b  P x} = 
       {x::real. x  a  x < b  PR_TAG P x}"
  "{x::real. x > a  x < b  P x} = 
       {x::real. x > a  x < b  PR_TAG P x}"
  "(x::real. a < x  x  b  f x < g x) = 
       (x::real. a < x  x  b  PR_TAG (λx. f x < g x) x)"
  "(x::real. a  x  x  b  f x < g x) = 
       (x::real. a  x  x  b  PR_TAG (λx. f x < g x) x)"
  "(x::real. a < x  x < b  f x < g x) = 
       (x::real. a < x  x < b  PR_TAG (λx. f x < g x) x)"
  "(x::real. a  x  x < b  f x < g x) = 
       (x::real. a  x  x < b  PR_TAG (λx. f x < g x) x)"
  "(x::real. a < x  x  b  P x) = 
       (x::real. a < x  x  b  ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real. a  x  x  b  P x) = 
       (x::real. a  x  x  b  ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real. a  x  x < b  P x) = 
       (x::real. a  x  x < b  ¬(PR_TAG (λx. ¬P x)) x)"
  "(x::real. a < x  x < b  P x) = 
       (x::real. a < x  x < b  ¬(PR_TAG (λx. ¬P x)) x)"
  by (simp_all add: PR_TAG_def)



lemma PR_TAG_intro_prio0:
  fixes P :: "real  bool" and f :: "real  real"
  shows
  "PR_TAG P = P'  PR_TAG (λx. ¬(¬P x)) = P'"
  "PR_TAG P = (λx. poly p x = 0); PR_TAG Q = (λx. poly q x = 0)
        PR_TAG (λx. P x  Q x) = (λx. poly (gcd p q) x = 0)" and
 " PR_TAG P = (λx. poly p x = 0); PR_TAG Q = (λx. poly q x = 0)
        PR_TAG (λx. P x  Q x) = (λx. poly (p*q) x = 0)" and

  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x = g x) = (λx. poly (p-q) x = 0)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x  g x) = (λx. poly (p-q) x  0)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x < g x) = (λx. poly (q-p) x > 0)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x  g x) = (λx. poly (q-p) x  0)"

  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. -f x) = (λx. poly (-p) x)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x + g x) = (λx. poly (p+q) x)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x - g x) = (λx. poly (p-q) x)"
  "PR_TAG f = (λx. poly p x); PR_TAG g = (λx. poly q x)
        PR_TAG (λx. f x * g x) = (λx. poly (p*q) x)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. (f x)^n) = (λx. poly (p^n) x)"
  "PR_TAG (λx. poly p x :: real) = (λx. poly p x)"
  "PR_TAG (λx. x::real) = (λx. poly [:0,1:] x)"
  "PR_TAG (λx. a::real) = (λx. poly [:a:] x)"
  by (simp_all add: PR_TAG_def poly_eq_0_iff_dvd field_simps)


lemma PR_TAG_intro_prio1:
  fixes f :: "real  real"
  shows
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x = 0) = (λx. poly p x = 0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x  0) = (λx. poly p x  0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. 0 = f x) = (λx. poly p x = 0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. 0  f x) = (λx. poly p x  0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x  0) = (λx. poly p x  0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x > 0) = (λx. poly p x > 0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x  0) = (λx. poly (-p) x  0)"
  "PR_TAG f = (λx. poly p x)  PR_TAG (λx. f x < 0) = (λx. poly (-p) x > 0)"
  "PR_TAG f = (λx. poly p x)  
       PR_TAG (λx. 0  f x) = (λx. poly (-p) x  0)"
  "PR_TAG f = (λx. poly p x)  
       PR_TAG (λx. 0 < f x) = (λx. poly (-p) x < 0)"
  "PR_TAG f = (λx. poly p x) 
        PR_TAG (λx. a * f x) = (λx. poly (smult a p) x)"
  "PR_TAG f = (λx. poly p x) 
        PR_TAG (λx. f x * a) = (λx. poly (smult a p) x)"
  "PR_TAG f = (λx. poly p x) 
        PR_TAG (λx. f x / a) = (λx. poly (smult (inverse a) p) x)"
  "PR_TAG (λx. x^n :: real) = (λx. poly (monom 1 n) x)"
by (simp_all add: PR_TAG_def field_simps poly_monom)

lemma PR_TAG_intro_prio2:
  "PR_TAG (λx. 1 / b) = (λx. inverse b)"
  "PR_TAG (λx. a / b) = (λx. a / b)"
  "PR_TAG (λx. a / b * x^n :: real) = (λx. poly (monom (a/b) n) x)"
  "PR_TAG (λx. x^n * a / b :: real) = (λx. poly (monom (a/b) n) x)"
  "PR_TAG (λx. a * x^n :: real) = (λx. poly (monom a n) x)"
  "PR_TAG (λx. x^n * a :: real) = (λx. poly (monom a n) x)"
  "PR_TAG (λx. x^n / a :: real) = (λx. poly (monom (inverse a) n) x)"
(* TODO: can this be done more efficiently? I should think so. *)
  "PR_TAG (λx. f x^(Suc (Suc 0)) :: real) = (λx. poly p x)
        PR_TAG (λx. f x * f x :: real) = (λx. poly p x)"
  "PR_TAG (λx. (f x)^Suc n :: real) = (λx. poly p x)
        PR_TAG (λx. (f x)^n * f x :: real) = (λx. poly p x)"
  "PR_TAG (λx. (f x)^Suc n :: real) = (λx. poly p x)
        PR_TAG (λx. f x * (f x)^n :: real) = (λx. poly p x)"
  "PR_TAG (λx. (f x)^(m+n) :: real) = (λx. poly p x)
        PR_TAG (λx. (f x)^m * (f x)^n :: real) = (λx. poly p x)"
by (simp_all add: PR_TAG_def field_simps poly_monom power_add)

lemma sturm_meta_spec: "(x::real. P x)  P x" by simp
lemma sturm_imp_conv: 
  "(a < x  x < b  c)  (a < x  x < b  c)"
  "(a  x  x < b  c)  (a  x  x < b  c)"
  "(a < x  x  b  c)  (a < x  x  b  c)"
  "(a  x  x  b  c)  (a  x  x  b  c)"
  "(x < b  a < x  c)  (a < x  x < b  c)"
  "(x < b  a  x  c)  (a  x  x < b  c)"
  "(x  b  a < x  c)  (a < x  x  b  c)"
  "(x  b  a  x  c)  (a  x  x  b  c)"
  by auto


subsection ‹Setup for the ``sturm'' method›

ML_file ‹sturm.ML›

method_setup sturm = ‹
  Scan.succeed (fn ctxt => SIMPLE_METHOD' (Sturm.sturm_tac ctxt true))

end

File ‹sturm.ML›

signature STURM =
sig
  val sturm_conv : Proof.context -> conv
  val sturm_tac :