# 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 x⇩_{0})" proof- { fix ε :: real assume "ε > 0" have fin: "finite {x. ¦x-x⇩_{0}¦ < ε ∧ x ≠ x⇩_{0}∧ poly p x = 0}" using poly_roots_finite[OF assms] by simp with ‹ε > 0›have "∃δ>0. δ≤ε ∧ (∀x. ¦x-x⇩_{0}¦ < δ ∧ x ≠ x⇩_{0}⟶ poly p x ≠ 0)" proof (induction "card {x. ¦x-x⇩_{0}¦ < ε ∧ x ≠ x⇩_{0}∧ poly p x = 0}" arbitrary: ε rule: less_induct) case (less ε) let ?A = "{x. ¦x - x⇩_{0}¦ < ε ∧ x ≠ x⇩_{0}∧ 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 - x⇩_{0}¦ < ε ∧ x ≠ x⇩_{0}∧ poly p x = 0} ≠ {}" by force then obtain x where x_props: "¦x - x⇩_{0}¦ < ε" "x ≠ x⇩_{0}" "poly p x = 0" by blast define ε' where "ε' = ¦x - x⇩_{0}¦ / 2" have "ε' > 0" "ε' < ε" unfolding ε'_def using x_props by simp_all from x_props(1,2) and ‹ε > 0› have "x ∉ {x'. ¦x' - x⇩_{0}¦ < ε' ∧ x' ≠ x⇩_{0}∧ poly p x' = 0}" (is "_ ∉ ?B") by (auto simp: ε'_def) moreover from x_props have "x ∈ {x. ¦x - x⇩_{0}¦ < ε ∧ x ≠ x⇩_{0}∧ 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 (x⇩_{0}:: real) ≠ 0" shows "eventually (λx. sgn (poly p x) = sgn (poly p x⇩_{0})) (at x⇩_{0})" proof - have cont: "isCont (λx. sgn (poly p x)) x⇩_{0}" by (rule isCont_sgn, rule poly_isCont, rule assms) then have "eventually (λx. ¦sgn (poly p x) - sgn (poly p x⇩_{0})¦ < 1) (at x⇩_{0})" 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. x≤l ⟹ poly p x ≠ 0" by (fastforce simp: l_def) from Max_ge[OF A] have u_props: "⋀x. x≥u ⟹ 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 = (∑i≤n. 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 = (∑i≤n. 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. ∑i≤n. coeff p i / x^(n-i)) ⤏ (∑i≤n. ?a i)) at_infinity" by (force intro!: tendsto_sum) also have "(∑i≤n. ?a i) = coeff p n" by (subst sum.delta, simp_all) finally show "((λx. ∑i≤n. 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 x⇩_{0}where "⋀x. x ≥ x⇩_{0}⟹ poly p x ≥ 1" by (fastforce simp: filterlim_at_top eventually_at_top_linorder less_eq_real_def) hence "⋀x. x ≥ x⇩_{0}⟹ sgn (poly p x) = 1" by force thus ?thesis by (simp only: eventually_at_top_linorder poly_inf_def, intro exI[of _ x⇩_{0}], simp add: True) next case False hence "?lc < 0" using ‹?lc ≠ 0› by linarith from poly_at_bot_at_top[OF deg this] obtain x⇩_{0}where "⋀x. x ≥ x⇩_{0}⟹ poly p x ≤ -1" by (fastforce simp: filterlim_at_bot eventually_at_top_linorder less_eq_real_def) hence "⋀x. x ≥ x⇩_{0}⟹ sgn (poly p x) = -1" by force thus ?thesis by (simp only: eventually_at_top_linorder poly_inf_def, intro exI[of _ x⇩_{0}], 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 x⇩_{0}where "⋀x. x ≤ x⇩_{0}⟹ 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 ≤ x⇩_{0}⟹ sgn (poly p x) = 1" by force thus ?thesis by (simp add: True eventually_at_bot_linorder poly_neg_inf_def, intro exI[of _ x⇩_{0}], simp add: lc_pos) next case False from poly_at_top_or_bot_at_bot[OF deg lc_pos] and False obtain x⇩_{0}where "⋀x. x ≤ x⇩_{0}⟹ poly p x ≤ -1" by (fastforce simp add: filterlim_at_bot eventually_at_bot_linorder less_eq_real_def) hence "⋀x. x ≤ x⇩_{0}⟹ sgn (poly p x) = -1" by force thus ?thesis by (simp add: False eventually_at_bot_linorder poly_neg_inf_def, intro exI[of _ x⇩_{0}], 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 x⇩_{0}where "⋀x. x ≤ x⇩_{0}⟹ poly p x ≤ -1" by (fastforce simp: filterlim_at_bot eventually_at_bot_linorder less_eq_real_def) hence "⋀x. x ≤ x⇩_{0}⟹ sgn (poly p x) = -1" by force thus ?thesis by (simp only: True eventually_at_bot_linorder poly_neg_inf_def, intro exI[of _ x⇩_{0}], simp add: lc_neg) next case False with poly_at_bot_or_top_at_bot[OF deg lc_neg] obtain x⇩_{0}where "⋀x. x ≤ x⇩_{0}⟹ poly p x ≥ 1" by (fastforce simp: filterlim_at_top eventually_at_bot_linorder less_eq_real_def) hence "⋀x. x ≤ x⇩_{0}⟹ sgn (poly p x) = 1" by force thus ?thesis by (simp only: False eventually_at_bot_linorder poly_neg_inf_def, intro exI[of _ x⇩_{0}], 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: "∀x≥u'. 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: "∀x≤l'. 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 x⇩_{0}where "∀x≥x⇩_{0}. sgn (poly p x) = poly_inf p" by (auto simp: eventually_at_top_linorder) hence "poly_inf p = sgn (poly p (max x⇩_{0}(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 x⇩_{0}where C: "∀x≥x⇩_{0}. sgn (poly p x) = poly_inf p" by (auto simp: eventually_at_top_linorder) hence "sgn (poly p (max x⇩_{0}(a+1))) = poly_inf p" by simp with A have D: "sgn (poly p (max x⇩_{0}(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 x⇩_{0}(a+1)))" by (auto simp: sgn_real_def split: if_split_asm) show False apply (cases x' "max x⇩_{0}(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" "∀x≥a. 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 "∀x≥a. 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 x⇩_{0}where "∀x≤x⇩_{0}. sgn (poly p x) = poly_neg_inf p" by (auto simp: eventually_at_bot_linorder) hence "poly_neg_inf p = sgn (poly p (min x⇩_{0}(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 x⇩_{0}where C: "∀x≤x⇩_{0}. sgn (poly p x) = poly_neg_inf p" by (auto simp: eventually_at_bot_linorder) hence "sgn (poly p (min x⇩_{0}(a - 1))) = poly_neg_inf p" by simp with A have D: "sgn (poly p (min x⇩_{0}(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 x⇩_{0}(a - 1)))" by (auto simp: sgn_real_def split: if_split_asm) show False apply (cases x' "min x⇩_{0}(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" "∀x≤a. 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 "∀x≤a. 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 ∧ x≤b ⟶ poly p x ≠ 0)" hence A: "b ≤ a ∨ 0 < poly p b" 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 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. a≤x ∧ x<b ⟶ poly p x ≠ 0)" hence A: "b ≤ a ∨ 0 < poly p a" 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 > 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. a≤x ∧ x≤b ⟶ poly p x ≠ 0)" hence A: "b < a ∨ 0 < poly p a" 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 > 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 (ps⇩_{1}@ [p] @ ps⇩_{2}) x = sign_changes (ps⇩_{1}@ [p]) x + sign_changes ([p] @ ps⇩_{2}) 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: "⋀x⇩_{0}. poly p x⇩_{0}= 0 ⟹ eventually (λx. sgn (poly (p * ps!1) x) = (if x > x⇩_{0}then 1 else -1)) (at x⇩_{0})" 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 "∀j≤i+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 "∀j≤i+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) x⇩_{0}≠ 0" defines "sign_changes' ≡ λps x. ∑ps'←split_sign_changes ps x. sign_changes ps' x" shows "sign_changes' ps x⇩_{0}= sign_changes ps x⇩_{0}" using assms(1) proof (induction x⇩_{0}rule: split_sign_changes_induct) case (3 p q r ps x⇩_{0}) hence "poly p x⇩_{0}≠ 0" by simp note IH = 3(2,3,4) show ?case proof (cases "poly q x⇩_{0}= 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 x⇩_{0}* poly p x⇩_{0}< 0" by simp with 3 have "poly r x⇩_{0}≠ 0" by force from sign_changes_distrib[OF this, of "[p,q]" ps] have "sign_changes (p#q#r#ps) x⇩_{0}= sign_changes ([p, q, r]) x⇩_{0}+ sign_changes (r # ps) x⇩_{0}" by simp also have "sign_changes (r#ps) x⇩_{0}= sign_changes' (r#ps) x⇩_{0}" using ‹poly q x⇩_{0}= 0› ‹poly p x⇩_{0}≠ 0› 3(5)‹poly r x⇩_{0}≠ 0› by (intro IH(1)[symmetric], simp_all) finally show ?thesis unfolding sign_changes'_def using True ‹poly p x⇩_{0}≠ 0› by simp next case False from sign_changes_distrib[OF this, of "[p]" "r#ps"] have "sign_changes (p#q#r#ps) x⇩_{0}= sign_changes ([p,q]) x⇩_{0}+ sign_changes (q#r#ps) x⇩_{0}" by simp also have "sign_changes (q#r#ps) x⇩_{0}= sign_changes' (q#r#ps) x⇩_{0}" using ‹poly q x⇩_{0}≠ 0› ‹poly p x⇩_{0}≠ 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) x⇩_{0}≠ 0" defines "sign_changes' ≡ λx⇩_{0}ps x. ∑ps'←split_sign_changes ps x⇩_{0}. sign_changes ps' x" shows "eventually (λx. sign_changes' x⇩_{0}ps x = sign_changes ps x) (at x⇩_{0})" proof (rule eventually_mono) show "eventually (λx. ∀p∈{p ∈ set ps. poly p x⇩_{0}≠ 0}. sgn (poly p x) = sgn (poly p x⇩_{0})) (at x⇩_{0})" by (rule eventually_ball_finite, auto intro: poly_neighbourhood_same_sign) next fix x show "(∀p∈{p ∈ set ps. poly p x⇩_{0}≠ 0}. sgn (poly p x) = sgn (poly p x⇩_{0})) ⟹ sign_changes' x⇩_{0}ps x = sign_changes ps x" proof - fix x assume nbh: "∀p∈{p ∈ set ps. poly p x⇩_{0}≠ 0}. sgn (poly p x) = sgn (poly p x⇩_{0})" thus "sign_changes' x⇩_{0}ps x = sign_changes ps x" using assms(1) proof (induction x⇩_{0}rule: split_sign_changes_induct) case (3 p q r ps x⇩_{0}) hence "poly p x⇩_{0}≠ 0" by simp note IH = 3(2,3,4) show ?case proof (cases "poly q x⇩_{0}= 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 x⇩_{0}* poly p x⇩_{0}< 0" by simp with 3 have "poly r x⇩_{0}≠ 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' x⇩_{0}(r#ps) x" using ‹poly q x⇩_{0}= 0› nbh ‹poly p x⇩_{0}≠ 0› 3(5)‹poly r x⇩_{0}≠ 0› by (intro IH(1)[symmetric], simp_all) finally show ?thesis unfolding sign_changes'_def using True ‹poly p x⇩_{0}≠ 0›by 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' x⇩_{0}(q#r#ps) x" using ‹poly q x⇩_{0}≠ 0› nbh ‹poly p x⇩_{0}≠ 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) x⇩_{0}≠ 0" and "ps' ∈ set (split_sign_changes ps x⇩_{0})" shows "eventually (λx. sign_changes ps' x = sign_changes ps' x⇩_{0}) (at x⇩_{0})" using assms proof (induction x⇩_{0}rule: split_sign_changes_induct) case (1 p x) thus ?case by (simp add: sign_changes_def) next case (2 p q x⇩_{0}) hence [simp]: "ps' = [p,q]" by simp from 2 have "poly p x⇩_{0}≠ 0" by simp from 2(1) interpret quasi_sturm_seq "[p,q]" . from poly_neighbourhood_same_sign[OF ‹poly p x⇩_{0}≠ 0›] have "eventually (λx. sgn (poly p x) = sgn (poly p x⇩_{0})) (at x⇩_{0})" . moreover from last_ps_sgn_const have sgn_q: "⋀x. sgn (poly q x) = sgn (poly q x⇩_{0})" by simp ultimately have A: "eventually (λx. ∀p∈set[p,q]. sgn (poly p x) = sgn (poly p x⇩_{0})) (at x⇩_{0})" by simp thus ?case by (force intro: eventually_mono[OF A] sign_changes_cong') next case (3 p q r ps'' x⇩_{0}) hence p_not_0: "poly p x⇩_{0}≠ 0" by simp note sturm = 3(1) note IH = 3(2,3) note ps''_props = 3(6) show ?case proof (cases "poly q x⇩_{0}= 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 x⇩_{0}* poly p x⇩_{0}< 0" by simp with p_not_0 have r_not_0: "poly r x⇩_{0}≠ 0" by force show ?thesis proof (cases "ps' ∈ set (split_sign_changes (r # ps'') x⇩_{0})") 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 x⇩_{0}* poly p x⇩_{0}< 0" by simp from p_not_0 sgn_r have A: "eventually (λx. sgn (poly p x) = sgn (poly p x⇩_{0}) ∧ sgn (poly r x) = sgn (poly r x⇩_{0})) (at x⇩_{0})" 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 x⇩_{0})" and B: "sgn (poly r x) = sgn (poly r x⇩_{0})" 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 x⇩_{0}≠ 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' x⇩_{0}" 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'') x⇩_{0})") 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]: "∀p∈set ps'. poly p x⇩_{0}≠ 0" using q_not_0 p_not_0 by simp show ?thesis proof (rule eventually_mono) fix x assume "∀p∈set ps'. sgn (poly p x) = sgn (poly p x⇩_{0})" thus "sign_changes ps' x = sign_changes ps' x⇩_{0}" by (rule sign_changes_cong') next show "eventually (λx. ∀p∈set ps'. sgn (poly p x) = sgn (poly p x⇩_{0})) (at x⇩_{0})" 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) x⇩_{0}≠ 0" shows "eventually (λx. sign_changes ps x = sign_changes ps x⇩_{0}) (at x⇩_{0})" proof- let ?pss = "split_sign_changes ps x⇩_{0}" 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' x⇩_{0}) (at x⇩_{0})" hence "eventually (λx. ?f pss x = ?f pss x⇩_{0}) (at x⇩_{0})" 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 x⇩_{0}) (at x⇩_{0})" 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 x⇩_{0}≠ 0 ⟹ eventually (λx. sign_changes ps x = sign_changes ps x⇩_{0}) (at x⇩_{0})" 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 x⇩_{0}= 0" "p ≠ 0" shows "eventually (λx. sign_changes ps x = sign_changes ps x⇩_{0}+ (if x<x⇩_{0}then 1 else 0)) (at x⇩_{0})" proof- from ps_first_two obtain q ps' where [simp]: "ps = p#q#ps'" . hence "ps!1 = q" by simp have "eventually (λx. x ≠ x⇩_{0}) (at x⇩_{0})" by (simp add: eventually_at, rule exI[of _ 1], simp) moreover from p_squarefree and assms(1) have "poly q x⇩_{0}≠ 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 x⇩_{0}≠ 0› have "eventually (λx. sign_changes (q#ps') x = sign_changes (q#ps') x⇩_{0}) (at x⇩_{0})" using hd_nonzero_imp_sign_changes_const[where x⇩_{0}=x⇩_{0}] by simp } moreover note poly_neighbourhood_without_roots[OF assms(2)] deriv[OF assms(1)] ultimately have A: "eventually (λx. x ≠ x⇩_{0}∧ poly p x ≠ 0 ∧ sgn (poly (p*ps!1) x) = (if x > x⇩_{0}then 1 else -1) ∧ sign_changes (q#ps') x = sign_changes (q#ps') x⇩_{0}) (at x⇩_{0})" 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 < x⇩_{0}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<x⇩_{0}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') x⇩_{0}= sign_changes ps x⇩_{0}" 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 x⇩_{2}where "x⇩_{2}= Min {x. x > a ∧ x ≤ b ∧ poly p x = 0}" from Min_in[OF fin] and True have x⇩_{2}_props: "x⇩_{2}> a" "x⇩_{2}≤ b" "poly p x⇩_{2}= 0" unfolding x⇩_{2}_def by blast+ from Min_le[OF fin] x⇩_{2}_props have x⇩_{2}_le: "⋀x'. ⟦x' > a; x' ≤ b; poly p x' = 0⟧ ⟹ x⇩_{2}≤ x'" unfolding x⇩_{2}_def by simp have left: "{x. a < x ∧ x ≤ x⇩_{2}∧ poly p x = 0} = {x⇩_{2}}" using x⇩_{2}_props x⇩_{2}_le by force hence [simp]: "card {x. a < x ∧ x ≤ x⇩_{2}∧ poly p x = 0} = 1" by simp from p_zero[OF ‹poly p x⇩_{2}= 0› ‹p ≠ 0›, unfolded eventually_at dist_real_def] guess ε .. hence ε_props: "ε > 0" "∀x. x ≠ x⇩_{2}∧ ¦x - x⇩_{2}¦ < ε ⟶ sign_changes ps x = sign_changes ps x⇩_{2}+ (if x < x⇩_{2}then 1 else 0)" by auto define x⇩_{1}where "x⇩_{1}= max (x⇩_{2}- ε / 2) a" have "¦x⇩_{1}- x⇩_{2}¦ < ε" using ‹ε > 0› x⇩_{2}_props by (simp add: x⇩_{1}_def) hence "sign_changes ps x⇩_{1}= (if x⇩_{1}< x⇩_{2}then sign_changes ps x⇩_{2}+ 1 else sign_changes ps x⇩_{2})" using ε_props(2) by (cases "x⇩_{1}= x⇩_{2}", auto) hence "sign_changes ps x⇩_{1}- sign_changes ps x⇩_{2}= 1" unfolding x⇩_{1}_def using x⇩_{2}_props ‹ε > 0› by simp also have "x⇩_{2}∉ {x. a < x ∧ x ≤ x⇩_{1}∧ poly p x = 0}" unfolding x⇩_{1}_def using ‹ε > 0› by force with left have "{x. a < x ∧ x ≤ x⇩_{1}∧ poly p x = 0} = {}" by force with less(1)[of a x⇩_{1}] have "sign_changes ps x⇩_{1}= sign_changes ps a" unfolding x⇩_{1}_def ‹ε > 0› by (force simp: card_greater_0) finally have signs_left: "sign_changes ps a - int (sign_changes ps x⇩_{2}) = 1" by simp have "{x. x > a ∧ x ≤ b ∧ poly p x = 0} = {x. a < x ∧ x ≤ x⇩_{2}∧ poly p x = 0} ∪ {x. x⇩_{2}< x ∧ x ≤ b ∧ poly p x = 0}" using x⇩_{2}_props by auto also note left finally have A: "card {x. x⇩_{2}< 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. x⇩_{2}< x ∧ x ≤ b ∧ poly p x = 0} < card {x. a < x ∧ x ≤ b ∧ poly p x = 0}" by simp from less(1)[OF this x⇩_{2}_props(2)] and A have signs_right: "sign_changes ps x⇩_{2}- 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 "∀j≤i+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 x⇩_{0})" assumes p_0: "poly p (x⇩_{0}::real) = 0" shows "eventually (λx. sgn (poly (p*q) x) = (if x > x⇩_{0}then 1 else -1)) (at x⇩_{0})" proof- have A: "eventually (λx. poly p x ≠ 0 ∧ poly q x ≠ 0 ∧ sgn (poly q x) = sgn (poly (pderiv p) x)) (at x⇩_{0})" 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 ≠ x⇩_{0}∧ ¦x - x⇩_{0}¦ < ε ⟶ 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 ≠ x⇩_{0}" "¦x - x⇩_{0}¦ < ε" 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 > x⇩_{0}then 1 else -1)" proof (cases "x ≥ x⇩_{0}") case True with ‹x ≠ x⇩_{0}› have "x > x⇩_{0}" by simp from poly_MVT[OF this, of p] guess ξ .. note ξ_props = this with ‹¦x - x⇩_{0}¦ < ε› ‹poly p x⇩_{0}= 0› ‹x > x⇩_{0}› ε_props have "¦ξ - x⇩_{0}¦ < ε" "sgn (poly p x) = sgn (x - x⇩_{0}) * sgn (poly q ξ)" by (auto simp add: q_pderiv sgn_mult) moreover from ξ_props ε_props ‹¦x - x⇩_{0}¦ < ε› 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 ≠ x⇩_{0}› ε_props ξ_props by (auto simp: sgn_mult sqr_pos) next case False hence "x < x⇩_{0}" by simp hence sgn: "sgn (x - x⇩_{0}) = -1" by simp from poly_MVT[OF ‹x < x⇩_{0}›, of p] guess ξ .. note ξ_props = this with ‹¦x - x⇩_{0}¦ < ε› ‹poly p x⇩_{0}= 0› ‹x < x⇩_{0}› ε_props have "¦ξ - x⇩_{0}¦ < ε" "poly p x = (x - x⇩_{0}) * poly (pderiv p) ξ" "poly p ξ ≠ 0" by (auto simp: field_simps) hence "sgn (poly p x) = sgn (x - x⇩_{0}) * sgn (poly q ξ)" using ε_props ξ_props by (auto simp: q_pderiv sgn_mult) moreover from ξ_props ε_props ‹¦x - x⇩_{0}¦ < ε› 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 ≠ x⇩_{0}› 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 (x⇩_{0}::real) = 0" shows "eventually (λx. sgn (poly (p * sturm p ! 1) x) = (if x > x⇩_{0}then 1 else -1)) (at x⇩_{0})" 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 "⋀x⇩_{0}. poly p x⇩_{0}= 0 ⟹ eventually (λx. sgn (poly (p*sturm p ! 1) x) = (if x > x⇩_{0}then 1 else -1)) (at x⇩_{0})" 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 "⋀x⇩_{0}. poly p' x⇩_{0}= 0 ⟹ eventually (λx. sgn (poly (p' * sturm_squarefree p ! 1) x) = (if x > x⇩_{0}then 1 else -1)) (at x⇩_{0})" . 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 "∀j≤i+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 x⇩_{0}:: real assume "poly ?p' x⇩_{0}= 0" hence "poly p x⇩_{0}= 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 x⇩_{0}= 0› have A: "eventually (λx. sgn (poly (p * pderiv p) x) = (if x⇩_{0}< x then 1 else -1)) (at x⇩_{0})" 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 x⇩_{0}< x then 1 else -1)) (at x⇩_{0})" 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 [x←map f (map (λq. q div d) ps). x ≠ 0]) - 1 = length (remdups_adj [x←map (λ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 [x←map f ps' . x ≠ 0]) - 1 = length (remdups_adj [x←map 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 :