Theory Sturm_Sequences.Misc_Polynomial

(* Author: Manuel Eberl <manuel@pruvisto.org> *)
theory Misc_Polynomial
imports "HOL-Computational_Algebra.Polynomial" "HOL-Computational_Algebra.Polynomial_Factorial" "Pure-ex.Guess"
begin

subsection ‹Analysis›

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

subsection ‹Polynomials›

subsubsection ‹General simplification lemmas›

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


subsubsection ‹Divisibility of polynomials›

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

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

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

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

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

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

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

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

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

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

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

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

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



subsubsection ‹Sign changes of a polynomial›

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

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


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



subsubsection ‹Limits of polynomials›

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


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

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


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

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

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

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

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

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



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

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


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




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

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



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

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

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

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

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

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


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

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





subsubsection ‹Signs of polynomials for sufficiently large values›

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

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


subsubsection ‹Positivity of polynomials›

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

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

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

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

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

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

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

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

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

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

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

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

end