# 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 = (∑k≤n. monom (of_nat (n choose k)) k)"
proof -
have "[:(1::'a::comm_ring_1), 1:]^n = (monom 1 1 + 1)^n"
monom_altdef pCons_one power_one_right smult_1_left)
also have "... = (∑k≤n. of_nat (n choose k) * monom 1 1 ^ k)"
apply (subst binomial_ring)
by force
also have "... = (∑k≤n. 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

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 = (∑i≤p. monom (coeff x i) i)"
by (simp add: poly_as_sum_of_monoms')
also have
"... = (∑i≤p. 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': "(∑i≤p. 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)

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))"
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 ‹n≤p› 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)"
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

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"

lemma reciprocal_sum:
fixes P :: "'a ⇒ 'b::comm_semiring_0 poly"
assumes hP: "⋀k. degree (P k) ≤ p"
shows "reciprocal_poly p (∑k∈A. P k) = (∑k∈A. 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) = (∑k∈F. 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))
= (∑k∈insert 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 "n≤p")
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 "¬ (∃x∈S. 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. x∈R ⟹ 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" "x∈S" "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 @ "
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
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)

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)"
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:])"
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
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 @ "
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
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)
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 (∑k≤p. 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)) *
(∑k≤p - 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

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

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
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"
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
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

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=0∨i=1∨i=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=0∨i=1∨i=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 "P≠0" 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)
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)
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) =
(∑x≤i+1. coeff f x * coeff g (i + 1 - x)) *
(∑x≤i+1. coeff f x * coeff g (i + 1 - x)) -
(∑x≤i+2. coeff f x * coeff g (i + 2 - x)) *
(∑x≤i. 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 "... =
(∑x≤i+1. ∑y≤i+1. coeff f x * coeff g (i + 1 - x) *
coeff f y * coeff g (i + 1 - y)) -
(∑x≤i+2. ∑y≤i. 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 "(∑x≤i + 1. ∑y≤i + 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 "(∑x≤i + 1. ∑y≤i + 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 "(∑x≤i + 2. ∑y≤i. 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 "(∑x≤i + 2. ∑y≤i. 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
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
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"
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)
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
apply (rule normal_of_no_zero_root)
apply (simp add: hg0)
using 1 apply force
using hg apply auto
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 * [: