Session Three_Circles

Theory RRI_Misc

section ‹Misc results about polynomials›

theory RRI_Misc imports 
  "HOL-Computational_Algebra.Computational_Algebra" 
  "Budan_Fourier.BF_Misc"
  "Polynomial_Interpolation.Ring_Hom_Poly"
begin

subsection ‹Misc›

declare pcompose_pCons[simp del]

lemma Setcompr_subset: "f P S. {f x | x. P x}  S = ( x. P x  f x  S)" 
  by blast

lemma map_cong':
  assumes "xs = map h ys" and "y. y  set ys  f (h y) = g y"
  shows "map f xs = map g ys"
  using assms map_replicate_trivial by simp

lemma nth_default_replicate_eq: 
    "nth_default dflt (replicate n x) i = (if i < n then x else dflt)"
  by (auto simp: nth_default_def)

lemma square_bounded_less: 
  fixes a b::"'a :: linordered_ring_strict"
  shows "-a < b  b < a  b*b < a*a"
  by (metis (no_types, lifting) leD leI minus_less_iff minus_mult_minus mult_strict_mono'
      neg_less_eq_nonneg neg_less_pos verit_minus_simplify(4) zero_le_mult_iff zero_le_square)

lemma square_bounded_le: 
  fixes a b::"'a :: linordered_ring_strict"
  shows "-a  b  b  a  b*b  a*a"
  by (metis le_less minus_mult_minus square_bounded_less)

context vector_space
begin

lemma card_le_dim_spanning:
  assumes BV: "B  V"
    and VB: "V  span B"
    and fB: "finite B"
    and dVB: "dim V  card B"
  shows "independent B"
proof -
  {
    fix a
    assume a: "a  B" "a  span (B - {a})"
    from a fB have c0: "card B  0"
      by auto
    from a fB have cb: "card (B - {a}) = card B - 1"
      by auto
    {
      fix x
      assume x: "x  V"
      from a have eq: "insert a (B - {a}) = B"
        by blast
      from x VB have x': "x  span B"
        by blast
      from span_trans[OF a(2), unfolded eq, OF x']
      have "x  span (B - {a})" .
    }
    then have th1: "V  span (B - {a})"
      by blast
    have th2: "finite (B - {a})"
      using fB by auto
    from dim_le_card[OF th1 th2]
    have c: "dim V  card (B - {a})" .
    from c c0 dVB cb have False by simp
  }
  then show ?thesis
    unfolding dependent_def by blast
qed

end

subsection ‹Misc results about polynomials›

lemma smult_power: "smult (x^n) (p^n) = smult x p ^ n"
  apply (induction n)
  subgoal by fastforce
  by (metis (no_types, hide_lams) mult_smult_left mult_smult_right 
      power_Suc smult_smult)

lemma reflect_poly_monom: "reflect_poly (monom n i) = monom n 0"
  apply (induction i)
  by (auto simp: coeffs_eq_iff coeffs_monom coeffs_reflect_poly)

lemma poly_eq_by_eval: 
  fixes P Q :: "'a::{comm_ring_1,ring_no_zero_divisors,ring_char_0} poly"
  assumes h: "x. poly P x = poly Q x" shows "P = Q"
proof -
  have "poly P = poly Q" using h by fast
  thus ?thesis by (auto simp: poly_eq_poly_eq_iff)
qed

lemma poly_binomial: 
  "[:(1::'a::comm_ring_1), 1:]^n = (kn. monom (of_nat (n choose k)) k)"
proof -
  have "[:(1::'a::comm_ring_1), 1:]^n = (monom 1 1 + 1)^n"
    by (metis (no_types, lifting) add.left_neutral add.right_neutral add_pCons
        monom_altdef pCons_one power_one_right smult_1_left)
  also have "... = (kn. of_nat (n choose k) * monom 1 1 ^ k)"
    apply (subst binomial_ring)
    by force
  also have "... = (kn. monom (of_nat (n choose k)) k)"
    by (auto simp: monom_altdef of_nat_poly)
  finally show ?thesis .
qed

lemma degree_0_iff: "degree P = 0  (a. P = [:a:])"
  by (meson degree_eq_zeroE degree_pCons_0)

interpretation poly_vs: vector_space smult
  by (simp add: vector_space_def smult_add_right smult_add_left)

lemma degree_subspace: "poly_vs.subspace {x. degree x  n}"
  by (auto simp: poly_vs.subspace_def degree_add_le)

lemma monom_span: 
  "poly_vs.span {monom 1 x | x. x  p} = {(x::'a::field poly). degree x  p}"
(is "?L = ?R")
proof
  show "?L  ?R"
  proof
    fix x assume "x  ?L"
    moreover have hfin: "finite {P. x  {..p}. P = monom 1 x}"
      by auto
    ultimately have 
      "x  range (λu. v{monom 1 x | x. x  {..p}}. smult (u v) v)"
      by (simp add: poly_vs.span_finite)
    hence " u. x = (v{monom 1 x | x. x  {..p}}. smult (u v) v)"
      by (auto simp: image_iff)
    then obtain u 
      where p': "x = (v{monom 1 x | x. x  {..p}}. smult (u v) v)"
      by blast
    have "v. v  {monom 1 x | x. x  {..p}}  degree (smult (u v) v)  p"
      by (auto simp add: degree_monom_eq)
    hence "degree x  p" using hfin
      apply (subst p')
      apply (rule degree_sum_le)
      by auto
    thus "x  {x. degree x  p}" by force
  qed
next
  show "?R  ?L"
  proof
    fix x assume "x  ?R"
    hence "degree x  p" by force
    hence "x = (ip. monom (coeff x i) i)"
      by (simp add: poly_as_sum_of_monoms')
    also have
      "... = (ip. smult (coeff x (degree (monom (1::'a) i))) (monom 1 i))"
      by (auto simp add: smult_monom degree_monom_eq)
    also have
      "... = (v{monom 1 x | x. x  {..p}}. smult ((λv. coeff x (degree v)) v) v)"
    proof (rule sum.reindex_cong)
      show "inj_on degree {monom (1::'a) x |x. x  {..p}}"
      proof
        fix x 
        assume "x  {monom (1::'a) x |x. x  {..p}}" 
        hence " a. x = monom 1 a" by blast
        then obtain a where hx: "x = monom 1 a" by blast
        fix y 
        assume "y  {monom (1::'a) x |x. x  {..p}}" 
        hence " b. y = monom 1 b" by blast
        then obtain b where hy: "y = monom 1 b" by blast
        assume "degree x = degree y"
        thus "x = y" using hx hy by (simp add: degree_monom_eq)
      qed
      show "{..p} = degree ` {monom (1::'a) x |x. x  {..p}}"
      proof
        show "{..p}  degree ` {monom (1::'a) x |x. x  {..p}}"
        proof
          fix x assume "x  {..p}"
          hence "monom (1::'a) x  {monom 1 x |x. x  {..p}}" by force
          moreover have "degree (monom (1::'a) x) = x"
            by (simp add: degree_monom_eq)
          ultimately show "x  degree ` {monom (1::'a) x |x. x  {..p}}" by auto
        qed
        show "degree ` {monom (1::'a) x |x. x  {..p}}  {..p}"
          by (auto simp add: degree_monom_eq)
      qed
    next
      fix y assume "y  {monom (1::'a) x |x. x  {..p}}"
      hence "z  {..p}. y = monom (1::'a) z" by blast
      then obtain z where "y = monom (1::'a) z" by blast
      thus 
        "smult (coeff x (degree (monom (1::'a) (degree y)))) (monom (1::'a) (degree y)) =
         smult (coeff x (degree y)) y"
        by (simp add: smult_monom degree_monom_eq)
    qed
    finally have "x = (v{monom 1 x | x. x  {..p}}. 
                      smult ((λv. coeff x (degree v)) v) v)" .
    thus "x  ?L" by (auto simp add: poly_vs.span_finite)
  qed
qed

lemma monom_independent: 
  "poly_vs.independent {monom (1::'a::field) x | x. x  p}"
proof (rule poly_vs.independent_if_scalars_zero)
  fix f::"'a poly  'a"
  assume h: "(x{monom 1 x |x. x  p}. smult (f x) x) = 0"
  have h': "(ip. monom (f (monom (1::'a) i)) i) = 
            (x{monom (1::'a) x |x. x  p}. smult (f x) x)"
  proof (rule sum.reindex_cong)
    show "inj_on degree {monom (1::'a) x |x. x  p}"
      by (smt (verit) degree_monom_eq inj_on_def mem_Collect_eq zero_neq_one)
    show "{..p} = degree ` {monom (1::'a) x |x. x  p}"
    proof
      show "{..p}  degree ` {monom (1::'a) x |x. x  p}"
      proof
        fix x assume "x  {..p}"
        then have "x = degree (monom (1::'a) x)  x  p"
          by (auto simp: degree_monom_eq)
        thus "x  degree ` {monom (1::'a) x |x. x  p}" 
          by blast
      qed
      show "degree ` {monom (1::'a) x |x. x  p}  {..p}"
        by (force simp: degree_monom_eq)
    qed
  qed (auto simp: degree_monom_eq smult_monom)

  fix x::"'a poly"
  assume "x  {monom 1 x |x. x  p}"
  then obtain y where "y  p" and "x = monom 1 y" by blast
  hence "f x = coeff (x{monom 1 x |x. x  p}. smult (f x) x) y"
    by (auto simp: coeff_sum h'[symmetric])
  thus "f x = 0"
    using h by auto
qed force

lemma dim_degree: "poly_vs.dim {x. degree x  n} = n + 1"
  using poly_vs.dim_eq_card_independent[OF monom_independent]
  by (auto simp: monom_span[symmetric] card_image image_Collect[symmetric]
      inj_on_def monom_eq_iff')

lemma degree_div: 
  fixes p q::"('a::idom_divide) poly" 
  assumes "q dvd p"
  shows "degree (p div q) = degree p - degree q" using assms
  by (metis (no_types, lifting) add_diff_cancel_left' degree_0 degree_mult_eq 
      diff_add_zero diff_zero div_by_0 dvd_div_eq_0_iff dvd_mult_div_cancel)

lemma lead_coeff_div: 
  fixes p q::"('a::{idom_divide, inverse}) poly" 
  assumes "q dvd p"
  shows "lead_coeff (p div q) = lead_coeff p / lead_coeff q" using assms
  by (smt (z3) div_by_0 dvd_div_mult_self lead_coeff_mult leading_coeff_0_iff
      nonzero_mult_div_cancel_right)

lemma complex_poly_eq: 
  "r = map_poly of_real (map_poly Re r) + smult 𝗂 (map_poly of_real (map_poly Im r))"
  by (auto simp: poly_eq_iff coeff_map_poly complex_eq)

lemma complex_poly_cong: 
  "(map_poly Re p = map_poly Re q  map_poly Im p = map_poly Im q) = (p = q)" 
  by (metis complex_poly_eq)

lemma map_poly_Im_of_real: "map_poly Im (map_poly of_real p) = 0"
  by (auto simp: poly_eq_iff coeff_map_poly)

lemma mult_map_poly_imp_map_poly: 
  assumes "map_poly complex_of_real q = r * map_poly complex_of_real p" 
          "p  0"
  shows "r = map_poly complex_of_real (map_poly Re r)"
proof -
  have h: "Im  (*) 𝗂  complex_of_real = id" by fastforce
  have "map_poly complex_of_real q = r * map_poly complex_of_real p" 
    using assms by blast
  also have "... = (map_poly of_real (map_poly Re r) + 
                    smult 𝗂 (map_poly of_real (map_poly Im r))) *
                   map_poly complex_of_real p"
    using complex_poly_eq by fastforce
  also have "... = map_poly of_real (map_poly Re r * p) + 
                   smult 𝗂 (map_poly of_real (map_poly Im r * p))"
    by (simp add: mult_poly_add_left)
  finally have "map_poly complex_of_real q = 
                map_poly of_real (map_poly Re r * p) + 
                smult 𝗂 (map_poly of_real (map_poly Im r * p))" .
  hence "0 = map_poly Im (map_poly of_real (map_poly Re r * p) + 
             smult 𝗂 (map_poly of_real (map_poly Im r * p)))"
    by (auto simp: complex_poly_cong[symmetric] map_poly_Im_of_real)
  also have "... = map_poly of_real (map_poly Im r * p)"
    apply (rule poly_eqI)
    by (auto simp: coeff_map_poly coeff_mult)
  finally have "0 = map_poly complex_of_real (map_poly Im r) *
                    map_poly complex_of_real p" 
    by auto
  hence "map_poly complex_of_real (map_poly Im r) = 0" using assms by fastforce
  thus ?thesis apply (subst complex_poly_eq) by auto
qed

lemma map_poly_dvd: 
  fixes p q::"real poly"
  assumes hdvd: "map_poly complex_of_real p dvd 
                    map_poly complex_of_real q" "q  0"
  shows "p dvd q" 
proof -
  from hdvd obtain r 
    where h:"map_poly complex_of_real q = r * map_poly complex_of_real p"
    by fastforce
  hence "r = map_poly complex_of_real (map_poly Re r)" 
    using mult_map_poly_imp_map_poly assms by force
  hence "map_poly complex_of_real q = map_poly complex_of_real (p * map_poly Re r)" 
    using h by auto
  hence "q = p * map_poly Re r" using of_real_poly_eq_iff by blast
  thus "p dvd q" by force
qed

lemma div_poly_eq_0: 
  fixes p q::"('a::idom_divide) poly" 
  assumes "q dvd p" "poly (p div q) x = 0" "q  0"
  shows "poly p x = 0"
  using assms by fastforce

lemma poly_map_poly_of_real_cnj: 
    "poly (map_poly of_real p) (cnj z) = cnj (poly (map_poly of_real p) z)"
  apply (induction p)
  by auto

text ‹
An induction rule on real polynomials, if $P \neq 0$ then either $(X-x)|P$ or 
$(X-z)(X-cnj z)|P$, we induct by dividing by these polynomials.
›
lemma real_poly_roots_induct: 
  fixes P::"real poly  bool" and p::"real poly"
  assumes IH_real: "p x. P p  P (p * [:-x, 1:])"
      and IH_complex: "p a b. b  0  P p 
               P (p * [: a*a + b*b, -2*a, 1 :])"
      and H0: "a. P [:a:]"
  defines "d  degree p"
  shows "P p"
  using d_def
proof (induction d arbitrary: p rule: less_induct)
  fix p::"real poly"
  assume IH: "(q. degree q < degree p  P q)"
  show "P p"
  proof (cases "0 = degree p")
    fix p::"real poly" assume "0 = degree p"
    hence " a. p = [:a:]"
      by (simp add: degree_0_iff)
    thus "P p" using H0 by blast
  next
    assume hdeg: "0  degree p" 
    hence "¬ constant (poly (map_poly of_real p))"
      by (metis (no_types, hide_lams) constant_def constant_degree
          of_real_eq_iff of_real_poly_map_poly)
    then obtain z::complex where h: "poly (map_poly of_real p) z = 0" 
      using fundamental_theorem_of_algebra by blast
    show "P p"
    proof cases
      assume "Im z = 0"
      hence "z = Re z" by (simp add: complex_is_Real_iff)
      moreover have "[:-z, 1:] dvd map_poly of_real p" 
        using h poly_eq_0_iff_dvd by blast
      ultimately have "[:-(Re z), 1:] dvd p"
        by (smt (z3) dvd_iff_poly_eq_0 h of_real_0 of_real_eq_iff of_real_poly_map_poly)
      hence 2:"P (p div [:-Re z, 1:])"
        apply (subst IH)
        using hdeg by (auto simp: degree_div)
      moreover have 1:"p = (p div [:- Re z, 1:]) * [:-Re z, 1:]"
        by (metis [:- Re z, 1:] dvd p dvd_div_mult_self)
      ultimately show "P p"
        apply (subst 1)
        by (rule IH_real[OF 2])
    next
      assume "Im z  0"
      hence hcnj: "cnj z  z" by (metis cnj.simps(2) neg_equal_zero)
      have h2: "poly (map_poly of_real p) (cnj z) = 0" 
        using h poly_map_poly_of_real_cnj by force
      have "[:-z, 1:] * [:-cnj z, 1:] dvd map_poly of_real p"
      proof (rule divides_mult)
        have "c. c dvd [:-z, 1:]  c dvd [:- cnj z, 1:]  is_unit c"
        proof -
          fix c
          assume h:"c dvd [:-z, 1:]" hence "degree c  1" using divides_degree by fastforce
          hence "degree c = 0  degree c = 1" by linarith
          thus "c dvd [:- cnj z, 1:]  is_unit c"
          proof
            assume "degree c = 0"
            moreover have "c  0" using h by fastforce
            ultimately show "is_unit c" by (simp add: is_unit_iff_degree)
          next
            assume hdeg: "degree c = 1"
            then obtain x where 1:"[:-z, 1:] = x*c" using h by fastforce
            hence "degree [:-z, 1:] = degree x + degree c"
              by (metis add.inverse_neutral degree_mult_eq mult_cancel_right
                  mult_poly_0_left pCons_eq_0_iff zero_neq_neg_one)
            hence "degree x = 0" using hdeg by auto
            then obtain x' where 2: "x = [:x':]" using degree_0_iff by auto
            assume "c dvd [:-cnj z, 1:]"
            then obtain y where 3: "[:-cnj z, 1:] = y*c" by fastforce
            hence "degree [:-cnj z, 1:] = degree y + degree c"
              by (metis add.inverse_neutral degree_mult_eq mult_cancel_right
                  mult_poly_0_left pCons_eq_0_iff zero_neq_neg_one)
            hence "degree y = 0" using hdeg by auto
            then obtain y' where 4: "y = [:y':]" using degree_0_iff by auto
            moreover from hdeg obtain a b where 5:"c = [:a, b:]" and 6: "b  0"
              by (meson degree_eq_oneE)
            from 1 2 5 6 have "x' = inverse b" by (auto simp: field_simps)
            moreover from 3 4 5 6 have "y' = inverse b" by (auto simp: field_simps)
            ultimately have "x = y" using 2 4 by presburger
            then have "z = cnj z" using 1 3 by (metis neg_equal_iff_equal pCons_eq_iff)
            thus "is_unit c" using hcnj by argo
          qed
        qed
        thus "coprime [:- z, 1:] [:- cnj z, 1:]" by (meson not_coprimeE)
        show "[:- z, 1:] dvd map_poly complex_of_real p"
          using h poly_eq_0_iff_dvd by auto
        show "[:- cnj z, 1:] dvd map_poly complex_of_real p"
          using h2 poly_eq_0_iff_dvd by blast
      qed
      moreover have "[:- z, 1:] * [:- cnj z, 1:] = 
                     map_poly of_real [:Re z*Re z + Im z*Im z, -2*Re z, 1:]"
        by (auto simp: complex_eqI)
      ultimately have hdvd: 
        "map_poly complex_of_real [:Re z*Re z + Im z*Im z, -2*Re z, 1:] dvd
         map_poly complex_of_real p"
        by force
      hence "[:Re z*Re z + Im z*Im z, -2*Re z, 1:] dvd p" using map_poly_dvd 
        by blast
      hence 2:"P (p div [:Re z*Re z + Im z*Im z, -2*Re z, 1:])"
        apply (subst IH)
        using hdeg by (auto simp: degree_div)
      moreover have 1:
        "p = (p div [:Re z*Re z + Im z*Im z, -2*Re z, 1:]) * 
                    [:Re z*Re z + Im z*Im z, -2*Re z, 1:]"
        apply (subst dvd_div_mult_self) 
        using [:Re z*Re z + Im z*Im z, -2*Re z, 1:] dvd p by auto
      ultimately show "P p"
        apply (subst 1)
        apply (rule IH_complex[of  "Im z" _ "Re z"])
        apply (meson ‹Im z  0)
        by blast
    qed
  qed
qed

subsection ‹The reciprocal polynomial›

definition reciprocal_poly :: "nat  'a::zero poly  'a poly"
  where "reciprocal_poly p P = 
          Poly (rev ((coeffs P) @ (replicate (p - degree P) 0)))"

lemma reciprocal_0: "reciprocal_poly p 0 = 0" by (simp add: reciprocal_poly_def)

lemma reciprocal_1: "reciprocal_poly p 1 = monom 1 p"
  by (simp add: reciprocal_poly_def monom_altdef Poly_append)

lemma coeff_reciprocal: 
  assumes hi: "i  p" and hP: "degree P  p"
  shows "coeff (reciprocal_poly p P) i = coeff P (p - i)"
proof cases
  assume "i < p - degree P"
  hence "degree P < p - i" using hP by linarith
  thus "coeff (reciprocal_poly p P) i = coeff P (p - i)"
    by (auto simp: reciprocal_poly_def nth_default_append coeff_eq_0)
next
  assume h: "¬i < p - degree P"
  show "coeff (reciprocal_poly p P) i = coeff P (p - i)"
  proof cases
    assume "P = 0"
    thus "coeff (reciprocal_poly p P) i = coeff P (p - i)" 
      by (simp add: reciprocal_0)
  next
    assume hP': "P  0"
    have "degree P  p - i" using h hP by linarith
    moreover hence "(i - (p - degree P)) < length (rev (coeffs P))" 
      using hP' hP hi by (auto simp: length_coeffs)
    thus "coeff (reciprocal_poly p P) i = coeff P (p - i)"
      by (auto simp: reciprocal_poly_def nth_default_append coeff_eq_0 hP hP' 
          nth_default_nth rev_nth calculation coeffs_nth length_coeffs_degree)
  qed
qed

lemma coeff_reciprocal_less: 
  assumes hn: "p < i" and hP: "degree P  p"
  shows "coeff (reciprocal_poly p P) i = 0"
proof cases
  assume "P = 0"
  thus ?thesis by (auto simp: reciprocal_0)
next
  assume "P  0"
  thus ?thesis
    using hn 
    by (auto simp: reciprocal_poly_def nth_default_append 
        nth_default_eq_dflt_iff hP length_coeffs)
qed

lemma reciprocal_monom: 
  assumes "n  p"
  shows "reciprocal_poly p (monom a n) = monom a (p-n)"
proof (cases "a=0")
  case True
  then show ?thesis by (simp add: reciprocal_0)
next
  case False
  with np show ?thesis 
    apply (rule_tac poly_eqI)
    by (metis coeff_monom coeff_reciprocal coeff_reciprocal_less 
        diff_diff_cancel diff_le_self lead_coeff_monom not_le_imp_less)
qed

lemma reciprocal_degree: "reciprocal_poly (degree P) P = reflect_poly P"
  by (auto simp add: reciprocal_poly_def reflect_poly_def)

lemma degree_reciprocal:
  fixes P :: "('a::zero) poly" 
  assumes hP: "degree P  p"
  shows "degree (reciprocal_poly p P)  p"
proof (auto simp add: reciprocal_poly_def)
  have "degree (reciprocal_poly p P)  
        length (replicate (p - degree P) (0::'a) @ rev (coeffs P))"
    by (metis degree_Poly reciprocal_poly_def rev_append rev_replicate)
  thus "degree (Poly (replicate (p - degree P) 0 @ rev (coeffs P)))  p" 
    by (smt Suc_le_mono add_Suc_right coeffs_Poly degree_0 hP le_SucE le_SucI 
        le_add_diff_inverse2 le_zero_eq length_append length_coeffs_degree
        length_replicate length_rev length_strip_while_le reciprocal_0
        reciprocal_poly_def rev_append rev_replicate)
qed

lemma reciprocal_0_iff: 
  assumes hP: "degree P  p" 
  shows "(reciprocal_poly p P = 0) = (P = 0)"
proof
  assume h: "reciprocal_poly p P = 0"
  show "P = 0"
  proof (rule poly_eqI)
    fix n
    show "coeff P n = coeff 0 n"
    proof cases
      assume hn: "n  p"
      hence "p - n  p" by auto
      hence "coeff (reciprocal_poly p P) (p - n) = coeff P n" 
        using hP hn by (auto simp: coeff_reciprocal)
      thus ?thesis using h by auto
    next
      assume hn: "¬ n  p"
      thus ?thesis using hP by (metis coeff_0 dual_order.trans le_degree)
    qed
  qed
next
  assume "P = 0"
  thus "reciprocal_poly p P = 0" using reciprocal_0 by fast
qed

lemma poly_reciprocal: 
  fixes P::"'a::field poly"
  assumes hp: "degree P  p" and hx: "x  0"
  shows "poly (reciprocal_poly p P) x = x^p * (poly P (inverse x))"
proof -
  have "poly (reciprocal_poly p P) x
      = poly ((Poly ((replicate (p - degree P) 0) @ rev (coeffs P)))) x"
    by (auto simp add: hx reflect_poly_def reciprocal_poly_def)
  also have "... = poly ((monom 1 (p - degree P)) * (reflect_poly P)) x"
    by (auto simp add: reflect_poly_def Poly_append)
  also have "... = x^(p - degree P) *  x ^ degree P * poly P (inverse x)"
    by (auto simp add: poly_reflect_poly_nz poly_monom hx)
  also have "... = x^p * poly P (inverse x)"
    by (auto simp add: hp power_add[symmetric])
  finally show ?thesis .
qed

lemma reciprocal_fcompose: 
  fixes P::"('a::{ring_char_0,field}) poly" 
  assumes hP: "degree P  p"
  shows "reciprocal_poly p P = monom 1 (p - degree P) * fcompose P 1 [:0, 1:]"
proof (rule poly_eq_by_eval, cases)
  fix x::'a
  assume hx: "x = 0"
  hence "poly (reciprocal_poly p P) x = coeff P p"
    using hP by (auto simp: poly_0_coeff_0 coeff_reciprocal)
  moreover have "poly (monom 1 (p - degree P) 
    * fcompose P 1 [:0, 1:]) x = coeff P p"
  proof cases
    assume "degree P = p"
    thus ?thesis
      apply (induction P arbitrary: p)
      using hx by (auto simp: poly_monom degree_0_iff fcompose_pCons)
  next
    assume "degree P  p"
    hence "degree P < p" using hP by auto
    thus ?thesis
      using hx by (auto simp: poly_monom coeff_eq_0)
  qed
  ultimately show "poly (reciprocal_poly p P) x = poly (monom 1 (p - degree P) * fcompose P 1 [:0, 1:]) x"
    by presburger
next
  fix x::'a assume "x  0"
  thus "poly (reciprocal_poly p P) x = 
        poly (monom 1 (p - degree P) * fcompose P 1 [:0, 1:]) x"
    using hP 
    by (auto simp: poly_reciprocal poly_fcompose inverse_eq_divide
        poly_monom power_diff)
qed

lemma reciprocal_reciprocal: 
  fixes P :: "'a::{field,ring_char_0} poly"
  assumes hP: "degree P  p"
  shows "reciprocal_poly p (reciprocal_poly p P) = P"
proof (rule poly_eq_by_eval)
  fix x
  show "poly (reciprocal_poly p (reciprocal_poly p P)) x = poly P x"
  proof cases
    assume "x = 0"
    thus "poly (reciprocal_poly p (reciprocal_poly p P)) x = poly P x" 
      using hP
      by (auto simp: poly_0_coeff_0 coeff_reciprocal degree_reciprocal)
  next
    assume hx: "x  0"
    hence "poly (reciprocal_poly p (reciprocal_poly p P)) x 
        = x ^ p * (inverse x ^ p * poly P x)" using hP
      by (auto simp: poly_reciprocal degree_reciprocal)
    thus "poly (reciprocal_poly p (reciprocal_poly p P)) x = poly P x" 
      using hP hx left_right_inverse_power right_inverse by auto
  qed
qed

lemma reciprocal_smult: 
  fixes P :: "'a::idom poly" 
  assumes h: "degree P  p"
  shows "reciprocal_poly p (smult n P) = smult n (reciprocal_poly p P)"
proof cases
  assume "n = 0"
  thus ?thesis by (auto simp add: reciprocal_poly_def)
next
  assume "n  0"
  thus ?thesis 
    by (auto simp add: reciprocal_poly_def smult_Poly coeffs_smult 
        rev_map[symmetric])
qed

lemma reciprocal_add: 
  fixes P Q :: "'a::comm_semiring_0 poly"
  assumes "degree P  p" and "degree Q  p" 
  shows "reciprocal_poly p (P + Q) = reciprocal_poly p P + reciprocal_poly p Q" 
(is "?L = ?R")
proof (rule poly_eqI, cases)
  fix n
  assume "n  p"
  then show "coeff ?L n = coeff ?R n"
    using assms by (auto simp: degree_add_le coeff_reciprocal)
next
  fix n assume "¬n  p"
  then show "coeff ?L n = coeff ?R n"
    using assms by (auto simp: degree_add_le coeff_reciprocal_less)
qed 

lemma reciprocal_diff: 
  fixes P Q :: "'a::comm_ring poly"
  assumes "degree P  p" and "degree Q  p" 
  shows "reciprocal_poly p (P - Q) = reciprocal_poly p P - reciprocal_poly p Q"
  by (metis (no_types, lifting) ab_group_add_class.ab_diff_conv_add_uminus assms
      add_diff_cancel degree_add_le degree_minus diff_add_cancel reciprocal_add)

lemma reciprocal_sum: 
  fixes P :: "'a  'b::comm_semiring_0 poly" 
  assumes hP: "k. degree (P k)  p"
  shows "reciprocal_poly p (kA. P k) = (kA. reciprocal_poly p (P k))"
proof (induct A rule: infinite_finite_induct)
  case (infinite A)
  then show ?case by (simp add: reciprocal_0)
next
  case empty
  then show ?case by (simp add: reciprocal_0)
next
  case (insert x F)
  assume "x  F"
     and "reciprocal_poly p (sum P F) = (kF. reciprocal_poly p (P k))"
     and "finite F"
  moreover hence "reciprocal_poly p (sum P (insert x F)) 
      = reciprocal_poly p (P x) + reciprocal_poly p (sum P F)"
    by (auto simp add: reciprocal_add hP degree_sum_le)
  ultimately show "reciprocal_poly p (sum P (insert x F)) 
      = (kinsert x F. reciprocal_poly p (P k))"
    by (auto simp: Groups_Big.comm_monoid_add_class.sum.insert_if)
qed

lemma reciprocal_mult: 
  fixes P Q::"'a::{ring_char_0,field} poly" 
  assumes "degree (P * Q)  p"
    and "degree P  p" and "degree Q  p"
  shows "monom 1 p * reciprocal_poly p (P * Q) = 
         reciprocal_poly p P * reciprocal_poly p Q"
proof (cases "P=0  Q=0")
  case True
  then show ?thesis using assms(1) 
    by (auto simp: reciprocal_fcompose fcompose_mult)
next
  case False
  then show ?thesis  
    using assms
    by (auto simp: degree_mult_eq mult_monom reciprocal_fcompose fcompose_mult)
qed

lemma reciprocal_reflect_poly: 
  fixes P::"'a::{ring_char_0,field} poly" 
  assumes hP: "degree P  p"
  shows "reciprocal_poly p P = monom 1 (p - degree P) * reflect_poly P"
proof (rule poly_eqI)
  fix n
  show "coeff (reciprocal_poly p P) n = 
    coeff (monom 1 (p - degree P) * reflect_poly P) n"
  proof cases
    assume "n  p" 
    thus ?thesis using hP
      by (auto simp: coeff_reciprocal coeff_monom_mult coeff_reflect_poly coeff_eq_0)
  next
    assume "¬ n  p"
    thus ?thesis using hP
      by (auto simp: coeff_reciprocal_less coeff_monom_mult coeff_reflect_poly)
  qed
qed

lemma map_poly_reciprocal: 
  assumes "degree P  p" and "f 0 = 0" 
  shows "map_poly f (reciprocal_poly p P)  = reciprocal_poly p (map_poly f P)"
proof (rule poly_eqI)
  fix n 
  show "coeff (map_poly f (reciprocal_poly p P)) n =
         coeff (reciprocal_poly p (map_poly f P)) n"
  proof (cases "np")
    case True
    then show ?thesis 
      apply (subst coeff_reciprocal[OF True])
      subgoal by (meson assms(1) assms(2) degree_map_poly_le le_trans)
      by (simp add: assms(1) assms(2) coeff_map_poly coeff_reciprocal)
  next
    case False
    then show ?thesis 
      by (metis assms(1) assms(2) coeff_map_poly coeff_reciprocal_less 
          degree_map_poly_le dual_order.trans leI)
  qed
qed

subsection ‹More about @{term proots_count}

lemma proots_count_monom: 
  assumes "0  A" 
  shows "proots_count (monom 1 d) A = 0"
  using assms by (auto simp: proots_count_def poly_monom)

lemma proots_count_reciprocal: 
  fixes P::"'a::{ring_char_0,field} poly"
  assumes hP: "degree P  p" and h0: "P  0" and h0': "0  A"
  shows "proots_count (reciprocal_poly p P) A = proots_count P {x. inverse x  A}"
proof -
  have "proots_count (reciprocal_poly p P) A = 
        proots_count (fcompose P 1 [:0, 1:]) A"
    apply (subst reciprocal_fcompose[OF hP], subst proots_count_times)
    subgoal using h0 by (metis hP reciprocal_0_iff reciprocal_fcompose)
    subgoal using h0' by (auto simp: proots_count_monom)
    done

  also have "... = proots_count P {x. inverse x  A}"
  proof (rule proots_fcompose_bij_eq[symmetric])
    show "bij_betw (λx. poly 1 x / poly [:0, 1:] x) A {x. inverse x  A}"
    proof (rule bij_betw_imageI)
      show "inj_on (λx. poly 1 x / poly [:0, 1:] x) A"
        by (simp add: inj_on_def)

      show "(λx. poly 1 x / poly [:0, 1:] x) ` A = {x. inverse x  A}"
      proof
        show "(λx. poly 1 x / poly [:0, 1:] x) ` A  {x. inverse x  A}"
          by force
        show "{x. inverse x  A}  (λx. poly 1 x / poly [:0, 1:] x) ` A"
        proof
          fix x assume "x  {x::'a. inverse x  A}"
          hence "x = poly 1 (inverse x) / poly [:0, 1:] (inverse x)  inverse x  A"
            by (auto simp: inverse_eq_divide)
          thus "x  (λx. poly 1 x / poly [:0, 1:] x) ` A" by blast
        qed
      qed
    qed
  next
    show "c. 1  smult c [:0, 1:]" 
      by (metis coeff_pCons_0 degree_1 lead_coeff_1 pCons_0_0 pcompose_0'
          pcompose_smult smult_0_right zero_neq_one)
  qed (auto simp: assms infinite_UNIV_char_0)
  finally show ?thesis by linarith
qed

lemma proots_count_reciprocal': 
  fixes P::"real poly"
  assumes hP: "degree P  p" and h0: "P  0"
  shows "proots_count P {x. 0 < x  x < 1} = 
         proots_count (reciprocal_poly p P) {x. 1 < x}"
proof (subst proots_count_reciprocal)
  show "proots_count P {x. 0 < x  x < 1} =
      proots_count P {x. inverse x  Collect ((<) 1)}"
    apply (rule arg_cong[of _ _ "proots_count P"])
    using one_less_inverse_iff by fastforce
qed (use assms in auto)

lemma proots_count_pos: 
  assumes "proots_count P S > 0" 
  shows "x  S. poly P x = 0"
proof (rule ccontr)
  assume "¬ (xS. poly P x = 0)"
  hence "x. x  S  poly P x  0" by blast
  hence "x. x  S  order x P = 0" using order_0I by blast
  hence "proots_count P S = 0" by (force simp: proots_count_def)
  thus False using assms by presburger
qed

lemma proots_count_of_root_set: 
  assumes "P  0" "R  S" and "x. xR  poly P x = 0"
  shows "proots_count P S  card R"
proof -
  have "card R  card (proots_within P S)"
    apply (rule card_mono)
    subgoal using assms by auto
    subgoal using assms(2) assms(3) by (auto simp: proots_within_def)
    done
  also have "...  proots_count P S"
    by (rule card_proots_within_leq[OF assms(1)])
  finally show ?thesis .
qed

lemma proots_count_of_root: assumes "P  0" "xS" "poly P x = 0"
  shows "proots_count P S > 0"
  using proots_count_of_root_set[of P "{x}" S] assms by force

subsection ‹More about @{term changes}

lemma changes_nonneg: "0  changes xs"
  apply (induction xs rule: changes.induct)
  by simp_all

lemma changes_replicate_0: shows "changes (replicate n 0) = 0"
  apply (induction n)
  by auto

lemma changes_append_replicate_0: "changes (xs @ replicate n 0) = changes xs"
proof (induction xs rule: changes.induct)
  case (2 uu)
  then show ?case
    apply (induction n)
    by auto
qed (auto simp: changes_replicate_0)

lemma changes_scale_Cons: 
  fixes xs:: "real list" assumes hs: "s > 0"
  shows "changes (s * x # xs) = changes (x # xs)"
  apply (induction xs rule: changes.induct)
  using assms by (auto simp: mult_less_0_iff zero_less_mult_iff)

lemma changes_scale: 
  fixes xs::"('a::linordered_idom) list"
  assumes hs: "i. i < n  s i > 0" and hn: "length xs  n"
  shows "changes [s i * (nth_default 0 xs i). i  [0..<n]] = changes xs"
using assms
proof (induction xs arbitrary: s n rule: changes.induct)
  case 1
  show ?case by (auto simp: map_replicate_const changes_replicate_0)
next
  case (2 uu)
  show ?case
  proof (cases n)
    case 0
    then show ?thesis by force
  next
    case (Suc m)
    hence "map (λi. s i * nth_default 0 [uu] i) [0..<n] = [s 0 * uu] @ replicate m 0"
    proof (induction m arbitrary: n)
      case (Suc m n)
      from Suc have "map (λi. s i * nth_default 0 [uu] i) [0..<Suc m] = 
                     [s 0 * uu] @ replicate m 0"
        by meson
      hence "map (λi. s i * nth_default 0 [uu] i) [0..<n] =
             [s 0 * uu] @ replicate m 0 @ [0]"
        using Suc by auto
      also have "... = [s 0 * uu] @ replicate (Suc m) 0"
        by (simp add: replicate_append_same)
      finally show ?case .
    qed fastforce
    then show ?thesis
      by (metis changes.simps(2) changes_append_replicate_0)
  qed
next
  case (3 a b xs s n)
  obtain m where hn: "n = m + 2"
    using 3(5)
    by (metis add_2_eq_Suc' diff_diff_cancel diff_le_self length_Suc_conv 
        nat_arith.suc1 ordered_cancel_comm_monoid_diff_class.add_diff_inverse)
  hence h: 
    "map (λi. s i * nth_default 0 (a # b # xs) i) [0..<n] = 
      s 0 * a # s 1 * b # map (λi. 
      (λ i. s (i+2)) i * nth_default 0 (xs) i) [0..<m]"
    apply (induction m arbitrary: n)
    by auto
  consider (neg)"a*b<0" | (nil)"b=0" | (pos)"¬a*b<0  ¬b=0" by linarith
  then show ?case
  proof (cases)
    case neg
    hence 
      "changes (map (λi. s i * nth_default 0 (a # b # xs) i) [0..<n]) =
        1 + changes (s 1 * b # map (λi. (λ i. s (i+2)) i 
                                        * nth_default 0 (xs) i) [0..<m])"
      apply (subst h)
      using 3(4)[of 0] 3(4)[of 1] hn
      by (metis (no_types, lifting) changes.simps(3) mult_less_0_iff pos2
          mult_pos_pos one_less_numeral_iff semiring_norm(76) trans_less_add2)
    also have 
      "... = 1 + changes (map (λi. s (Suc i) * nth_default 0 (b # xs) i) [0..<Suc m])"
      apply (rule arg_cong[of _ _ "λ x. 1 + changes x"])
      apply (induction m)
      by auto
    also have "... = changes (a # b # xs)"
      apply (subst 3(1)[OF neg])
      using 3 neg hn by auto
    finally show ?thesis .
  next
    case nil
    hence "changes (map (λi. s i * nth_default 0 (a # b # xs) i) [0..<n]) =
           changes (s 0 * a # map (λi. (λ i. s (i+2)) i * nth_default 0 (xs) i) [0..<m])"
      apply (subst h)
      using 3(4)[of 0] 3(4)[of 1] hn
      by auto
    also have 
      "... = changes (map 
               (λi. s (if i = 0 then 0 else Suc i) * nth_default 0 (a # xs) i) 
              [0..<Suc m])"
      apply (rule arg_cong[of _ _ "λ x. changes x"])
      apply (induction m)
      by auto
    also have "... = changes (a # b # xs)"
      apply (subst 3(2))
      using 3 nil hn by auto
    finally show ?thesis .
  next
    case pos
    hence "changes (map (λi. s i * nth_default 0 (a # b # xs) i) [0..<n]) =
           changes (s 1 * b # map (λi. (λ i. s (i+2)) i * nth_default 0 (xs) i) [0..<m])"
      apply (subst h)
      using 3(4)[of 0] 3(4)[of 1] hn
      by (metis (no_types, lifting) changes.simps(3) divisors_zero 
          mult_less_0_iff nat_1_add_1 not_square_less_zero one_less_numeral_iff
          semiring_norm(76) trans_less_add2 zero_less_mult_pos zero_less_two)
    also have
      "... = changes (map (λi. s (Suc i) * nth_default 0 (b # xs) i) [0..<Suc m])"
      apply (rule arg_cong[of _ _ "λ x. changes x"])
      apply (induction m)
      by auto
    also have "... = changes (a # b # xs)"
      apply (subst 3(3))
      using 3 pos hn by auto
    finally show ?thesis .
  qed
qed

lemma changes_scale_const: fixes xs::"'a::linordered_idom list" 
  assumes hs: "s  0"
  shows "changes (map ((*) s) xs) = changes xs"
  apply (induction xs rule: changes.induct)
    apply (simp, force)
  using hs by (auto simp: mult_less_0_iff zero_less_mult_iff)

lemma changes_snoc: fixes xs::"'a::linordered_idom list" 
  shows "changes (xs @ [b, a]) = (if a * b < 0 then 1 + changes (xs @ [b])
         else if b = 0 then changes (xs @ [a]) else changes (xs @ [b]))"
  apply (induction xs rule: changes.induct)
  subgoal by (force simp: mult_less_0_iff)
  subgoal by (force simp: mult_less_0_iff)
  subgoal by force
  done

lemma changes_rev: fixes xs:: "'a::linordered_idom list" 
  shows "changes (rev xs) = changes xs"
  apply (induction xs rule: changes.induct)
  by (auto simp: changes_snoc)

lemma changes_rev_about: fixes xs:: "'a::linordered_idom list" 
  shows "changes (replicate (p - length xs) 0 @ rev xs) = changes xs"
proof (induction p)
  case (Suc p)
  then show ?case
  proof cases
    assume "¬Suc p  length xs"
    hence "Suc p - length xs = Suc (p - length xs)" by linarith
    thus ?case using Suc.IH changes_rev by auto
  qed (auto simp: changes_rev)
qed (auto simp: changes_rev)

lemma changes_add_between: 
  assumes "a  x" and "x  b"
  shows "changes (as @ [a, b] @ bs) = changes (as @ [a, x, b] @ bs)"
proof (induction as rule: changes.induct)
  case 1
  then show ?case using assms 
    apply (induction bs)
    by (auto simp: mult_less_0_iff)
next
  case (2 c)
  then show ?case 
    apply (induction bs)
    using assms by (auto simp: mult_less_0_iff)
next
  case (3 y z as)
  then show ?case
    using assms by (auto simp: mult_less_0_iff)
qed

lemma changes_all_nonneg: assumes "i. nth_default 0 xs i  0" shows "changes xs = 0"
  using assms
proof (induction xs rule: changes.induct)
  case (3 x1 x2 xs)
  moreover assume "(i. 0  nth_default 0 (x1 # x2 # xs) i)"
  moreover hence "(i. 0  nth_default 0 (x1 # xs) i)" 
    and "(i. 0  nth_default 0 (x2 # xs) i)"
    and "x1 * x2  0"
  proof -
    fix i
    assume h:"(i. 0  nth_default 0 (x1 # x2 # xs) i)"
    show "0  nth_default 0 (x1 # xs) i"
    proof (cases i)
      case 0
      then show ?thesis using h[of 0] by force
    next
      case (Suc nat)
      then show ?thesis using h[of "Suc (Suc nat)"] by force
    qed
    show "0  nth_default 0 (x2 # xs) i" using h[of "Suc i"] by simp
    show "x1*x2  0" using h[of 0] h[of 1] by simp
  qed
  ultimately show ?case by auto
qed auto

lemma changes_pCons: "changes (coeffs (pCons 0 f)) = changes (coeffs f)"
  by (auto simp: cCons_def)

lemma changes_increasing: 
  assumes "i. i < length xs - 1  xs ! (i + 1)  xs ! i" 
    and "length xs > 1"
    and "hd xs < 0" 
    and "last xs > 0"
  shows "changes xs = 1"
  using assms
proof (induction xs rule:changes.induct)
  case (3 x y xs)
  consider (neg)"x*y<0" | (nil)"y=0" | (pos)"¬x*y<0  ¬y=0" by linarith
  then show ?case
  proof cases
    case neg
    have "changes (y # xs) = 0"
    proof (rule changes_all_nonneg)
      fix i
      show "0  nth_default 0 (y # xs) i"
      proof (cases "i < length (y # xs)")
        case True
        then show ?thesis using 3(4)[of i] 
          apply (induction i)
          subgoal using 3(6) neg by (fastforce simp: mult_less_0_iff)
          subgoal using 3(4) by (auto simp: nth_default_def)
          done
      next
        case False
        then show ?thesis by (simp add: nth_default_def)
      qed
    qed
    thus "changes (x # y # xs) = 1"
      using neg by force
  next
    case nil
    hence "xs  []" using 3(7) by force
    have h: "i. i < length (x # xs) - 1  (x # xs) ! i  (x # xs) ! (i + 1)"
    proof -
      fix i assume "i < length (x # xs) - 1"
      thus "(x # xs) ! i  (x # xs) ! (i + 1)"
        apply (cases "i = 0")
        subgoal using 3(4)[of 0] 3(4)[of 1] xs  [] by force
        using 3(4)[of "i+1"] by simp
    qed
    have "changes (x # xs) = 1"
      apply (rule 3(2))
      using nil h xs  [] 3(6) 3(7) by auto
    thus ?thesis
      using nil by force
  next
    case pos
    hence "xs  []" using 3(6) 3(7) by (fastforce simp: mult_less_0_iff)
    have "changes (y # xs) = 1"
    proof (rule 3(3))
      show "¬ x * y < 0" "y  0"
        using pos by auto
      show "i. i < length (y # xs) - 1 
         (y # xs) ! i  (y # xs) ! (i + 1)"
        using 3(4) by force
      show "1 < length (y # xs)"
        using xs  [] by force
      show "hd (y # xs) < 0"
        using 3(6) pos by (force simp: mult_less_0_iff)
      show "0 < last (y # xs)"
        using 3(7) by force
    qed
    thus ?thesis using pos by auto
  qed
qed auto

end

Theory Bernstein_01

section ‹Bernstein Polynomials over the interval [0, 1]›

theory Bernstein_01
  imports "HOL-Computational_Algebra.Computational_Algebra" 
    "Budan_Fourier.Budan_Fourier"
    "RRI_Misc"
begin

text ‹
The theorem of three circles is a statement about the Bernstein coefficients of a polynomial, the
coefficients when a polynomial is expressed as a sum of Bernstein polynomials. These coefficients
behave nicely under translations and rescaling and are the coefficients of a particular polynomial
in the [0, 1] case. We shall define the [0, 1] case now and consider the general case later,
deriving all the results by rescaling.
›

subsection ‹Definition and basic results›

definition Bernstein_Poly_01 :: "nat  nat  real poly" where
  "Bernstein_Poly_01 j p = (monom (p choose j) j) 
                              * (monom 1 (p-j) p [:1, -1:])"

lemma degree_Bernstein: 
  assumes hb: "j  p" 
  shows "degree (Bernstein_Poly_01 j p) = p"
proof -
  have ha: "monom (p choose j) j  (0::real poly)" using hb by force
  have hb: "monom 1 (p-j) p [:1, -1:]  (0::real poly)"
  proof
    assume "monom 1 (p-j) p [:1, -1:] = (0::real poly)"
    hence "lead_coeff (monom 1 (p - j) p [:1, -1:]) = (0::real)"
      apply (subst leading_coeff_0_iff)
      by simp
    moreover have "lead_coeff (monom (1::real) (p - j) 
        p [:1, -1:]) = (((- 1) ^ (p - j))::real)"
      by (subst lead_coeff_comp, auto simp: degree_monom_eq)
    ultimately show "False" by auto
  qed
  from ha hb show ?thesis
    by (auto simp add: Bernstein_Poly_01_def degree_mult_eq 
          degree_monom_eq degree_pcompose)
qed

lemma coeff_gt: 
  assumes hb: "j > p" 
  shows "Bernstein_Poly_01 j p = 0"
  by (simp add: hb Bernstein_Poly_01_def)

lemma degree_Bernstein_le: "degree (Bernstein_Poly_01 j p)  p"
  apply (cases "j  p")
  by (simp_all add: degree_Bernstein coeff_gt)

lemma poly_Bernstein_nonneg: 
  assumes "x  0" and "1  x" 
  shows "poly (Bernstein_Poly_01 j p) x  0"
  using assms by (simp add: poly_monom poly_pcompose Bernstein_Poly_01_def)

lemma Bernstein_symmetry: 
  assumes "j  p"
  shows "(Bernstein_Poly_01 j p) p [:1, -1:] = Bernstein_Poly_01 (p-j) p"
proof -
  have "(Bernstein_Poly_01 j p) p [:1, -1:]
         = ((monom (p choose j) j) * (monom 1 (p-j) p [:1, -1:])) p [:1, -1:]"
    by (simp add: Bernstein_Poly_01_def)
  also have "... = (monom (p choose (p-j)) j * 
                    (monom 1 (p-j) p [:1, -1:])) p [:1, -1:]" 
    by (fastforce simp: binomial_symmetric[OF assms])
  also have "... = monom (p choose (p-j)) j p [:1, -1:] * 
                   (monom 1 (p-j)) p ([:1, -1:] p [:1, -1:])"
    by (force simp: pcompose_mult pcompose_assoc)
  also have "... = (monom (p choose (p-j)) j p [:1, -1:]) * monom 1 (p-j)"
    by (force simp: pcompose_pCons)
  also have "... = smult (p choose (p-j)) (monom 1 j p [:1, -1:]) 
                    * monom 1 (p-j)" 
    by (simp add: assms smult_monom pcompose_smult[symmetric])
  also have "... = (monom 1 j p [:1, -1:]) * monom (p choose (p-j)) (p-j)"
    apply (subst mult_smult_left)
    apply (subst mult_smult_right[symmetric])
    apply (subst smult_monom)
    by force
  also have "... = Bernstein_Poly_01 (p-j) p" using assms
    by (auto simp: Bernstein_Poly_01_def)
  finally show ?thesis .
qed

subsection @{term Bernstein_Poly_01} and @{term reciprocal_poly}

lemma Bernstein_reciprocal: 
  "reciprocal_poly p (Bernstein_Poly_01 i p) 
    = smult (p choose i) ([:-1, 1:]^(p-i))"
proof cases
  assume "i  p"
  hence "reciprocal_poly p (Bernstein_Poly_01 i p) = 
         reciprocal_poly (degree (Bernstein_Poly_01 i p)) (Bernstein_Poly_01 i p)"
    by (auto simp: degree_Bernstein)
  also have "... = reflect_poly (Bernstein_Poly_01 i p)"
    by (rule reciprocal_degree)
  also have "... = smult (p choose i) ([:-1, 1:]^(p-i))"
    by (auto simp: Bernstein_Poly_01_def reflect_poly_simps monom_altdef
         pcompose_pCons reflect_poly_pCons' hom_distribs)
  finally show ?thesis .
next
  assume h:"¬ i  p"
  hence "reciprocal_poly p (Bernstein_Poly_01 i p) = (0::real poly)"
    by (auto simp: coeff_gt reciprocal_poly_def)
  also have "... = smult (p choose i) ([:-1, 1:]^(p - i))" using h
    by fastforce
  finally show ?thesis .
qed

lemma Bernstein_reciprocal_translate: 
  "reciprocal_poly p (Bernstein_Poly_01 i p) p [:1, 1:] = 
   monom (p choose i) (p - i)"
  by (auto simp: Bernstein_reciprocal pcompose_smult pcompose_pCons monom_altdef hom_distribs)

lemma coeff_Bernstein_sum_01: fixes b::"nat  real" assumes hi: "p  i"
  shows 
    "coeff (reciprocal_poly p 
            (x = 0..p. smult (b x) (Bernstein_Poly_01 x p)) p [:1, 1:]) 
      (p - i) = (p choose i) * (b i)" (is "?L = ?R")
proof -
  define P where "P  (x = 0..p. (smult (b x) (Bernstein_Poly_01 x p)))"

  have "x. degree (smult (b x) (Bernstein_Poly_01 x p))  p"
  proof -
    fix x
    show "degree (smult (b x) (Bernstein_Poly_01 x p))  p"
      apply (cases "x  p")
      by (auto simp: degree_Bernstein coeff_gt)
  qed
  hence "reciprocal_poly p P = 
         (x = 0..p. reciprocal_poly p (smult (b x) (Bernstein_Poly_01 x p)))"
    apply (subst P_def)
    apply (rule reciprocal_sum)
    by presburger
  also have
    "... = (x = 0..p. (smult (b x * (p choose x)) ([:-1, 1:]^(p-x))))"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "reciprocal_poly p (smult (b x) (Bernstein_Poly_01 x p)) =
          smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))"
      by (auto simp add: reciprocal_smult degree_Bernstein Bernstein_reciprocal)
  qed (simp)
  finally have 
    "reciprocal_poly p P = 
     (x = 0..p. (smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))))" .
  hence 
    "(reciprocal_poly p P) p [:1, 1:] = 
     (x = 0..p. (smult ((b x) * (p choose x)) ([:-1, 1:]^(p-x))) p [:1, 1:])"
    by (simp add: pcompose_sum pcompose_add)
  also have "... = (x = 0..p. (monom ((b x) * (p choose x)) (p - x)))"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "smult (b x * (p choose x)) ([:- 1, 1:] ^ (p - x)) p [:1, 1:] =
          monom (b x * (p choose x)) (p - x)"
      by (simp add: hom_distribs pcompose_smult pcompose_pCons monom_altdef)
  qed (simp)
  finally have "(reciprocal_poly p P) p [:1, 1:] = 
                (x = 0..p. (monom ((b x) * (p choose x)) (p - x)))" .
  hence "?L = (x = 0..p. if p - x = p - i then b x * real (p choose x) else 0)"
    by (auto simp add: P_def coeff_sum)
  also have "... = (x = 0..p. if x = i then b x * real (p choose x) else 0)"
  proof (rule sum.cong)
    fix x assume "x  {0..p}"
    hence "x  p" by simp
    thus "(if p - x = p - i then b x * real (p choose x) else 0) =
          (if x = i then b x * real (p choose x) else 0)" using hi
      by (auto simp add: leI)
  qed (simp)
  also have "... = ?R" by simp
  finally show ?thesis .
qed

lemma Bernstein_sum_01: assumes hP: "degree P  p"
  shows 
  "P = (j = 0..p. smult 
     (inverse (real (p choose j)) * 
      coeff (reciprocal_poly p P p [:1, 1:]) (p-j))
   (Bernstein_Poly_01 j p))"
proof -
  define Q where "Q  reciprocal_poly p P p [:1, 1:]"
  from hP Q_def have hQ: "degree Q  p" 
    by (auto simp: degree_reciprocal degree_pcompose)
  have "reciprocal_poly p (j = 0..p. 
        smult (inverse (real (p choose j)) * coeff Q (p-j)) 
        (Bernstein_Poly_01 j p)) p [:1, 1:] = Q"
  proof (rule poly_eqI)
    fix n
    show "coeff (reciprocal_poly p (j = 0..p. 
          smult (inverse (real (p choose j)) * coeff Q (p-j))
          (Bernstein_Poly_01 j p)) p [:1, 1:]) n = coeff Q n" 
      (is "?L = ?R")
    proof cases
      assume hn: "n  p"
      hence "?L = coeff (reciprocal_poly p (j = 0..p. 
             smult (inverse (real (p choose j)) * coeff Q (p-j)) 
             (Bernstein_Poly_01 j p)) p [:1, 1:]) (p - (p - n))"
        by force
      also have "... = (p choose (p-n)) * 
                       (inverse (real (p choose (p-n))) * 
                        coeff Q (p-(p-n)))"
        apply (subst coeff_Bernstein_sum_01)
        by auto
      also have "... = ?R" using hn
        by fastforce
      finally show "?L = ?R" .
    next
      assume hn: "¬ n  p"
      have "degree (j = 0..p.
            smult (inverse (real (p choose j)) * coeff Q (p - j))
            (Bernstein_Poly_01 j p))  p"
      proof (rule degree_sum_le)
        fix q assume "q  {0..p}"
        hence "q  p" by fastforce
        thus "degree (smult (inverse (real (p choose q)) * 
              coeff Q (p - q)) (Bernstein_Poly_01 q p))  p"
          by (auto simp add: degree_Bernstein degree_smult_le)
      qed simp
      hence "degree (reciprocal_poly p (j = 0..p.
            smult (inverse (real (p choose j)) * coeff Q (p - j)) 
            (Bernstein_Poly_01 j p)) p [:1, 1:])  p"
        by (auto simp add: degree_pcompose degree_reciprocal)
      hence "?L = 0" using hn by (auto simp add: coeff_eq_0)
      thus "?L = ?R" using hQ hn by (simp add: coeff_eq_0)
    qed
  qed
  hence "reciprocal_poly p P p [:1, 1:] = 
         reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j))
         (Bernstein_Poly_01 j p)) p [:1, 1:]" 
    by (auto simp: degree_reciprocal degree_pcompose Q_def)
  hence "reciprocal_poly p P p ([:1, 1:] p [:-1, 1:]) =
         reciprocal_poly p (j = 0..p. smult (inverse (real (p choose j)) * 
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) 
         (Bernstein_Poly_01 j p)) p ([:1, 1:] p [:-1, 1:])"
    by (auto simp: pcompose_assoc)
  hence "reciprocal_poly p P = reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p))" 
    by (auto simp: pcompose_pCons)
  hence "reciprocal_poly p (reciprocal_poly p P) = 
         reciprocal_poly p (reciprocal_poly p (j = 0..p. 
         smult (inverse (real (p choose j)) *
         coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p)))"
    by argo
  thus "P = (j = 0..p. smult (inverse (real (p choose j)) * 
        coeff (reciprocal_poly p P p [:1, 1:]) (p-j)) (Bernstein_Poly_01 j p))"
    using hP by (auto simp: reciprocal_reciprocal degree_sum_le degree_smult_le 
                 degree_Bernstein degree_add_le)
qed

lemma Bernstein_Poly_01_span1: 
  assumes hP: "degree P  p"
  shows "P  poly_vs.span {Bernstein_Poly_01 x p | x. x  p}"
proof -
  have "Bernstein_Poly_01 x p
          poly_vs.span {Bernstein_Poly_01 x p |x. x  p}"
    if "x  {0..p}" for x
  proof -
    have "n. Bernstein_Poly_01 x p = Bernstein_Poly_01 n p  n  p"
      using that by force
    then show 
      "Bernstein_Poly_01 x p  poly_vs.span {Bernstein_Poly_01 n p |n. n  p}"
      by (simp add: poly_vs.span_base)
  qed
  thus ?thesis
    apply (subst Bernstein_sum_01[OF hP])
    apply (rule poly_vs.span_sum)
    apply (rule poly_vs.span_scale)
    by blast
qed

lemma Bernstein_Poly_01_span:
  "poly_vs.span {Bernstein_Poly_01 x p | x. x  p} 
      = {x. degree x  p}"
  apply (subst monom_span[symmetric])
  apply (subst poly_vs.span_eq)
  by (auto simp: monom_span degree_Bernstein_le
      Bernstein_Poly_01_span1 degree_monom_eq)

subsection ‹Bernstein coefficients and changes›

definition Bernstein_coeffs_01 :: "nat  real poly  real list" where 
  "Bernstein_coeffs_01 p P = 
   [(inverse (real (p choose j)) * 
    coeff (reciprocal_poly p P p [:1, 1:]) (p-j)). j  [0..<(p+1)]]"

lemma length_Bernstein_coeffs_01: "length (Bernstein_coeffs_01 p P) = p + 1"
  by (auto simp: Bernstein_coeffs_01_def)

lemma nth_default_Bernstein_coeffs_01: assumes "degree P  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p P) i = 
         inverse (p choose i) * coeff (reciprocal_poly p P p [:1, 1:]) (p-i)"
  apply (cases "p = i")
  using assms by (auto simp: Bernstein_coeffs_01_def nth_default_append
                  nth_default_Cons Nitpick.case_nat_unfold binomial_eq_0)

lemma Bernstein_coeffs_01_sum: assumes "degree P  p"
  shows "P = (j = 0..p. smult (nth_default 0 (Bernstein_coeffs_01 p P) j) 
             (Bernstein_Poly_01 j p))"
  apply (subst nth_default_Bernstein_coeffs_01[OF assms])
  apply (subst Bernstein_sum_01[OF assms])
  by argo

definition Bernstein_changes_01 :: "nat  real poly  int" where
  "Bernstein_changes_01 p P = nat (changes (Bernstein_coeffs_01 p P))"

lemma Bernstein_changes_01_def': 
  "Bernstein_changes_01 p P = nat (changes [(inverse (real (p choose j)) * 
     coeff (reciprocal_poly p P p [:1, 1:]) (p-j)). j  [0..<p + 1]])"
  by (simp add: Bernstein_changes_01_def Bernstein_coeffs_01_def)

lemma Bernstein_changes_01_eq_changes: 
  assumes hP: "degree P  p"
  shows "Bernstein_changes_01 p P = 
         changes (coeffs ((reciprocal_poly p P) p [:1, 1:]))"
proof (subst Bernstein_changes_01_def')
  have h: 
    "map (λj. inverse (real (p choose j)) * 
     coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1] = 
     map (λj. inverse (real (p choose j)) * 
     nth_default 0 [nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:]))
                    (p - j). j  [0..<p + 1]] j) [0..<p + 1]"
  proof (rule map_cong)
    fix x
    assume "x  set [0..<p+1]"
    hence hx: "x  p" by fastforce
    moreover have 1:
      "length (map (λj. nth_default 0 
       (coeffs (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1])  Suc p"
      by force
    moreover have "length (coeffs (reciprocal_poly p P p [:1, 1:]))  Suc p"
    proof (cases "P=0")
      case False
      then have "reciprocal_poly p P p [:1, 1:]  0" 
        using hP by (simp add: Missing_Polynomial.pcompose_eq_0 reciprocal_0_iff)
      moreover have "Suc (degree (reciprocal_poly p P p [:1, 1:]))  Suc p"
        using hP by (auto simp: degree_pcompose degree_reciprocal)
      ultimately show ?thesis 
        using length_coeffs_degree by force
    qed (auto simp: reciprocal_0)
    ultimately have h: 
      "nth_default 0 (map (λj. nth_default 0 (coeffs 
       (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) x =
       nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) (p - x)"
      (is "?L = ?R")
    proof -
      have "?L = (map (λj. nth_default 0 (coeffs
            (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) ! x"
        using hx by (auto simp: nth_default_nth)
      also have "... =  nth_default 0 
          (coeffs (reciprocal_poly p P p [:1, 1:])) (p - [0..<p + 1] ! x)"
        apply (subst nth_map)
        using hx by auto
      also have "... = ?R"
        apply (subst nth_upt)
        using hx by auto
      finally show ?thesis .
    qed
    show "inverse (real (p choose x)) *
          coeff (reciprocal_poly p P p [:1, 1:]) (p - x) =
          inverse (real (p choose x)) *
          nth_default 0 (map (λj. nth_default 0 
          (coeffs (reciprocal_poly p P p [:1, 1:])) (p - j)) [0..<p + 1]) x"
      apply (subst h)
      apply (subst nth_default_coeffs_eq)
      by blast
  qed auto

  have 1: 
    "rev (map (λj. nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) 
     (p - j)) [0..<p + 1]) = map (λj. nth_default 0 (coeffs 
     (reciprocal_poly p P p [:1, 1:])) j) [0..<p + 1]"
  proof (subst rev_map, rule map_cong')
    have "q. (q  p  rev [q-p..<q+1] = map ((-) q) [0..<p+1])"
    proof (induction p)
      case 0
      then show ?case by simp
    next
      case (Suc p)
      have IH: "q. (q  p  rev [q-p..<q+1] = map ((-) q) [0..<p+1])"
        using Suc.IH by blast
      show ?case
      proof
        assume hq: "Suc p  q"
        then have h: "rev [q - p..<q + 1] = map ((-) (q)) [0..<p + 1]"
          apply (subst IH)
          using hq by auto
        have "[q - Suc p..<q + 1] = (q - Suc p) # [q - p..<q + 1]"
          by (simp add: Suc_diff_Suc Suc_le_lessD hq upt_conv_Cons)
        hence "rev [q - Suc p..<q + 1] = rev [q - p..<q + 1] @ [q - Suc p]"
          by force
        also have "... = map ((-) (q)) [0..<p + 1] @ [q - Suc p]"
          using h by blast
        also have "... = map ((-) q) [0..<Suc p + 1]"
          by force
        finally show "rev [q - Suc p..<q + 1] = map ((-) q) [0..<Suc p + 1]" .
      qed
    qed
    thus "rev [0..<p + 1] = map ((-) p) [0..<p + 1]"
      by force
  next
    fix y
    assume "y  set [0..<p + 1]"
    hence "y  p" by fastforce
    thus "nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) (p - (p - y)) =
          nth_default 0 (coeffs (reciprocal_poly p P p [:1, 1:])) y"
      by fastforce
  qed

  have 2: " f. f  0  degree f  p 
           map (nth_default 0 (coeffs f)) [0..<p + 1] = 
           coeffs f @ replicate (p - degree f) 0"
  proof (induction p)
    case 0
    then show ?case by (auto simp: degree_0_iff)
  next
    fix f
    case (Suc p)
    hence IH: "(f  0 
                degree f  p 
                map (nth_default 0 (coeffs f)) [0..<p + 1] =
                coeffs f @ replicate (p - degree f) 0)" by blast
    then show ?case
    proof (cases)
      assume h': "Suc p = degree f"
      hence h: "[0..<Suc p + 1] = [0..<length (coeffs f)]" 
        by (metis add_is_0 degree_0 length_coeffs plus_1_eq_Suc zero_neq_one)
      thus ?thesis
        apply (subst h)
        apply (subst map_nth_default)
        using h' by fastforce
    next
      assume h': "Suc p  degree f"
      show ?thesis
      proof
        assume hf: "f  0"
        show "degree f  Suc p 
            map (nth_default 0 (coeffs f)) [0..<Suc p + 1] =
            coeffs f @ replicate (Suc p - degree f) 0"
        proof
          assume "degree f  Suc p"
          hence 1: "degree f  p" using h' by fastforce
          hence 2: "map (nth_default 0 (coeffs f)) [0..<p + 1] =
                  coeffs f @ replicate (p - degree f) 0" using IH hf by blast
          have "map (nth_default 0 (coeffs f)) [0..<Suc p + 1] = 
                map (nth_default 0 (coeffs f)) [0..<p + 1] @
                     [nth_default 0 (coeffs f) (Suc p)]"
            by fastforce
          also have
            "... = coeffs f @ replicate (p - degree f) 0 @ [coeff f (Suc p)]"
            using 2 
            by (auto simp: nth_default_coeffs_eq)
          also have "... = coeffs f @ replicate (p - degree f) 0 @ [0]"
            using ‹degree f  Suc p h' le_antisym le_degree by blast
          also have "... = coeffs f @ replicate (Suc p - degree f) 0" using 1
            by (simp add: Suc_diff_le replicate_app_Cons_same)
          finally show "map (nth_default 0 (coeffs f)) [0..<Suc p + 1] =
                coeffs f @ replicate (Suc p - degree f) 0" .
        qed
      qed
    qed
  qed
  
  thus "int (nat (changes (map (λj. inverse (real (p choose j)) *
        coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1]))) =
        changes (coeffs (reciprocal_poly p P p [:1, 1:]))"
  proof cases
    assume hP: "P = 0"
    show "int (nat (changes (map (λj. inverse (real (p choose j)) *
          coeff (reciprocal_poly p P p [:1, 1:]) (p - j)) [0..<p + 1]))) =
          changes (coeffs (reciprocal_poly p P p [:1, 1:]))" (is "?L = ?R")
    proof -
      have "?L = int (nat (changes (map (λj. 0::real) [0..<p+1])))"
        using hP by (auto simp: reciprocal_0 changes_nonneg)
      also have "... = 0"
        apply (induction p)
        by (auto simp: map_replicate_trivial changes_nonneg
            replicate_app_Cons_same)
      also have "0 = changes ([]::real list)" by simp
      also have "... = ?R" using hP by (auto simp: reciprocal_0)
      finally show ?thesis .
    qed
  next
    assume hP': "P  0"
    thus ?thesis
      apply (subst h)
      apply (subst changes_scale)
        apply auto[2]
      apply (subst changes_rev[symmetric])
      apply (subst 1)
      apply (subst 2)
      apply (simp add: pcompose_eq_0 hP reciprocal_0_iff)
      using assms apply (auto simp: degree_reciprocal)[1]
      by (auto simp: changes_append_replicate_0 changes_nonneg)
  qed
qed

lemma Bernstein_changes_01_test: fixes P::"real poly"
  assumes hP: "degree P  p" and h0: "P  0"
  shows "proots_count P {x. 0 < x  x < 1}  Bernstein_changes_01 p P 
        even (Bernstein_changes_01 p P - proots_count P {x. 0 < x  x < 1})"
proof -
  let ?Q = "(reciprocal_poly p P) p [:1, 1:]"

  have 1: "changes (coeffs ?Q)  proots_count ?Q {x. 0 < x}  
        even (changes (coeffs ?Q) - proots_count ?Q {x. 0 < x})"
    apply (rule descartes_sign)
    by (simp add: Missing_Polynomial.pcompose_eq_0 h0 hP reciprocal_0_iff)
  
  have "((+) (1::real) ` Collect ((<) (0::real))) = {x. (1::real)<x}"
  proof
    show "{x::real. 1 < x}  (+) 1 ` Collect ((<) 0)"
    proof
      fix x::real assume "x  {x. 1 < x}"
      hence "1 < x" by simp
      hence "-1 + x  Collect ((<) 0)" by auto
      hence "1 + (-1 + x)  (+) 1 ` Collect ((<) 0)" by blast
      thus "x  (+) 1 ` Collect ((<) 0)" by argo
    qed
  qed auto
  hence 2:  "proots_count P {x. 0 < x  x < 1} = proots_count ?Q {x. 0 < x}"
    using assms
    by (auto simp: proots_pcompose reciprocal_0_iff proots_count_reciprocal')
  
  show ?thesis
    apply (subst Bernstein_changes_01_eq_changes[OF hP])
    apply (subst Bernstein_changes_01_eq_changes[OF hP])
    apply (subst 2)
    apply (subst 2)
    by (rule 1)
qed

subsection ‹Expression as a Bernstein sum›

lemma Bernstein_coeffs_01_0: "Bernstein_coeffs_01 p 0 = replicate (p+1) 0"
  by (auto simp: Bernstein_coeffs_01_def reciprocal_0 map_replicate_trivial
      replicate_append_same)

lemma Bernstein_coeffs_01_1: "Bernstein_coeffs_01 p 1 = replicate (p+1) 1"
proof -
  have "Bernstein_coeffs_01 p 1 =
     map (λj. inverse (real (p choose j)) *
     coeff (kp. smult (real (p choose k)) ([:0, 1:] ^ k)) (p - j)) [0..<(p+1)]"
    by (auto simp: Bernstein_coeffs_01_def reciprocal_1 monom_altdef
        hom_distribs pcompose_pCons poly_0_coeff_0[symmetric] poly_binomial)
  also have "... = map (λj. inverse (real (p choose j)) * 
             real (p choose (p - j))) [0..<(p+1)]"
    by (auto simp: monom_altdef[symmetric] coeff_sum binomial)
  also have "... = map (λj. 1) [0..<(p+1)]"
    apply (rule map_cong)
    subgoal by argo
    subgoal apply (subst binomial_symmetric)
      by auto
    done
  also have "... = replicate (p+1) 1"
    by (auto simp: map_replicate_trivial replicate_append_same)
  finally show ?thesis .
qed

lemma Bernstein_coeffs_01_x: assumes "p  0"
  shows "Bernstein_coeffs_01 p (monom 1 1) = [i/p. i  [0..<(p+1)]]"
proof -
  have 
    "Bernstein_coeffs_01 p (monom 1 1) = map (λj. inverse (real (p choose j)) *
     coeff (monom 1 (p - Suc 0) p [:1, 1:]) (p - j)) [0..<(p+1)]"
    using assms by (auto simp: Bernstein_coeffs_01_def reciprocal_monom)
  also have 
    "... = map (λj. inverse (real (p choose j)) *
     (kp - Suc 0. coeff (monom (real (p -  1 choose k)) k) (p - j))) [0..<(p+1)]"
    by (auto simp: monom_altdef hom_distribs pcompose_pCons poly_binomial coeff_sum)
  also have"... = map (λj. inverse (real (p choose j)) *
            real (p -  1 choose (p - j))) [0..<(p+1)]"
    by auto
  also have "... = map (λj. j/p) [0..<(p+1)]"
  proof (rule map_cong)
    fix x assume "x  set [0..<(p+1)]"
    hence "x  p" by force
    thus "inverse (real (p choose x)) * real (p - 1 choose (p - x)) =
          real x / real p"
    proof (cases "x = 0")
      show "x = 0  ?thesis"
        using assms by fastforce
      assume 1: "x  p" and 2: "x  0"
      hence "p - x  p - 1" by force
      hence "(p - 1 choose (p - x)) = (p - 1 choose (x - 1))"
        apply (subst binomial_symmetric)
        using 1 2 by auto
      hence "x * (p choose x) = p * (p - 1 choose (p - x))"
         using 2 times_binomial_minus1_eq by simp
       hence "real x * real (p choose x) = real p * real (p - 1 choose (p - x))"
         by (metis of_nat_mult)
       thus ?thesis using 1 2
         by (auto simp: divide_simps)
    qed
  qed blast
  finally show ?thesis .
qed

lemma Bernstein_coeffs_01_add: 
  assumes "degree P  p" and "degree Q  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p (P + Q)) i = 
    nth_default 0 (Bernstein_coeffs_01 p P) i +
    nth_default 0 (Bernstein_coeffs_01 p Q) i"
  using assms by (auto simp: nth_default_Bernstein_coeffs_01 degree_add_le
                    reciprocal_add pcompose_add algebra_simps)

lemma Bernstein_coeffs_01_smult: 
  assumes "degree P  p"
  shows "nth_default 0 (Bernstein_coeffs_01 p (smult a P)) i =
          a * nth_default 0 (Bernstein_coeffs_01 p P) i"
  using assms
  by (auto simp: nth_default_Bernstein_coeffs_01 reciprocal_smult
      pcompose_smult)

end

Theory Bernstein

section ‹Bernstein Polynomials over any finite interval›

theory Bernstein
  imports "Bernstein_01"
begin

subsection ‹Definition and relation to Bernstein Polynomials over [0, 1]›

definition Bernstein_Poly :: "nat  nat  real  real  real poly" where
  "Bernstein_Poly j p c d = smult ((p choose j)/(d - c)^p)
      (((monom 1 j) p [:-c, 1:]) * (monom 1 (p-j) p [:d, -1:]))"

lemma Bernstein_Poly_altdef: 
  assumes "c  d" and "j  p"
  shows "Bernstein_Poly j p c d = smult (p choose j) 
            ([:-c/(d-c), 1/(d-c):]^j * [:d/(d-c), -1/(d-c):]^(p-j))" 
    (is "?L = ?R")
proof -
  have "?L = smult (p choose j) (smult ((1/(d - c))^j)
        (smult ((1/(d - c))^(p-j)) ([:-c, 1:]^j * [:d, -1:]^(p-j))))"
    using assms by (auto simp: Bernstein_Poly_def monom_altdef hom_distribs
                    pcompose_pCons smult_eq_iff field_simps power_add[symmetric])
  also have "... = ?R"
    apply (subst mult_smult_right[symmetric])
    apply (subst mult_smult_left[symmetric])
    apply (subst smult_power)
    apply (subst smult_power)
    by auto
  finally show ?thesis .
qed

lemma Bernstein_Poly_nonneg: 
  assumes "c  x" and "x  d"
  shows "poly (Bernstein_Poly j p c d) x  0"
  using assms by (auto simp: Bernstein_Poly_def poly_pcompose poly_monom)

lemma Bernstein_Poly_01: "Bernstein_Poly j p 0 1 = Bernstein_Poly_01 j p"
  by (auto simp: Bernstein_Poly_def Bernstein_Poly_01_def monom_altdef)

lemma Bernstein_Poly_rescale: 
  assumes "a  b"
  shows "Bernstein_Poly j p c d p [:a, 1:] p [:0, b-a:] 
            = Bernstein_Poly j p ((c-a)/(b-a)) ((d-a)/(b-a))"
  (is "?L = ?R")
proof -
  have "?R = smult (real (p choose j) 
      / ((d - a) / (b - a) - (c - a) / (b - a)) ^ p)
      ([:- ((c - a) / (b - a)), 1:] ^ j 
      * [:(d - a) / (b - a), - 1:] ^ (p - j))"
    by (auto simp: Bernstein_Poly_def monom_altdef hom_distribs 
        pcompose_pCons) 
  also have "... = smult (real (p choose j) / ((d - c) / (b - a)) ^ p)
              ([:- ((c - a) / (b - a)), 1:] ^ j * [:(d - a) / (b - a), - 1:] 
            ^ (p - j))"
    by argo
  also have "... = smult (real (p choose j) / (d - c) ^ p) 
      (smult ((b - a) ^ (p - j)) (smult ((b - a) ^ j)
      ([:- ((c - a) / (b - a)), 1:] ^ j * [:(d - a) / (b - a), - 1:] 
      ^ (p - j))))"
    by (auto simp: power_add[symmetric] power_divide)
  also have "... = smult (real (p choose j) / (d - c) ^ p)
              ([:- (c - a), b - a:] ^ j * [:d - a, -(b - a):] ^ (p - j))"
    apply (subst mult_smult_left[symmetric])
    apply (subst mult_smult_right[symmetric])
    using assms by (auto simp: smult_power)
  also have "... = ?L"
    using assms 
    by (auto simp: Bernstein_Poly_def monom_altdef pcompose_mult 
        pcompose_smult hom_distribs pcompose_pCons)
  finally show ?thesis by presburger
qed

lemma Bernstein_Poly_rescale_01: 
  assumes "c  d"
  shows "Bernstein_Poly j p c d p [:c, 1:] p [:0, d-c:] 
          = Bernstein_Poly_01 j p"
  apply (subst Bernstein_Poly_rescale)
  using assms by (auto simp: Bernstein_Poly_01)

lemma Bernstein_Poly_eq_rescale_01: 
  assumes "c  d"
  shows "Bernstein_Poly j p c d = Bernstein_Poly_01 j p 
            p [:0, 1/(d-c):] p [:-c, 1:]"
  apply (subst Bernstein_Poly_rescale_01[symmetric])
  using assms by (auto simp: pcompose_pCons pcompose_assoc[symmetric])

lemma coeff_Bernstein_sum: 
  fixes b::"nat  real" and p::nat and c d::real
  defines "P  (j = 0..p. (smult (b j) (Bernstein_Poly j p c d)))"
  assumes "i  p" and "c  d"
  shows "coeff ((reciprocal_poly p (P p [:c, 1:] 
      p [:0, d-c:])) p [:1, 1:]) (p - i) = (p choose i) * (b i)"
proof -
  have h: "P p [:c, 1:] p [:0, d-c:] 
      = (j = 0..p. (smult (b j) (Bernstein_Poly_01 j p)))"
    using assms 
    by (auto simp: P_def pcompose_sum pcompose_smult 
          pcompose_add Bernstein_Poly_rescale_01)
  then show ?thesis
    using coeff_Bernstein_sum_01 assms by simp
qed

lemma Bernstein_sum: 
  assumes "c  d" and "degree P  p"
  shows "P = (j = 0..p. smult (inverse (real (p choose j))
     * coeff (reciprocal_poly p (P p [:c, 1:] p [:0, d-c:]) 
      p [:1, 1:]) (p-j)) (Bernstein_Poly j p c d))"
  apply (subst Bernstein_Poly_eq_rescale_01)
  subgoal using assms by blast
  subgoal 
    apply (subst pcompose_smult[symmetric])
    apply (subst pcompose_sum[symmetric])
    apply (subst pcompose_smult[symmetric])
    apply (subst pcompose_sum[symmetric])
    apply (subst Bernstein_sum_01[symmetric])
    using assms by (auto simp: degree_pcompose pcompose_assoc[symmetric] 
        pcompose_pCons)
  done

lemma Bernstein_Poly_span1: 
  assumes "c  d" and "degree P  p"
  shows "P  poly_vs.span {Bernstein_Poly x p c d | x. x  p}"
proof (subst Bernstein_sum[OF assms], rule poly_vs.span_sum)
  fix x :: nat
  assume "x  {0..p}"
  then have "n. Bernstein_Poly x p c d = Bernstein_Poly n p c d  n  p"
    by auto
  then have 
    "Bernstein_Poly x p c d  poly_vs.span {Bernstein_Poly n p c d |n. n  p}"
    by (simp add: poly_vs.span_base)
  thus "smult (inverse (real (p choose x)) *
        coeff (reciprocal_poly p (P p [:c, 1:] p [:0, d - c:]) p [:1, 1:])
        (p - x)) (Bernstein_Poly x p c d)
          poly_vs.span {Bernstein_Poly x p c d |x. x  p}"
    by (rule poly_vs.span_scale)
qed

lemma Bernstein_Poly_span: 
  assumes "c  d" 
  shows "poly_vs.span {Bernstein_Poly x p c d | x. x  p} = {x. degree x  p}"
proof (subst Bernstein_Poly_01_span[symmetric], subst poly_vs.span_eq, rule conjI)
  show "{Bernstein_Poly x p c d |x. x  p}
       poly_vs.span {Bernstein_Poly_01 x p |x. x  p}"
    apply (subst Setcompr_subset)
    apply (rule allI, rule impI)
    apply (rule Bernstein_Poly_01_span1)
    using assms by (auto simp: degree_Bernstein_le Bernstein_Poly_eq_rescale_01
                    degree_pcompose)

  show "{Bernstein_Poly_01 x p |x. x  p}
       poly_vs.span {Bernstein_Poly x p c d |x. x  p}"
    apply (subst Setcompr_subset)
    apply (rule allI, rule impI)
    apply (rule Bernstein_Poly_span1)
    using assms by (auto simp: degree_Bernstein_le)
qed

lemma Bernstein_Poly_independent: assumes "c  d" 
  shows "poly_vs.independent {Bernstein_Poly x p c d | x. x  {..p}}"
proof (rule poly_vs.card_le_dim_spanning)
  show "{Bernstein_Poly x p c d |x. x  {.. p}}  {x. degree x  p}"
    using assms 
    by (auto simp: degree_Bernstein Bernstein_Poly_eq_rescale_01 degree_pcompose)
  show "{x. degree x  p}  poly_vs.span {Bernstein_Poly x p c d |x. x  {..p}}"
    using assms by (auto simp: Bernstein_Poly_span1)
  show "finite {Bernstein_Poly x p c d |x. x  {..p}}" by fastforce
  show "card {Bernstein_Poly x p c d |x. x  {..p}}  poly_vs.dim {x. degree x  p}"
    apply (rule le_trans)
     apply (subst image_Collect[symmetric], rule card_image_le, force)
    by (force simp: dim_degree)
qed

subsection ‹Bernstein coefficients and changes over any interval›

definition Bernstein_coeffs ::
  "nat  real  real  real poly  real list" where 
  "Bernstein_coeffs p c d P = 
    [(inverse (real (p choose j)) * 
      coeff (reciprocal_poly p (P p [:c, 1:] p [:0, d-c:]) p [:1, 1:]) (p-j)). 
     j  [0..<(p+1)]]"

lemma Bernstein_coeffs_eq_rescale: assumes "c  d"
  shows "Bernstein_coeffs p c d P = Bernstein_coeffs_01 p (P p [:c, 1:] p [:0, d-c:])"
  using assms by (auto simp: pcompose_pCons pcompose_assoc[symmetric]
                  Bernstein_coeffs_def Bernstein_coeffs_01_def)

lemma nth_default_Bernstein_coeffs: assumes "degree P  p"
  shows "nth_default 0 (Bernstein_coeffs p c d P) i =
         inverse (p choose i) * coeff
         (reciprocal_poly p (P p [:c, 1:] p [:0, d-c:]) p [:1, 1:]) (p-i)"
  apply (cases "p = i")
  using assms by (auto simp: Bernstein_coeffs_def nth_default_append
                  nth_default_Cons Nitpick.case_nat_unfold binomial_eq_0)

lemma Bernstein_coeffs_sum: assumes "c  d" and hP: "degree P  p"
  shows "P = (j = 0..p. smult (nth_default 0 (Bernstein_coeffs p c d P) j)
         (Bernstein_Poly j p c d))"
  apply (subst nth_default_Bernstein_coeffs[OF hP])
  apply (subst Bernstein_sum[OF assms])
  by argo

definition Bernstein_changes :: "nat  real  real  real poly  int" where
  "Bernstein_changes p c d P = nat (changes (Bernstein_coeffs p c d P))"

lemma Bernstein_changes_eq_rescale: assumes "c  d" and "degree P  p"
  shows "Bernstein_changes p c d P =
         Bernstein_changes_01 p (P p [:c, 1:] p [:0, d-c:])"
  using assms by (auto simp: Bernstein_coeffs_eq_rescale Bernstein_changes_def
                  Bernstein_changes_01_def)

text ‹This is related and mostly equivalent to previous Descartes test \cite{li2019counting}›
lemma Bernstein_changes_test: 
  fixes P::"real poly"
  assumes "degree P  p" and "P  0" and "c < d"
  shows "proots_count P {x. c < x  x < d}  Bernstein_changes p c d P 
        even (Bernstein_changes p c d P - proots_count P {x. c < x  x < d})"
proof -
  define Q where "Q=P p [:c, 1:] p [:0, d - c:]"

  have "int (proots_count Q {x. 0 < x  x < 1}) 
           Bernstein_changes_01 p Q 
              even (Bernstein_changes_01 p Q - 
                  int (proots_count Q {x. 0 < x  x < 1}))"
    unfolding Q_def
    apply (rule Bernstein_changes_01_test)
    subgoal using assms by fastforce
    subgoal using assms by (auto simp: pcompose_eq_0)
    done
  moreover have "proots_count P {x. c < x  x < d} =
            proots_count Q {x. 0 < x  x < 1}"
    unfolding Q_def
  proof (subst proots_pcompose)
    have "poly [:c, 1:] ` poly [:0, d - c:] ` {x. 0 < x  x < 1} =
        {x. c < x  x < d}" (is "?L = ?R")
    proof
      have "c + x * (d - c) < d" if "x < 1" for x
      proof - 
        have "x * (d - c) < 1 * (d - c)"
          using c < d that by force
        then show ?thesis by fastforce
      qed
      then show "?L  ?R"
        using assms by auto
    next
      show "?R  ?L"
      proof
        fix x::real assume "x  ?R"
        hence "c < x" and "x < d" by auto
        thus "x  ?L"
        proof (subst image_eqI)
          show "x = poly [:c, 1:] (x - c)" by force
          assume "c < x" and "x < d"
          thus "x - c  poly [:0, d - c:] ` {x. 0 < x  x < 1}"
          proof (subst image_eqI)
            show "x - c = poly [:0, d - c:] ((x - c)/(d - c))"
              using assms by fastforce
            assume "c < x" and "x < d"
            thus "(x - c) / (d - c)  {x. 0 < x  x < 1}"
              by auto
          qed fast
        qed fast
      qed
    qed
    then show "proots_count P {x. c < x  x < d} =
        proots_count (P p [:c, 1:]) 
        (poly [:0, d - c:] ` {x. 0 < x  x < 1})"
      using assms by (auto simp:proots_pcompose)
    show "P p [:c, 1:]  0"
      by (simp add: pcompose_eq_0 assms(2))
    show "degree [:0, d - c:] = 1"
      using assms by auto
  qed
  moreover have " Bernstein_changes p c d P = Bernstein_changes_01 p Q"
    unfolding Q_def
    apply (rule Bernstein_changes_eq_rescale)
    using assms by auto
  ultimately show ?thesis by auto
qed

subsection ‹The control polygon of a polynomial›

definition control_points ::
  "nat  real  real  real poly  (real × real) list"
where
  "control_points p c d P = 
   [(((real i)*d + (real (p - i))*c)/p, 
      nth_default 0 (Bernstein_coeffs p c d P) i).
      i  [0..<(p+1)]]"

lemma line_above: 
  fixes a b c d :: real and p :: nat and P :: "real poly"
  assumes hline: "i. i  p  a * (((real i)*d + (real (p - i))*c)/p) + b 
                  nth_default 0 (Bernstein_coeffs p c d P) i"
  and hp: "p  0" and hcd: "c  d" and hP: "degree P  p"
  shows "x. c  x  x  d  a*x + b  poly P x"
proof -
  fix x
  assume hc: "c  x" and  hd: "x  d"

  have bern_eq:"Bernstein_coeffs p c d [:b, a:] =
           [a*(real i * d + real (p - i) * c)/p + b. i  [0..<(p+1)]]"
  proof -
    have "Bernstein_coeffs p c d [:b, a:] = map (nth_default 0
          (Bernstein_coeffs_01 p ([:b, a:] p [:c, 1:] p [:0, d - c:])))
         [0..<p+1]"
      apply (subst Bernstein_coeffs_eq_rescale["OF" hcd])
      apply (subst map_nth_default[symmetric])
      apply (subst length_Bernstein_coeffs_01)
      by blast
    also have 
      "... = map (λi. a * (real i * d + real (p - i) * c) / real p + b) [0..<p + 1]"
    proof (rule map_cong)
      fix x assume hx: "x  set [0..<p + 1]"
      have "nth_default 0 (Bernstein_coeffs_01 p
            ([:b, a:] p [:c, 1:] p [:0, d - c:])) x =
            nth_default 0 (Bernstein_coeffs_01 p
            (smult (b + a*c) 1 + smult (a*(d - c)) (monom 1 1))) x"
      proof-
        have "[:b, a:] p [:c, 1:] p [:0, d - c:] =
                  smult (b + a*c) 1 + smult (a*(d - c)) (monom 1 1)"
          by (simp add: monom_altdef pcompose_pCons)
        then show ?thesis by auto
      qed
      also have "... = 
          nth_default 0 (Bernstein_coeffs_01 p (smult (b + a * c) 1)) x +
          nth_default 0 (Bernstein_coeffs_01 p (smult (a * (d - c)) (monom 1 1))) x"
        apply (subst Bernstein_coeffs_01_add)
        using hp by (auto simp: degree_monom_eq)
      also have "...  =
            (b + a*c) * nth_default 0 (Bernstein_coeffs_01 p 1) x +
            (a*(d - c)) * nth_default 0 (Bernstein_coeffs_01 p (monom 1 1)) x"
        apply (subst Bernstein_coeffs_01_smult)
        using hp by (auto simp: Bernstein_coeffs_01_smult degree_monom_eq)
      also have "... =
          (b + a * c) * (if x < p + 1 then 1 else 0) +
           a * (d - c) * (real (nth_default 0 [0..<p + 1] x) / real p)" 
        apply (subst Bernstein_coeffs_01_1, subst Bernstein_coeffs_01_x[OF hp])
        apply (subst nth_default_replicate_eq, subst nth_default_map_eq[of _ 0])
        by auto
      also have "... =
              (b + a * c) * (if x < p + 1 then 1 else 0) +
              a * (d - c) * (real ([0..<p + 1] ! x) / real p)"
        apply (subst nth_default_nth)
        using hx by auto
      also have "... = (b + a * c) * (if x < p + 1 then 1 else 0) +
              a * (d - c) * (real (0 + x) / real p)"
        apply (subst nth_upt)
        using hx by auto
      also have "... = a * (real x * d + real (p - x) * c) / real p + b"
        apply (subst of_nat_diff)
        using hx hp by (auto simp: field_simps)
      finally show "nth_default 0 (Bernstein_coeffs_01 p
                    ([:b, a:] p [:c, 1:] p [:0, d - c:])) x =
                    a * (real x * d + real (p - x) * c) / real p + b" .
    qed blast
    finally show ?thesis .
  qed

  have nth_default_geq:"nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i 
           nth_default 0 (Bernstein_coeffs p c d P) i" for i
  proof -
    show "nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i 
          nth_default 0 (Bernstein_coeffs p c d P) i"
    proof cases
      define p1 where "p1  p+1"
      assume h: "i  p"
      hence "nth_default 0 (Bernstein_coeffs p c d P) i 
             a * (((real i)*d + (real (p - i))*c)/p) + b"
        by (rule hline)
      also have "... = nth_default 0 (map (λi. a * (real i * d 
          + real (p - i) * c) / real p + b) [0..<p + 1]) i"
        apply (subst p1_def[symmetric])
        using h apply (auto simp: nth_default_def)
        by (auto simp: p1_def)
      also have "... = nth_default 0 (Bernstein_coeffs p c d [:b, a:]) i"
        using bern_eq by simp
      finally show ?thesis .
    next
      assume h: "¬i  p"
      thus ?thesis
        using assms 
        by (auto simp: nth_default_def Bernstein_coeffs_eq_rescale
                        length_Bernstein_coeffs_01)
    qed
  qed
  
  have "poly P x = (k = 0..p.
        poly (smult (nth_default 0 (Bernstein_coeffs p c d P) k)
        (Bernstein_Poly k p c d)) x)"
    apply (subst Bernstein_coeffs_sum[OF hcd hP])
    by (rule poly_sum)
  also have "...  (k = 0..p.
      poly (smult (nth_default 0 (Bernstein_coeffs p c d [:b, a:]) k)
        (Bernstein_Poly k p c d)) x)"
    apply (rule sum_mono)
    using mult_right_mono[OF nth_default_geq] Bernstein_Poly_nonneg[OF hc hd]
    by auto
  also have "... = poly [:b, a:] x"
    apply (subst(2) Bernstein_coeffs_sum[of c d "[:b, a:]" p])
    using assms apply auto[2]
    by (rule poly_sum[symmetric])
  also have "... = a*x + b" by force
  finally show "poly P x  a*x + b" .
qed

end

Theory Normal_Poly

section ‹Normal Polynomials›

theory Normal_Poly
  imports "RRI_Misc"
begin

text ‹
Here we define normal polynomials as defined in
  Basu, S., Pollack, R., Roy, M.-F.: Algorithms in Real Algebraic Geometry. 
  Springer Berlin Heidelberg, Berlin, Heidelberg (2016).
›

definition normal_poly :: "('a::{comm_ring_1,ord}) poly  bool" where
"normal_poly p 
  (p  0) 
  ( i. 0  coeff p i) 
  ( i. coeff p i * coeff p (i+2)  (coeff p (i+1))^2) 
  ( i j k. i  j  j  k  0 < coeff p i 
       0 < coeff p k  0 < coeff p j)"

lemma normal_non_zero: "normal_poly p  p  0" 
  using normal_poly_def by blast

lemma normal_coeff_nonneg: "normal_poly p  0  coeff p i" 
  using normal_poly_def by metis

lemma normal_poly_coeff_mult: 
    "normal_poly p  coeff p i * coeff p (i+2)  (coeff p (i+1))^2" 
  using normal_poly_def by blast

lemma normal_poly_pos_interval: 
    "normal_poly p  i  j  j  k  0 < coeff p i  0 < coeff p k 
       0 < coeff p j" 
  using normal_poly_def by blast

lemma normal_polyI:
  assumes "(p  0)"
      and "( i. 0  coeff p i)"
      and "( i. coeff p i * coeff p (i+2)  (coeff p (i+1))^2)"
      and "( i j k. i  j  j  k  0 < coeff p i  0 < coeff p k  0 < coeff p j)"
    shows "normal_poly p"
  using assms by (force simp: normal_poly_def)

lemma linear_normal_iff: 
  fixes x::real 
  shows "normal_poly [:-x, 1:]  x  0"
proof
  assume "normal_poly [:-x, 1:]"
  thus "x  0" using normal_coeff_nonneg[of "[:-x, 1:]" 0] by auto
next
  assume "x  0"
  then have "0  coeff [:- x, 1:] i" for i
    by (cases i) (simp_all add: pCons_one)
  moreover have "0 < coeff [:- x, 1:] j"
    if "i  j" "j  k" "0 < coeff [:- x, 1:] i" 
        "0 < coeff [:- x, 1:] k" for i j k
    apply (cases "k=0  i=0")
    subgoal using that 
      by (smt (z3) bot_nat_0.extremum_uniqueI degree_pCons_eq_if 
          le_antisym le_degree not_less_eq_eq)
    subgoal using that 
      by (smt (z3) One_nat_def degree_pCons_eq_if le_degree less_one
          not_le one_neq_zero pCons_one verit_la_disequality)
    done
  ultimately show "normal_poly [:-x, 1:]"
    unfolding normal_poly_def by auto
qed

lemma quadratic_normal_iff: 
  fixes z::complex 
  shows "normal_poly [:(cmod z)2, -2*Re z, 1:] 
           Re z  0  4*(Re z)^2  (cmod z)^2"
proof
  assume "normal_poly [:(cmod z)2, - 2 * Re z, 1:]"
  hence "-2*Re z  0  (cmod z)^2  0  (-2*Re z)^2  (cmod z)^2"
    using normal_coeff_nonneg[of _ 1] normal_poly_coeff_mult[of _ 0] 
    by fastforce
  thus "Re z  0  4*(Re z)^2  (cmod z)^2"
    by auto
next
  assume asm:"Re z  0  4*(Re z)^2  (cmod z)^2"
  define P where "P=[:(cmod z)2, - 2 * Re z, 1:]"

  have "0  coeff P i" for i 
    unfolding P_def using asm
    apply (cases "i=0i=1i=2")
    by (auto simp:numeral_2_eq_2 coeff_eq_0)
  moreover have "coeff P i * coeff P (i + 2)  (coeff P (i + 1))2" for i
    apply (cases "i=0i=1i=2")
    using asm 
    unfolding P_def by (auto simp:coeff_eq_0)
  moreover have "0 < coeff P j"
    if "0 < coeff P k" "0 < coeff P i" "j  k" "i  j"
    for i j k
    using that unfolding P_def 
    apply (cases "k=0  k=1  k=2")
    subgoal using asm
      by (smt (z3) One_nat_def Suc_1 bot_nat_0.extremum_uniqueI 
          coeff_pCons_0 coeff_pCons_Suc le_Suc_eq 
          zero_less_power2)
    subgoal by (auto simp:coeff_eq_0)
    done
  moreover have "P0" unfolding P_def by auto
  ultimately show "normal_poly P"
    unfolding normal_poly_def by blast
qed

lemma normal_of_no_zero_root: 
  fixes f::"real poly" 
  assumes hzero: "poly f 0  0" and hdeg: "i  degree f" 
    and hnorm: "normal_poly f"
  shows "0 < coeff f i"
proof -
  have "coeff f 0 > 0" using hzero normal_coeff_nonneg[OF hnorm]
    by (metis eq_iff not_le_imp_less poly_0_coeff_0)
  moreover have "coeff f (degree f) > 0" using normal_coeff_nonneg[OF hnorm] normal_non_zero[OF hnorm]
    by (meson dual_order.irrefl eq_iff eq_zero_or_degree_less not_le_imp_less)
  moreover have "0  i" by simp
  ultimately show "0 < coeff f i" using hdeg normal_poly_pos_interval[OF hnorm] by blast
qed

lemma normal_divide_x: 
  fixes f::"real poly" 
  assumes hnorm: "normal_poly (f*[:0,1:])"
  shows "normal_poly f"
proof (rule normal_polyI)
  show "f  0"
    using normal_non_zero[OF hnorm] by auto
next
  fix i
  show "0  coeff f i"
    using normal_coeff_nonneg[OF hnorm, of "Suc i"] by (simp add: coeff_pCons)
next
  fix i
  show "coeff f i * coeff f (i + 2)  (coeff f (i + 1))2"
    using normal_poly_coeff_mult[OF hnorm, of "Suc i"] by (simp add: coeff_pCons)
next
  fix i j k
  show "i  j  j  k  0 < coeff f i  0 < coeff f k  0 < coeff f j"
    using normal_poly_pos_interval[of _ "Suc i" "Suc j" "Suc k", OF hnorm]
    by (simp add: coeff_pCons)
qed

lemma normal_mult_x: 
  fixes f::"real poly" 
  assumes hnorm: "normal_poly f"
  shows "normal_poly (f * [:0, 1:])"
proof (rule normal_polyI)
  show "f * [:0, 1:]  0"
    using normal_non_zero[OF hnorm] by auto
next
  fix i
  show "0  coeff (f * [:0, 1:]) i"
    using normal_coeff_nonneg[OF hnorm, of "i-1"] by (cases i, auto simp: coeff_pCons)
next
  fix i
  show "coeff (f * [:0, 1:]) i * coeff (f * [:0, 1:]) (i + 2)  (coeff (f * [:0, 1:]) (i + 1))2"
    using normal_poly_coeff_mult[OF hnorm, of "i-1"] by (cases i, auto simp: coeff_pCons)
next
  fix i j k
  show "i  j  j  k  0 < coeff (f * [:0, 1:]) i  0 < coeff (f * [:0, 1:]) k  0 < coeff (f * [:0, 1:]) j"
    using normal_poly_pos_interval[of _ "i-1" "j-1" "k-1", OF hnorm]
    apply (cases i, force)
    apply (cases j, force)
    apply (cases k, force)
    by (auto simp: coeff_pCons)
qed

lemma normal_poly_general_coeff_mult: 
  fixes f::"real poly" 
  assumes "normal_poly f" and "h  j"
  shows "coeff f (h+1) * coeff f (j+1)  coeff f h * coeff f (j+2)"
using assms proof (induction j)
  case 0
  then show ?case
    using normal_poly_coeff_mult by (auto simp: power2_eq_square)[1]
next
  case (Suc j)
  then show ?case
  proof (cases "h = Suc j")
    assume "h = Suc j" "normal_poly f"
    thus ?thesis
      using normal_poly_coeff_mult by (auto simp: power2_eq_square)
  next
    assume "(normal_poly f 
       h  j  coeff f h * coeff f (j + 2)  coeff f (h + 1) * coeff f (j + 1))"
      "normal_poly f" and h: "h  Suc j" "h  Suc j"
    hence IH: "coeff f h * coeff f (j + 2)  coeff f (h + 1) * coeff f (j + 1)"
      by linarith
    show ?thesis
    proof (cases "coeff f (Suc j + 1) = 0", cases "coeff f (Suc j + 2) = 0")
      show "coeff f (Suc j + 1) = 0  coeff f (Suc j + 2) = 0 
        coeff f h * coeff f (Suc j + 2)  coeff f (h + 1) * coeff f (Suc j + 1)"
        by (metis assms(1) mult_zero_right normal_coeff_nonneg)
    next
      assume 1: "coeff f (Suc j + 1) = 0" "coeff f (Suc j + 2)  0"
      hence "coeff f (Suc j + 2) > 0" "¬coeff f (Suc j + 1) > 0" 
        using normal_coeff_nonneg[of f "Suc j + 2"] assms(1) by auto
      hence "coeff f h > 0  False"
        using normal_poly_pos_interval[of f h "Suc j + 1" "Suc j + 2"] assms(1) h by force
      hence "coeff f h = 0"
        using normal_coeff_nonneg[OF assms(1)] less_eq_real_def by auto
      thus "coeff f h * coeff f (Suc j + 2)  coeff f (h + 1) * coeff f (Suc j + 1)"
        using 1 by fastforce
    next
      assume 1: "coeff f (Suc j + 1)  0"
      show "coeff f h * coeff f (Suc j + 2)  coeff f (h + 1) * coeff f (Suc j + 1)"
      proof (cases "coeff f (Suc j) = 0")
        assume 2: "coeff f (Suc j) = 0"
        hence "coeff f (Suc j + 1) > 0" "¬coeff f (Suc j) > 0" 
          using normal_coeff_nonneg[of f "Suc j + 1"] assms(1) 1 by auto
        hence "coeff f h > 0  False"
          using normal_poly_pos_interval[of f h "Suc j" "Suc j + 1"] assms(1) h by force
        hence "coeff f h = 0"
          using normal_coeff_nonneg[OF assms(1)] less_eq_real_def by auto
        thus "coeff f h * coeff f (Suc j + 2)  coeff f (h + 1) * coeff f (Suc j + 1)"
          by (simp add: assms(1) normal_coeff_nonneg)
      next
        assume 2: "coeff f (Suc j)  0"
        from normal_poly_coeff_mult[OF assms(1), of "Suc j"] normal_coeff_nonneg[OF assms(1), of "Suc j"]
          normal_coeff_nonneg[OF assms(1), of "Suc (Suc j)"] 1 2 
        have 3: "coeff f (Suc j + 1) / coeff f (Suc j)  coeff f (Suc j + 2) / coeff f (Suc j + 1)"
          by (auto simp: power2_eq_square divide_simps algebra_simps)
        have "(coeff f h * coeff f (j + 2)) * (coeff f (Suc j + 2) / coeff f (Suc j + 1))  (coeff f (h + 1) * coeff f (j + 1)) * (coeff f (Suc j + 1) / coeff f (Suc j))"
          apply (rule mult_mono[OF IH])
          using 3 by (simp_all add: assms(1) normal_coeff_nonneg)
        thus "coeff f h * coeff f (Suc j + 2)  coeff f (h + 1) * coeff f (Suc j + 1)"
          using 1 2 by fastforce
      qed
    qed
  qed
qed

lemma normal_mult: 
  fixes f g::"real poly"
  assumes hf: "normal_poly f" and hg: "normal_poly g"
  defines "df  degree f" and "dg  degree g"
  shows "normal_poly (f*g)"
using df_def hf proof (induction df arbitrary: f)
text ‹We shall first show that without loss of generality we may assume poly f 0 ≠ 0›,
      this is done by induction on the degree, if 0 is a root then we derive the result from f/[:0,1:]›.›
  fix f::"real poly" fix i::nat
  assume "0 = degree f" and hf: "normal_poly f"
  then obtain a where "f = [:a:]" using degree_0_iff by auto
  then show "normal_poly (f*g)"
    apply (subst normal_polyI)
    subgoal using normal_non_zero[OF hf] normal_non_zero[OF hg] by auto
    subgoal 
      using normal_coeff_nonneg[of _ 0, OF hf] normal_coeff_nonneg[OF hg] 
      by simp
    subgoal 
      using normal_coeff_nonneg[of _ 0, OF hf] normal_poly_coeff_mult[OF hg] 
      by (auto simp: algebra_simps power2_eq_square mult_left_mono)[1]
    subgoal 
      using normal_non_zero[OF hf] normal_coeff_nonneg[of _ 0, OF hf] normal_poly_pos_interval[OF hg]
      by (simp add: zero_less_mult_iff)
    subgoal by simp
    done
next
  case (Suc df)
  then show ?case
  proof (cases "poly f 0 = 0")
    assume "poly f 0 = 0" and hf:"normal_poly f"
    moreover then obtain f' where hdiv: "f = f'*[:0,1:]"
      by (smt (verit) dvdE mult.commute poly_eq_0_iff_dvd) 
    ultimately have hf': "normal_poly f'" using normal_divide_x by blast
    assume "Suc df = degree f"
    hence "degree f' = df" using hdiv normal_non_zero[OF hf'] by (auto simp: degree_mult_eq)
    moreover assume "f. df = degree f  normal_poly f  normal_poly (f * g)"
    ultimately have "normal_poly (f'*g)" using hf' by blast
    thus "normal_poly (f*g)" using hdiv normal_mult_x by fastforce
  next
    assume hf: "normal_poly f" and hf0: "poly f 0  0"
    define dg where "dg  degree g"
    show "normal_poly (f * g)"
    using dg_def hg proof (induction dg arbitrary: g)
      text ‹Similarly we may assume poly g 0 ≠ 0›.›
      fix g::"real poly" fix i::nat
      assume "0 = degree g" and hg: "normal_poly g"
      then obtain a where "g = [:a:]" using degree_0_iff by auto
      then show "normal_poly (f*g)"
        apply (subst normal_polyI)
        subgoal 
          using normal_non_zero[OF hg] normal_non_zero[OF hf] by auto
        subgoal 
          using normal_coeff_nonneg[of _ 0, OF hg] normal_coeff_nonneg[OF hf] 
          by simp
        subgoal 
          using normal_coeff_nonneg[of _ 0, OF hg] normal_poly_coeff_mult[OF hf] 
          by (auto simp: algebra_simps power2_eq_square mult_left_mono)
        subgoal 
          using normal_non_zero[OF hf] normal_coeff_nonneg[of _ 0, OF hg] 
            normal_poly_pos_interval[OF hf]
          by (simp add: zero_less_mult_iff)
        by simp
    next
      case (Suc dg)
      then show ?case
      proof (cases "poly g 0 = 0")
        assume "poly g 0 = 0" and hg:"normal_poly g"
        moreover then obtain g' where hdiv: "g = g'*[:0,1:]"
          by (smt (verit) dvdE mult.commute poly_eq_0_iff_dvd) 
        ultimately have hg': "normal_poly g'" using normal_divide_x by blast
        assume "Suc dg = degree g"
        hence "degree g' = dg" using hdiv normal_non_zero[OF hg'] by (auto simp: degree_mult_eq)
        moreover assume "g. dg = degree g  normal_poly g  normal_poly (f * g)"
        ultimately have "normal_poly (f*g')" using hg' by blast
        thus "normal_poly (f*g)" using hdiv normal_mult_x by fastforce
      next
        text ‹It now remains to show that $(fg)_i \geq 0$. This follows by decomposing $\{(h, j) \in
              \mathbf{Z}^2 | h > j\} = \{(h, j) \in \mathbf{Z}^2 | h \leq j\} \cup \{(h, h - 1) \in 
              \mathbf{Z}^2 | h \in \mathbf{Z}\}$.
              Note in order to avoid working with infinite sums over integers all these sets are
              bounded, which adds some complexity compared to the proof of lemma 2.55 in
              Basu, S., Pollack, R., Roy, M.-F.: Algorithms in Real Algebraic Geometry.
              Springer Berlin Heidelberg, Berlin, Heidelberg (2016).›
        assume hg0: "poly g 0  0" and hg: "normal_poly g"
        have "f * g  0" using hf hg by (simp add: normal_non_zero Suc.prems)
        moreover have "i. coeff (f*g) i  0"
          apply (subst coeff_mult, rule sum_nonneg, rule mult_nonneg_nonneg)
          using normal_coeff_nonneg[OF hf] normal_coeff_nonneg[OF hg] by auto
        moreover have "
          coeff (f*g) i * coeff (f*g) (i+2)  (coeff (f*g) (i+1))^2" for i
        proof -

          text ‹$(fg)_{i+1}^2 - (fg)_i(fg)_{i+2} = \left(\sum_x f_xg_{i+1-x}\right)^2 - 
                \left(\sum_x f_xg_{i+2-x}\right)\left(\sum_x f_xg_{i-x}\right)$›
          have "(coeff (f*g) (i+1))^2 - coeff (f*g) i * coeff (f*g) (i+2) = 
              (xi+1. coeff f x * coeff g (i + 1 - x)) *
              (xi+1. coeff f x * coeff g (i + 1 - x)) -
              (xi+2. coeff f x * coeff g (i + 2 - x)) *
              (xi. coeff f x * coeff g (i - x))"
            by (auto simp: coeff_mult power2_eq_square algebra_simps)
          
          text ‹$\dots = \sum_{x, y} f_xg_{i+1-x}f_yg_{i+1-y} - \sum_{x, y} f_xg_{i+2-x}f_yg_{i-y}$›
          also have "... =
              (xi+1. yi+1. coeff f x * coeff g (i + 1 - x) *
                                  coeff f y * coeff g (i + 1 - y)) -
              (xi+2. yi. coeff f x * coeff g (i + 2 - x) *
                                coeff f y * coeff g (i - y))"
            by (subst sum_product, subst sum_product, auto simp: algebra_simps)

          text ‹$\dots = \sum_{h \leq j} f_hg_{i+1-h}f_jg_{i+1-j} + \sum_{h>j} f_hg_{i+1-h}f_jg_{i+1-j}
                       - \sum_{h \leq j} f_hg_{i+2-h}f_jg_{i-j} - \sum_{h>j} f_hg_{i+2-h}f_jg_{i-j}$›
          also have "... =
              ((h, j){(h, j). i+1  j  j  h}. 
                coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j)) +
              ((h, j){(h, j). i+1  h  h > j}. 
                coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j)) -
             (((h, j){(h, j). i  j  j  h}. 
                coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j)) +
              ((h, j){(h, j). i + 2  h  h > j  i  j}.
                coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j)))"
          proof -
            have "(xi + 1. yi + 1. coeff f x * coeff g (i + 1 - x) * coeff f y * coeff g (i + 1 - y)) =
              ((h, j){(h, j). j  i + 1  h  j}.
                 coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j)) +
              ((h, j){(h, j). h  i + 1  j < h}.
                 coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
            proof (subst sum.union_disjoint[symmetric])
              have H:"{(h, j). j  i + 1  h  j}  {..i+1} × {..i+1}"
                     "{(h, j). h  i + 1  j < h}  {..i+1} × {..i+1}"
                     "finite ({..i+1} × {..i+1})"
                by (fastforce, fastforce, fastforce)
              show "finite {(h, j). j  i + 1  h  j}"
                apply (rule finite_subset) using H by (blast, blast)
              show "finite {(h, j). h  i + 1  j < h}"
                apply (rule finite_subset) using H by (blast, blast)
              show "{(h, j). j  i + 1  h  j}  {(h, j). h  i + 1  j < h} = {}"
                by fastforce
              show "(xi + 1. yi + 1. coeff f x * coeff g (i + 1 - x) * coeff f y * coeff g (i + 1 - y)) =
                  ((h, j){(h, j). j  i + 1  h  j}  {(h, j). h  i + 1  j < h}.
                    coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
                apply (subst sum.cartesian_product, rule sum.cong)
                apply force by blast
            qed
            moreover have "(xi + 2. yi. coeff f x * coeff g (i + 2 - x) * coeff f y * coeff g (i - y)) =
               ((h, j){(h, j). j  i  h  j}.
                  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j)) +
               ((h, j){(h, j). i + 2  h  h > j  i  j}.
                  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
            proof (subst sum.union_disjoint[symmetric])
              have H:"{(h, j). j  i  h  j}  {..i+2} × {..i}"
                     "{(h, j). i + 2  h  h > j  i  j}  {..i+2} × {..i}"
                     "finite ({..i+2} × {..i})"
                by (fastforce, fastforce, fastforce)
              show "finite {(h, j). j  i  h  j}"
                apply (rule finite_subset) using H by (blast, blast)
              show "finite {(h, j). i + 2  h  h > j  i  j}"
                apply (rule finite_subset) using H by (blast, blast)
              show "{(h, j). j  i  h  j}  {(h, j). i + 2  h  h > j  i  j} = {}"
                by fastforce
              show "(xi + 2. yi. coeff f x * coeff g (i + 2 - x) * coeff f y * coeff g (i - y)) =
                  ((h, j){(h, j). j  i  h  j}  {(h, j). i + 2  h  h > j  i  j}.
                     coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
                apply (subst sum.cartesian_product, rule sum.cong)
                apply force by blast
            qed
            ultimately show ?thesis by presburger
          qed

          text ‹$\dots = \sum_{h \leq j} f_hg_{i+1-h}f_jg_{i+1-j} + \sum_{h \leq j} f_{j+1}g_{i-j}f_{h-2}g_{i+2-h} 
                       + \sum_h f_hg_{i+1-h}f_{h-1}g_{i+2-h} - \sum_{h \leq j} f_hg_{i+2-h}f_jg_{i-j}
                       - \sum_{h \leq j} f_{j+1}g_{i+1-j}f_{h-2}g_{i+1-h} - \sum_h f_hg_{i+2-h}f_{h-1}g_{i+1-h}$›
          also have "... =
              ((h, j){(h, j). j  i + 1  h  j}.
                 coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j)) +
              ((h, j){(h, j). j  i  h  j  0 < h}.
                 coeff f (j+1) * coeff g (i - j) * coeff f (h-1) * coeff g (i + 2 - h)) +
              (h{1..i+1}.
                 coeff f h * coeff g (i + 1 - h) * coeff f (h-1) * coeff g (i + 2 - h)) -
              (((h, j){(h, j). j  i  h  j}.
                  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j)) +
               ((h, j){(h, j). j  i + 1  h  j  0 < h}.
                  coeff f (j+1) * coeff g (i + 1 - j) * coeff f (h-1) * coeff g (i + 1 - h)) +
               (h{1..i+1}.
                  coeff f h * coeff g (i + 2 - h) * coeff f (h-1) * coeff g (i + 1 - h)))"
          proof -
            have "((h, j){(h, j). h  i + 1  j < h}.
                   coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j)) = 
                  ((h, j){(h, j). j  i  h  j  0 < h}.
                   coeff f (j + 1) * coeff g (i - j) * coeff f (h - 1) * coeff g (i + 2 - h)) +
                  (h = 1..i + 1. coeff f h * coeff g (i + 1 - h) * coeff f (h - 1) * coeff g (i + 2 - h))"
            proof -
              have 1: "((h, j){(h, j). j  i  h  j  0 < h}.
                     coeff f (j + 1) * coeff g (i - j) * coeff f (h - 1) * coeff g (i + 2 - h)) =
                    ((h, j){(h, j). h  i + 1  j < h  h  j + 1}.
                     coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
              proof (rule sum.reindex_cong)
                show "{(h, j). j  i  h  j  0 < h} = (λ(h, j). (j+1, h-1)) ` {(h, j). h  i + 1  j < h  h  j + 1}"
                proof
                  show "(λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 1  j < h  h  j + 1}  {(h, j). j  i  h  j  0 < h}"
                    by fastforce
                  show "{(h, j). j  i  h  j  0 < h}  (λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 1  j < h  h  j + 1}"
                  proof
                    fix x
                    assume "x  {(h, j). j  i  h  j  0 < h}"
                    then obtain h j where "x = (h, j)" "j  i" "h  j" "0 < h" by blast
                    hence "j + 1  i + 1  h - 1 < j + 1  j + 1  h - 1 + 1  x = ((h - 1) + 1, (j + 1) - 1)"
                      by auto
                    thus "x  (λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 1  j < h  h  j + 1}"
                      by (auto simp: image_iff)
                  qed
                qed
                show "inj_on (λ(h, j). (j + 1, h - 1)) {(h, j). h  i + 1  j < h  h  j + 1}"
                proof
                  fix x y::"nat×nat"
                  assume "x  {(h, j). h  i + 1  j < h  h  j + 1}" "y  {(h, j). h  i + 1  j < h  h  j + 1}"
                  thus "(case x of (h, j)  (j + 1, h - 1)) = (case y of (h, j)  (j + 1, h - 1))  x = y"
                    by auto
                qed
                show "x. x  {(h, j). h  i + 1  j < h  h  j + 1} 
                 (case case x of (h, j)  (j + 1, h - 1) of
                  (h, j)  coeff f (j + 1) * coeff g (i - j) * coeff f (h - 1) * coeff g (i + 2 - h)) =
                 (case x of (h, j)  coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
                  by fastforce
              qed
              have 2: "(h = 1..i + 1. coeff f h * coeff g (i + 1 - h) * coeff f (h - 1) * coeff g (i + 2 - h)) =
                    ((h, j){(h, j). h  i + 1  j < h  h = j + 1}.
                     coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
              proof (rule sum.reindex_cong)
                show "{1..i + 1} = fst ` {(h, j). h  i + 1  j < h  h = j + 1}"
                proof
                  show "{1..i + 1}  fst ` {(h, j). h  i + 1  j < h  h = j + 1}"
                  proof
                    fix x
                    assume "x  {1..i + 1}"
                    hence "x  i + 1  x - 1 < x  x = x - 1 + 1  x = fst (x, x-1)"
                      by auto
                    thus "x  fst ` {(h, j). h  i + 1  j < h  h = j + 1}"
                      by blast
                  qed
                  show "fst ` {(h, j). h  i + 1  j < h  h = j + 1}  {1..i + 1}"
                    by force
                qed
                show "inj_on fst {(h, j). h  i + 1  j < h  h = j + 1}"
                proof
                  fix x y
                  assume "x  {(h, j). h  i + 1  j < h  h = j + 1}"
                         "y  {(h, j). h  i + 1  j < h  h = j + 1}"
                  hence "x = (fst x, fst x - 1)" "y = (fst y, fst y - 1)" "fst x > 0" "fst y > 0"
                    by auto
                  thus "fst x = fst y  x = y" by presburger
                qed
                show "x. x  {(h, j). h  i + 1  j < h  h = j + 1} 
                   coeff f (fst x) * coeff g (i + 1 - fst x) * coeff f (fst x - 1) * coeff g (i + 2 - fst x) =
                   (case x of (h, j)  coeff f h * coeff g (i + 1 - h) * coeff f j * coeff g (i + 1 - j))"
                  by fastforce
              qed
              have H: "{(h, j). h  i + 1  j < h  h  j + 1}  {0..i+1}×{0..i+1}"
                      "{(h, j). h  i + 1  j < h  h = j + 1}  {0..i+1}×{0..i+1}"
                      "finite ({0..i+1}×{0..i+1})"
                by (fastforce, fastforce, fastforce)
              have "finite {(h, j). h  i + 1  j < h  h  j + 1}"
                   "finite {(h, j). h  i + 1  j < h  h = j + 1}"
                apply (rule finite_subset) using H apply (simp, simp)
                apply (rule finite_subset) using H apply (simp, simp)
                done
              thus ?thesis 
                apply (subst 1, subst 2, subst sum.union_disjoint[symmetric])
                   apply auto[3]
                apply (rule sum.cong)
                by auto
            qed
            moreover have "((h, j){(h, j). h  i + 2  j < h  j  i}.
                  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j)) = 
               ((h, j){(h, j). j  i + 1  h  j  0 < h}.
                  coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) +
               (h = 1..i + 1. coeff f h * coeff g (i + 2 - h) * coeff f (h - 1) * coeff g (i + 1 - h))"
            proof -
              have 1: "((h, j){(h, j). j  i + 1  h  j  0 < h}.
                     coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) =
                    ((h, j){(h, j). h  i + 2  j < h  j  i  h  j + 1}.
                     coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
              proof (rule sum.reindex_cong)
                show "{(h, j). j  i + 1  h  j  0 < h} = (λ(h, j). (j+1, h-1)) ` {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                proof
                  show "(λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 2  j < h  j  i  h  j + 1}  {(h, j). j  i + 1  h  j  0 < h}"
                    by fastforce
                  show "{(h, j). j  i + 1  h  j  0 < h}  (λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                  proof
                    fix x
                    assume "x  {(h, j). j  i + 1  h  j  0 < h}"
                    then obtain h j where "x = (h, j)" "j  i + 1" "h  j" "0 < h" by blast
                    hence "j + 1  i + 2  h - 1 < j + 1  h - 1  i  j + 1  h - 1 + 1  x = ((h - 1) + 1, (j + 1) - 1)"
                      by auto
                    thus "x  (λ(h, j). (j + 1, h - 1)) ` {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                      by (auto simp: image_iff)
                  qed
                qed
                show "inj_on (λ(h, j). (j + 1, h - 1)) {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                proof
                  fix x y::"nat×nat"
                  assume "x  {(h, j). h  i + 2  j < h  j  i  h  j + 1}" "y  {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                  thus "(case x of (h, j)  (j + 1, h - 1)) = (case y of (h, j)  (j + 1, h - 1))  x = y"
                    by auto
                qed
                show "x. x  {(h, j). h  i + 2  j < h  j  i  h  j + 1} 
                   (case case x of (h, j)  (j + 1, h - 1) of
                    (h, j)  coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) =
                   (case x of (h, j)  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
                  by fastforce
              qed
              have 2: "(h = 1..i + 1. coeff f h * coeff g (i + 2 - h) * coeff f (h - 1) * coeff g (i + 1 - h)) =
                    ((h, j){(h, j). h  i + 2  j < h  j  i  h = j + 1}.
                     coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
              proof (rule sum.reindex_cong)
                show "{1..i + 1} = fst ` {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                proof
                  show "{1..i + 1}  fst ` {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                  proof
                    fix x
                    assume "x  {1..i + 1}"
                    hence "x  i + 2  x - 1 < x  x - 1  i  x = x - 1 + 1  x = fst (x, x-1)"
                      by auto
                    thus "x  fst ` {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                      by blast
                  qed
                  show "fst ` {(h, j). h  i + 2  j < h  j  i  h = j + 1}  {1..i + 1}"
                    by force
                qed
                show "inj_on fst {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                proof
                  fix x y
                  assume "x  {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                         "y  {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                  hence "x = (fst x, fst x - 1)" "y = (fst y, fst y - 1)" "fst x > 0" "fst y > 0"
                    by auto
                  thus "fst x = fst y  x = y" by presburger
                qed
                show "x. x  {(h, j). h  i + 2  j < h  j  i  h = j + 1} 
                   coeff f (fst x) * coeff g (i + 2 - fst x) * coeff f (fst x - 1) * coeff g (i + 1 - fst x) =
                   (case x of (h, j)  coeff f h * coeff g (i + 2 - h) * coeff f j * coeff g (i - j))"
                  by fastforce
              qed
              have H: "{(h, j). h  i + 2  j < h  j  i  h  j + 1}  {0..i+2}×{0..i}"
                      "{(h, j). h  i + 2  j < h  j  i  h = j + 1}  {0..i+2}×{0..i}"
                      "finite ({0..i+2}×{0..i})"
                by (fastforce, fastforce, fastforce)
              have "finite {(h, j). h  i + 2  j < h  j  i  h  j + 1}"
                   "finite {(h, j). h  i + 2  j < h  j  i  h = j + 1}"
                apply (rule finite_subset) using H apply (simp, simp)
                apply (rule finite_subset) using H apply (simp, simp)
                done
              thus ?thesis 
                apply (subst 1, subst 2, subst sum.union_disjoint[symmetric])
                   apply auto[3]
                apply (rule sum.cong)
                by auto
            qed
            ultimately show ?thesis
              by algebra
          qed

          text ‹$\dots = \sum_{h \leq j} f_hf_j\left(g_{i+1-h}g_{i+1-j} - g_{i+2-h}g_{i-j}\right) + 
                         \sum_{h \leq j} f_{j+1}f_{h-1}\left(g_{i-j}g_{i+2-h} - g_{i+1-j}f_jg_{i+1-h}\right)$

                Note we have to also consider the edge cases caused by making these sums finite.›
          also have "... =
              ((h, j){(h, j). j = i + 1  h  j}.
                 coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
              ((h, j){(h, j). j  i  h  j}.
                 coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j))) +
              ((h, j){(h, j). j  i  h  j  0 < h}.
                 coeff f (j+1) * coeff f (h-1) * (coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h))) -
              (((h, j){(h, j). j = i + 1  h  j  0 < h}.
                  coeff f (j+1) * coeff g (i + 1 - j) * coeff f (h-1) * coeff g (i + 1 - h)))" (is "?L = ?R")
          proof -
            have "?R = 
              ((h, j){(h, j). j = i + 1  h  j}.
                 coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
              (((h, j){(h, j). j  i  h  j}.
                 coeff f h * coeff f j * coeff g (i + 1 - h) * coeff g (i + 1 - j)) -
              ((h, j){(h, j). j  i  h  j}.
                 coeff f h * coeff f j * coeff g (i + 2 - h) * coeff g (i - j))) +
              (((h, j){(h, j). j  i  h  j  0 < h}.
                 coeff f (j+1) * coeff f (h-1) * coeff g (i - j) * coeff g (i + 2 - h)) - 
              ((h, j){(h, j). j  i  h  j  0 < h}.
                 coeff f (j+1) * coeff f (h-1) * coeff g (i + 1 - j) * coeff g (i + 1 - h))) -
              (((h, j){(h, j). j = i + 1  h  j  0 < h}.
                  coeff f (j+1) * coeff g (i + 1 - j) * coeff f (h-1) * coeff g (i + 1 - h)))"
              apply (subst sum_subtractf[symmetric], subst sum_subtractf[symmetric])
              by (auto simp: algebra_simps split_beta)
            also have "... =
                (((h, j){(h, j). j = i + 1  h  j}.
                   coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
                ((h, j){(h, j). j  i  h  j}.
                   coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j)))) -
                 ((h, j){(h, j). j  i  h  j}.
                    coeff f h * coeff f j * coeff g (i + 2 - h) * coeff g (i - j)) +
                ((h, j){(h, j). j  i  h  j  0 < h}.
                    coeff f (j + 1) * coeff f (h - 1) * coeff g (i - j) * coeff g (i + 2 - h)) -
                (((h, j){(h, j). j  i  h  j  0 < h}.
                   coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) +
                ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                   coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)))"
              by (auto simp: algebra_simps)
            also have "... = ?L"
            proof -
              have "((h, j){(h, j). j = i + 1  h  j}.
                       coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
                    ((h, j){(h, j). j  i  h  j}.
                       coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) =
                    ((h, j){(h, j). j  i + 1  h  j}.
                       coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j)))"
              proof (subst sum.union_disjoint[symmetric])
                have "{(h, j). j = i + 1  h  j}  {..i + 1} × {..i + 1}"
                     "{(h, j). j  i  h  j}  {..i + 1} × {..i + 1}"
                  by (fastforce, fastforce)
                thus "finite {(h, j). j = i + 1  h  j}" "finite {(h, j). j  i  h  j}"
                  by (auto simp: finite_subset)
                show "{(h, j). j = i + 1  h  j}  {(h, j). j  i  h  j} = {}"
                  by fastforce
              qed (rule sum.cong, auto)
              moreover have "((h, j){(h, j). j  i  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) +
                 ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) =
                 ((h, j){(h, j). j  i + 1  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h))"
              proof (subst sum.union_disjoint[symmetric])
                have "{(h, j). j  i  h  j  0 < h}  {..i + 1} × {..i + 1}"
                     "{(h, j). j = i + 1  h  j  0 < h}  {..i + 1} × {..i + 1}"
                  by (fastforce, fastforce)
                thus "finite {(h, j). j  i  h  j  0 < h}" "finite {(h, j). j = i + 1  h  j  0 < h}"
                  by (auto simp: finite_subset)
                show "{(h, j). j  i  h  j  0 < h}  {(h, j). j = i + 1  h  j  0 < h} = {}"
                  by fastforce
              qed (rule sum.cong, auto)
              ultimately show ?thesis 
                by (auto simp: algebra_simps)
            qed
            finally show ?thesis by presburger
          qed

          text ‹$\dots = \sum_{h \leq j} \left(f_hf_j - f_{j+1}f_{h-1}\right)
                                         \left(g_{i+1-h}g_{i+1-j} - g_{i+2-h}g_{i-j}\right)$›
          also have "... = 
              ((h, j){(h, j). j  i  h  j  0 < h}.
                 -(coeff f h * coeff f j - coeff f (j+1) * coeff f (h-1)) * (coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h))) +
              ((h, j){(h, j). j  i  h  j  h = 0}.
                 coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j))) +
              ((h, j){(h, j). j = i + 1  h  j}.
                 coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
              (((h, j){(h, j). j = i + 1  h  j  0 < h}.
                  coeff f (j+1) * coeff g (i + 1 - j) * coeff f (h-1) * coeff g (i + 1 - h)))" (is "?L = ?R")
          proof -
            have "((h, j){(h, j). j  i  h  j}.
               coeff f h * coeff f j *
               (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j))) =
             ((h, j){(h, j). j  i  h  j  0 < h}.
               coeff f h * coeff f j *
               (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j))) +
             ((h, j){(h, j). j  i  h  j  0 = h}.
               coeff f h * coeff f j *
               (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j)))"
            proof (subst sum.union_disjoint[symmetric])
              have "{(h, j). j  i  h  j  0 < h}  {..i}×{..i}" "{(h, j). j  i  h  j  0 = h}  {..i}×{..i}"
                by (force, force)
              thus "finite {(h, j). j  i  h  j  0 < h}" "finite {(h, j). j  i  h  j  0 = h}"
                by (auto simp: finite_subset)
              show "{(h, j). j  i  h  j  0 < h}  {(h, j). j  i  h  j  0 = h} = {}"
                by fast
            qed (rule sum.cong, auto)
            
            moreover have "((h, j){(h, j). j  i  h  j  0 < h}.
               (-coeff f h * coeff f j + coeff f (j + 1) * coeff f (h - 1)) *
               (coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h))) =
                ((h, j){(h, j). j  i  h  j  0 < h}.
                   coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j))) +
                ((h, j){(h, j). j  i  h  j  0 < h}.
                   coeff f (j + 1) * coeff f (h - 1) *
                   (coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h)))"
              by (subst sum.distrib[symmetric], rule sum.cong, fast, auto simp: algebra_simps)

            ultimately show ?thesis
              by (auto simp: algebra_simps)
          qed

          text ‹$\dots \geq 0$ by normal_poly_general_coeff_mult›
          also have "...  0"
          proof -
            have "0  ((h, j){(h, j). j  i  h  j  0 < h}.
                        - (coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)) *
                        (coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h)))"
            proof (rule sum_nonneg)
              fix x assume "x  {(h, j). j  i  h  j  0 < h}"
              then obtain h j where H: "x = (h, j)" "j  i" "h  j" "0 < h" by fast
              hence "h - 1  j - 1" by force
              hence 1: "coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)  0"
                using normal_poly_general_coeff_mult[OF hf, of "h-1" "j-1"] H
                by (auto simp: algebra_simps)
              from H have "i - j  i - h" by force
              hence 2: "coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h)  0"
                using normal_poly_general_coeff_mult[OF hg, of "i - j" "i - h"] H
                by (smt (verit, del_insts) Nat.add_diff_assoc2 le_trans)
              show "0  (case x of
               (h, j) 
                 - (coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)) *
                 (coeff g (i - j) * coeff g (i + 2 - h) -
                  coeff g (i + 1 - j) * coeff g (i + 1 - h)))"
                apply (subst H(1), subst split, rule mult_nonpos_nonpos, subst neg_le_0_iff_le)
                subgoal using 1 by blast
                subgoal using 2 by blast
                done
            qed
            moreover have "0  ((h, j){(h, j). j  i  h  j  h = 0}.
                        coeff f h * coeff f j *
                        (coeff g (i + 1 - h) * coeff g (i + 1 - j) - coeff g (i + 2 - h) * coeff g (i - j)))"
            proof (rule sum_nonneg)
              fix x assume "x  {(h, j). j  i  h  j  h = 0}"
              then obtain h j where H: "x = (h, j)" "j  i" "h  j" "h = 0" by fast
              have 1: "coeff f h * coeff f j  0"
                by (simp add: hf normal_coeff_nonneg)
              from H have "i - j  i - h" by force
              hence 2: "coeff g (i - j) * coeff g (i + 2 - h) - coeff g (i + 1 - j) * coeff g (i + 1 - h)  0"
                using normal_poly_general_coeff_mult[OF hg, of "i - j" "i - h"] H
                by (smt (verit, del_insts) Nat.add_diff_assoc2 le_trans)
              show "0  (case x of
               (h, j) 
                 coeff f h * coeff f j *
                 (coeff g (i + 1 - h) * coeff g (i + 1 - j) -
                  coeff g (i + 2 - h) * coeff g (i - j)))"
                apply (subst H(1), subst split, rule mult_nonneg_nonneg)
                subgoal using 1 by blast
                subgoal using 2 by argo
                done
            qed
            moreover have "0  ((h, j){(h, j). j = i + 1  h  j}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
                 ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h))"
            proof -
              have "((h, j){(h, j). j = i + 1  h  j}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
                 ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) = 
                 ((h, j){(h, j). j = i + 1  h  j  h = 0}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
                 ((h, j){(h, j). j = i + 1  h  j  0 < h}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
                 ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                    coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h))"
              proof (subst sum.union_disjoint[symmetric])
                have "{(h, j). j = i + 1  h  j  h = 0} = {(0, i + 1)}"
                     "{(h, j). j = i + 1  h  j  0 < h} = {1..i+1} × {i + 1}"
                  by (fastforce, force)
                thus "finite {(h, j). j = i + 1  h  j  h = 0}"
                     "finite {(h, j). j = i + 1  h  j  0 < h}"
                  by auto
                show "{(h, j). j = i + 1  h  j  h = 0}  {(h, j). j = i + 1  h  j  0 < h} = {}" 
                  by fastforce
                have "{(h, j). j = i + 1  h  j  h = 0}  {(h, j). j = i + 1  h  j  0 < h} = {(h, j). j = i + 1  h  j}"
                  by fastforce
                thus "((h, j){(h, j). j = i + 1  h  j}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
                      ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                         coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h)) =
                      ((h, j){(h, j). j = i + 1  h  j  h = 0}  {(h, j). j = i + 1  h  j  0 < h}.
                         coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) -
                      ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                         coeff f (j + 1) * coeff g (i + 1 - j) * coeff f (h - 1) * coeff g (i + 1 - h))"
                  by presburger
              qed
              also have "... =
                ((h, j){(h, j). j = i + 1  h  j  h = 0}. coeff f h * coeff f j * (coeff g (i + 1 - h) * coeff g (i + 1 - j))) +
                ((h, j){(h, j). j = i + 1  h  j  0 < h}. 
                   (coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)) * (coeff g (i + 1 - h) * coeff g (i + 1 - j)))"
                by (subst add_diff_eq[symmetric], subst sum_subtractf[symmetric], subst add_left_cancel, rule sum.cong, auto simp: algebra_simps)
              also have "...  0"
              proof (rule add_nonneg_nonneg)
                show "0  ((h, j){(h, j). j = i + 1  h  j  0 < h}.
                        (coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)) *
                        (coeff g (i + 1 - h) * coeff g (i + 1 - j)))"
                proof (rule sum_nonneg)
                  fix x assume "x  {(h, j). j = i + 1  h  j  0 < h}"
                  then obtain h j where H: "x = (h, j)" "j = i + 1" "h  j" "0 < h" by fast
                  hence "h - 1  j - 1" by force
                  hence 1: "coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)  0"
                    using normal_poly_general_coeff_mult[OF hf, of "h-1" "j-1"] H
                    by (auto simp: algebra_simps)
                  hence 2: "0  coeff g (i + 1 - h) * coeff g (i + 1 - j)"
                    by (meson hg mult_nonneg_nonneg normal_coeff_nonneg)
                  show "0  (case x of
                   (h, j) 
                     (coeff f h * coeff f j - coeff f (j + 1) * coeff f (h - 1)) *
                     (coeff g (i + 1 - h) * coeff g (i + 1 - j)))"
                    apply (subst H(1), subst split, rule mult_nonneg_nonneg)
                    subgoal using 1 by blast
                    subgoal using 2 by blast
                    done
                qed
              qed (rule sum_nonneg, auto simp: hf hg normal_coeff_nonneg)[1]
              finally show ?thesis .
            qed
            ultimately show ?thesis by auto
          qed
          finally show "coeff (f * g) i * coeff (f * g) (i + 2)  (coeff (f * g) (i + 1))2" by (auto simp: power2_eq_square)
        qed
        moreover have "i j k. i  j  j  k  0 < coeff (f*g) i  0 < coeff (f*g) k  0 < coeff (f*g) j"
        proof -
          fix j k
          assume "0 < coeff (f * g) k"
          hence "k  degree (f * g)" using le_degree by force
          moreover assume "j  k"
          ultimately have "j  degree (f * g)" by auto
          hence 1: "j  degree f + degree g"
            by (simp add: degree_mult_eq hf hg normal_non_zero)
          show "0 < coeff (f * g) j"
            apply (subst coeff_mult, rule sum_pos2[of _ "min j (degree f)"], simp, simp)
            apply (rule mult_pos_pos, rule normal_of_no_zero_root, simp add: hf0, simp)
            using hf apply auto[1]
             apply (rule normal_of_no_zero_root)
               apply (simp add: hg0)
            using 1 apply force
            using hg apply auto[1]
            by (simp add: hf hg normal_coeff_nonneg)
        qed
        ultimately show "normal_poly (f*g)" 
          by (rule normal_polyI)
      qed
    qed
  qed
qed

lemma normal_poly_of_roots: 
  fixes p::"real poly"
  assumes "z. poly (map_poly complex_of_real p) z = 0 
         Re z  0  4*(Re z)^2  (cmod z)^2"
      and "lead_coeff p = 1"
  shows "normal_poly p"
  using assms
proof (induction p rule: real_poly_roots_induct)
  fix p::"real poly" and x::real
  assume "lead_coeff (p * [:- x, 1:]) = 1"
  hence 1: "lead_coeff p = 1"
    by (metis coeff_degree_mult lead_coeff_pCons(1) mult_cancel_left1 pCons_one zero_neq_one)
  assume h: "(z. poly (map_poly complex_of_real (p * [:- x, 1:])) z = 0 
                 Re z  0  (cmod z)2  4 * (Re z)2)"
  hence 2: "(z. poly (map_poly complex_of_real p) z = 0 
                 Re z  0  (cmod z)2  4 * (Re z)2)"
    by (metis (no_types, hide_lams) four_x_squared mult.commute mult_cancel_left1 of_real_poly_map_mult poly_mult)
  have 3: "normal_poly [:-x, 1:]"
    apply (subst linear_normal_iff, 
        subst Re_complex_of_real[symmetric], rule conjunct1)
    by (rule h[of x], subst of_real_poly_map_poly[symmetric], force)
  assume "(z. poly (map_poly complex_of_real p) z = 0
              Re z  0  (cmod z)2  4 * (Re z)2) 
            lead_coeff p = 1  normal_poly p"
  hence "normal_poly p" using 1 2 by fast
  then show "normal_poly (p * [: