Theory Three_Circles
section ‹Proof of the theorem of three circles›
theory Three_Circles
imports "Bernstein" "Normal_Poly"
begin
text ‹
The theorem of three circles is a result in real algebraic geometry about the number of real roots
in an interval. It says if the number of roots in certain circles in the complex plane are zero or
one then the number of roots in the circles is equal to the sign changes of the Bernstein
coefficients on that interval for which the circles intersect the real line. This can then be used
to determine if an interval has a real root in the bisection procedure, which is more efficient than
Descartes' rule of signs.
The proof here follows Theorem 10.50 in
Basu, S., Pollack, R., Roy, M.-F.: Algorithms in Real Algebraic Geometry.
Springer Berlin Heidelberg, Berlin, Heidelberg (2016).
This theorem has also been fomalised in Coq \<^cite>‹"zsido2014theorem"›. The relationship
between this theorem and root isolation has been elaborated in Eigenwillig's PhD
thesis \<^cite>‹"eigenwillig2008real"›.
›
subsection ‹No sign changes case›
declare degree_pcompose[simp del]
corollary descartes_sign_zero:
fixes p::"real poly"
assumes "⋀x::complex. poly (map_poly of_real p) x = 0 ⟹ Re x ≤ 0"
and "lead_coeff p = 1"
shows "coeff p i ≥ 0"
using assms
proof (induction p arbitrary: i rule: real_poly_roots_induct)
case (1 p x)
interpret map_poly_idom_hom complex_of_real ..
have h: "⋀ i. 0 ≤ coeff p i"
apply (rule 1(1))
using 1(2) apply (metis lambda_zero of_real_poly_map_mult poly_mult)
using 1(3) apply (metis lead_coeff_1 lead_coeff_mult lead_coeff_pCons(1)
mult_cancel_right2 pCons_one zero_neq_one)
done
have "x ≤ 0"
apply (subst Re_complex_of_real[symmetric])
apply (rule 1(2))
apply (subst hom_mult)
by (auto)
thus ?case
apply (cases i)
subgoal using h[of i] h[of "i-1"]
by (fastforce simp: coeff_pCons mult_nonneg_nonpos2)
subgoal using h[of i] h[of "i-1"] mult_left_mono_neg
by (fastforce simp: coeff_pCons)
done
next
case (2 p a b)
interpret map_poly_idom_hom complex_of_real ..
have h: "⋀ i. 0 ≤ coeff p i"
apply (rule 2(2))
using 2(3) apply (metis lambda_zero of_real_poly_map_mult poly_mult)
using 2(4) apply (metis lead_coeff_1 lead_coeff_mult lead_coeff_pCons(1)
mult_cancel_right2 pCons_one zero_neq_one)
done
have "Re (a + b * 𝗂) ≤ 0"
apply (rule 2(3))
apply (subst hom_mult)
by (auto simp: algebra_simps)
hence 1: "0 ≤ - 2 * a" by fastforce
have 2: "0 ≤ a * a + b * b" by fastforce
have "⋀ x. 0 ≤ coeff [:a * a + b * b, - 2 * a, 1:] x"
proof -
fix x
show "0 ≤ coeff [:a * a + b * b, - 2 * a, 1:] x"
using 2 apply (cases "x = 0", fastforce)
using 1 apply (cases "x = 1", fastforce)
apply (cases "x = 2", fastforce simp: coeff_pCons)
by (auto simp: coeff_eq_0)
qed
thus ?case
apply (subst mult.commute, subst coeff_mult)
apply (rule sum_nonneg, rule mult_nonneg_nonneg[OF _ h])
by auto
next
case (3 a)
then show ?case
by (smt (z3) bot_nat_0.extremum_uniqueI degree_1 le_degree
lead_coeff_pCons(2) pCons_one)
qed
definition circle_01_diam :: "complex set" where
"circle_01_diam =
{x. cmod (x - (of_nat 1 :: complex)/(of_nat 2)) < (real 1)/(real 2)}"
lemma pos_real_map:
"{x::complex. 1 / x ∈ (λx. x + 1) ` {x. 0 < Re x}} = circle_01_diam"
proof
show "{x. 1 / x ∈ (λx. x + 1) ` {x. 0 < Re x}} ⊆ circle_01_diam"
proof
fix x assume "x ∈ {x. 1 / x ∈ (λx. x + 1) ` {x. 0 < Re x}}"
then obtain y where h: "1 / x = y + 1" and hy: "0 < Re y" by blast
hence hy': "y = 1 / x - 1" by fastforce
hence hy'': "y + 1 ≠ 0" using h hy by fastforce
hence "x = 1 / (y + 1)" using h
by (metis div_by_1 divide_divide_eq_right mult.left_neutral)
have "¦Re y - 1¦ < ¦Re y + 1¦" using hy by simp
hence "cmod (y - 1) < cmod (y + 1)"
by (smt (z3) cmod_Re_le_iff minus_complex.simps(1) minus_complex.simps(2)
one_complex.simps plus_complex.simps(1) plus_complex.simps(2))
hence "cmod ((y - 1)/(y + 1)) < 1"
by (smt (verit) divide_less_eq_1_pos nonzero_norm_divide zero_less_norm_iff)
thus "x ∈ circle_01_diam" using hy' hy''
by (auto simp: field_simps norm_minus_commute circle_01_diam_def)
qed
show "circle_01_diam ⊆ {x. 1 / x ∈ (λx. x + 1) ` {x. 0 < Re x}}"
proof
fix x assume "x ∈ circle_01_diam"
hence "cmod (x - 1 / 2) * 2 < 1" by (auto simp: circle_01_diam_def)
hence h: "x ≠ 0" and "cmod (x - 1 / 2) * cmod 2 < 1" by auto
hence "cmod (2*x - 1) < 1"
by (smt (verit) dbl_simps(3) dbl_simps(5) div_self times_divide_eq_left
left_diff_distrib_numeral mult.commute mult_numeral_1
norm_eq_zero norm_mult norm_numeral norm_one numeral_One)
hence "cmod (((1/x - 1) - 1)/((1/x - 1) + 1)) < 1"
by (auto simp: divide_simps norm_minus_commute)
hence "cmod (((1/x - 1) - 1)/ cmod ((1/x - 1) + 1)) < 1"
by (metis (no_types, lifting) abs_norm_cancel norm_divide norm_of_real)
hence "cmod ((1/x - 1) - 1) < cmod ((1/x - 1) + 1)" using h
by (smt (verit) diff_add_cancel divide_eq_0_iff divide_less_eq_1_pos
norm_divide norm_of_real zero_less_norm_iff zero_neq_one)
hence "¦Re (1/x - 1) - 1¦ < ¦Re (1/x - 1) + 1¦"
by (smt (z3) cmod_Re_le_iff minus_complex.simps(1) minus_complex.simps(2)
one_complex.simps plus_complex.simps(1) plus_complex.simps(2))
hence "0 < Re (1/x - 1)" by linarith
moreover have "1 / x = (1/x - 1) + 1" by simp
ultimately have "0 < Re (1/x - 1) ∧ 1 / x = (1/x - 1) + 1" by blast
hence "∃xa. 0 < Re xa ∧ 1 / x = xa + 1" by blast
thus "x ∈ {x. 1 / x ∈ (λx. x + 1) ` {x. 0 < Re x}}" by blast
qed
qed
lemma one_circle_01: fixes P::"real poly" assumes hP: "degree P ≤ p" and "P ≠ 0"
and "proots_count (map_poly of_real P) circle_01_diam = 0"
shows "Bernstein_changes_01 p P = 0"
proof -
let ?Q = "(reciprocal_poly p P) ∘⇩p [:1, 1:]"
have hQ: "?Q ≠ 0"
using assms
by (simp add: Missing_Polynomial.pcompose_eq_0 reciprocal_0_iff)
hence 1: "changes (coeffs ?Q) ≥ proots_count ?Q {x. 0 < x} ∧
even (changes (coeffs ?Q) - proots_count ?Q {x. 0 < x})"
by (rule descartes_sign)
have hdeg: "degree (map_poly complex_of_real P) ≤ p"
by (rule le_trans, rule degree_map_poly_le, auto simp: assms)
have hx: "⋀x. 1 + x = 0 ⟹ 0 < Re x ⟹ False"
proof -
fix x::complex assume "1 + x = 0"
hence "x = -1" by algebra
thus "0 < Re x ⟹ False" by auto
qed
have 2: "proots_count (map_poly of_real P) circle_01_diam =
proots_count (map_poly of_real ?Q) {x. 0 < Re x}"
apply (subst pos_real_map[symmetric])
apply (subst of_real_hom.map_poly_pcompose)
apply (subst map_poly_reciprocal) using assms apply auto[2]
apply (subst proots_pcompose)
using assms apply (auto simp: reciprocal_0_iff degree_map_poly)[2]
apply (subst proots_count_reciprocal)
using assms apply (auto simp: degree_map_poly inverse_eq_divide)[2]
using hx apply fastforce
by (auto simp: inverse_eq_divide algebra_simps)
hence 3:"proots_count (map_poly of_real ?Q) {x. 0 < Re x} = 0"
using assms(3) by presburger
hence "⋀x::complex.
poly (map_poly of_real (smult (inverse (lead_coeff ?Q)) ?Q)) x = 0 ⟹
Re x ≤ 0"
proof cases
fix x::complex show "Re x ≤ 0 ⟹ Re x ≤ 0" by fast
assume "¬Re x ≤ 0" hence h:"0 < Re x" by simp
assume "poly (map_poly of_real (smult (inverse (lead_coeff ?Q)) ?Q)) x = 0"
hence h2:"poly (map_poly of_real ?Q) x = 0" by fastforce
hence "order x (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])) > 0"
using assms by (fastforce simp: order_root pcompose_eq_0 reciprocal_0_iff)
hence "proots_count (map_poly of_real ?Q) {x. 0 < Re x} ≠ 0"
proof -
have h3: "finite {x. poly (map_poly complex_of_real
(reciprocal_poly p P ∘⇩p [:1, 1:])) x = 0}"
apply (rule poly_roots_finite)
using assms by (fastforce simp: order_root pcompose_eq_0 reciprocal_0_iff)
have "0 < order x (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:]))"
using h2 assms by (fastforce simp: order_root pcompose_eq_0 reciprocal_0_iff)
also have "... ≤ (∑r∈{x. 0 < Re x ∧
poly (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])) x =
0}.
order r (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])))"
apply (rule member_le_sum) using h h2 h3 by auto
finally have
"0 < (∑r∈{x. 0 < Re x ∧
poly (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])) x = 0}.
order r (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])))" .
thus
"0 < order x (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:])) ⟹
proots_count (map_poly complex_of_real (reciprocal_poly p P ∘⇩p [:1, 1:]))
{x. 0 < Re x} ≠ 0 "
by (auto simp: proots_count_def proots_within_def)
qed
thus "Re x ≤ 0" using 3 by blast
qed
hence "⋀i. coeff (smult (inverse (lead_coeff ?Q)) ?Q) i ≥ 0"
apply (frule descartes_sign_zero)
using assms by (fastforce simp: pcompose_eq_0 reciprocal_0_iff)
hence "changes (coeffs (smult (inverse (lead_coeff ?Q)) ?Q)) = 0"
by (subst changes_all_nonneg, auto simp: nth_default_coeffs_eq)
hence "changes (coeffs ?Q) = 0"
using hQ by (auto simp: coeffs_smult changes_scale_const)
thus ?thesis
apply (subst Bernstein_changes_01_eq_changes["OF" hP])
by blast
qed
definition circle_diam :: "real ⇒ real ⇒ complex set" where
"circle_diam l r = {x. cmod ((x - l) - (r - l)/2) < (r - l)/2}"
lemma circle_diam_rescale: assumes "l < r"
shows "circle_diam l r = (λ x . (x*(r - l) + l)) ` circle_01_diam"
proof
show "circle_diam l r ⊆ (λx. x * (complex_of_real r - complex_of_real l) +
complex_of_real l) ` circle_01_diam"
proof
fix x assume "x ∈ circle_diam l r"
hence "cmod ((x - l) - (r - l)/2) < (r - l)/2" by (auto simp: circle_diam_def)
hence "cmod ((r - l) * ((x - l)/(r - l) - 1/2)) < (r - l)/2" using assms
by (subst right_diff_distrib, fastforce)
hence "(r - l) * cmod ((x - l)/(r - l) - 1/2) < (r - l) * 1/2"
apply (subst(2) abs_of_pos[symmetric])
subgoal using assms by argo
subgoal
apply (subst norm_scaleR[symmetric])
by (simp add: scaleR_conv_of_real)
done
hence "cmod ((x - l)/(r - l) - 1/2) < 1/2"
apply (subst mult_less_cancel_left_pos[of "r-l",symmetric])
using assms by auto
hence
"cmod ((x-l)/(r-l) - 1 / 2) * 2 < 1 ∧
x = (x-l)/(r-l) * (complex_of_real r - complex_of_real l) + complex_of_real l"
by force
thus "x ∈ (λx. x * (complex_of_real r - complex_of_real l) + complex_of_real l) `
circle_01_diam"
by (force simp: circle_01_diam_def)
qed
show "(λx. x * (complex_of_real r - complex_of_real l) + complex_of_real l) `
circle_01_diam ⊆ circle_diam l r"
proof
fix x::complex
assume
"x ∈ (λx. x * (complex_of_real r - complex_of_real l) + complex_of_real l) `
circle_01_diam"
then obtain y::complex where "x = y * (r - l) + l" "cmod (y - 1/2) < 1/2"
by (fastforce simp: circle_01_diam_def)
moreover hence "y = (x - l) / (r - l)" using assms by force
ultimately have "cmod ((x - l) / (r - l) - 1/2) < 1/2" by presburger
hence "(r - l) * (cmod ((x - l) / (r - l) - 1/2)) < (r - l) * (1/2)"
apply (subst mult_less_cancel_left_pos)
using assms by auto
hence "cmod ((x - l) - (r - l)/2) < (r - l)/2"
apply (subst(asm) (2) abs_of_pos[symmetric])
using assms apply argo
apply (subst(asm) norm_scaleR[symmetric])
by (smt (verit, del_insts)
‹x = y * complex_of_real (r - l) + complex_of_real l›
‹y = (x - complex_of_real l) / complex_of_real (r - l)›
add_diff_cancel divide_divide_eq_right divide_numeral_1 mult.commute
of_real_1 of_real_add of_real_divide one_add_one scaleR_conv_of_real
scale_right_diff_distrib times_divide_eq_right)
thus "x ∈ circle_diam l r"
by (force simp: circle_diam_def)
qed
qed
lemma one_circle: fixes P::"real poly" assumes "l < r"
and "proots_count (map_poly of_real P) (circle_diam l r) = 0"
and "P ≠ 0"
and "degree P ≤ p"
shows "Bernstein_changes p l r P = 0"
proof (subst Bernstein_changes_eq_rescale)
show "l ≠ r" using assms(1) by force
show "degree P ≤ p" using assms(4) by blast
show "Bernstein_changes_01 p (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:]) = 0"
proof (rule one_circle_01)
show "degree (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:]) ≤ p"
using assms(4) by (force simp: degree_pcompose)
show "P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:] ≠ 0"
using assms by (smt (z3) degree_0_iff gr_zeroI pCons_eq_0_iff pCons_eq_iff
pcompose_eq_0)
have "proots_count (map_poly of_real P) (circle_diam l r) =
proots_count (map_poly complex_of_real (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:]))
circle_01_diam"
apply (subst of_real_hom.map_poly_pcompose)
apply (subst proots_pcompose)
apply (metis assms(3) degree_eq_zeroE of_real_poly_eq_0_iff
pCons_eq_iff pCons_one pcompose_eq_0 zero_neq_one)
using assms(1) apply fastforce
apply (subst of_real_hom.map_poly_pcompose)
apply (subst proots_pcompose)
apply (auto simp: assms(3))[2]
apply (subst circle_diam_rescale[OF assms(1)])
apply (rule arg_cong[of _ _ "proots_count (map_poly complex_of_real P)"])
by fastforce
thus "proots_count (map_poly complex_of_real (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:]))
circle_01_diam = 0"
using assms(2) by presburger
qed
qed
subsection ‹One sign change case›
definition upper_circle_01 :: "complex set" where
"upper_circle_01 = {x. cmod (x - (1/2 + sqrt(3)/6 * 𝗂)) < sqrt 3 / 3}"
lemma upper_circle_map:
"{x::complex. 1 / x ∈ (λx. x + 1) ` {x. Im x < sqrt 3 * Re x}} = upper_circle_01"
proof
show "{x::complex. 1 / x ∈ (λx. x + 1) ` {x. Im x < sqrt 3 * Re x}} ⊆ upper_circle_01"
proof
fix x
assume "x ∈ {x. 1 / x ∈ (λx. x + 1) ` {x. Im x < sqrt 3 * Re x}}"
then obtain y where "1 / x = y + 1" and h: "Im y < sqrt 3 * Re y" by fastforce
hence hy: "y = 1/x - 1" by simp
hence hx: "x = 1/(y+1)" by auto
from h have hy1: "y ≠ -1" by fastforce
hence hx0: "x ≠ 0" using hy by fastforce
from h have "0 < Re ((𝗂 + sqrt 3) * y)" by fastforce
hence "cmod ((𝗂 + sqrt 3) * y - 1) < cmod ((𝗂 + sqrt 3) * y + 1)"
by (auto simp: cmod_def power2_eq_square algebra_simps)
hence 1: "cmod (((𝗂 + sqrt 3) * y - 1)/((𝗂 + sqrt 3) * y + 1)) < 1"
by (auto simp: norm_divide divide_simps)
also have "cmod (((𝗂 + sqrt 3) * y - 1)/((𝗂 + sqrt 3) * y + 1)) =
cmod (((𝗂 + sqrt 3) * y * x - x)/((𝗂 + sqrt 3) * y * x + x))"
apply (subst mult_divide_mult_cancel_right[symmetric, OF hx0])
apply (subst ring_distribs(2)[of _ _ x])
apply (subst left_diff_distrib[of _ _ x])
by simp
also have "... = cmod
(((-1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂)) /
(( 1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂)))"
by (auto simp: hy algebra_simps hx0)
also have
"... = cmod ((-1 - complex_of_real (sqrt 3) - 𝗂) * x +
(complex_of_real (sqrt 3) + 𝗂)) /
cmod (( 1 - complex_of_real (sqrt 3) - 𝗂) * x +
(complex_of_real (sqrt 3) + 𝗂))"
by (auto simp: norm_divide)
finally have
"cmod ((-1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂))
< cmod ((1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂))"
proof (subst divide_less_eq_1_pos[symmetric], subst zero_less_norm_iff)
show "(1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂) ≠ 0"
proof
have "-𝗂 + 1 ≠ complex_of_real (sqrt 3)" by (auto simp: complex_eq_iff)
moreover assume
"(1 - complex_of_real (sqrt 3) - 𝗂) * x + (complex_of_real (sqrt 3) + 𝗂) = 0"
ultimately have
"x = (-complex_of_real (sqrt 3) - 𝗂)/(1 - complex_of_real (sqrt 3) - 𝗂)"
by (auto simp: divide_simps algebra_simps)
thus False
using h by (auto simp: hy field_simps Im_divide Re_divide)
qed
qed
hence "cmod (x - (1/2 + sqrt(3)/6 * 𝗂)) < sqrt 3 / 3"
apply (subst(3) abs_of_pos[symmetric, of 3]) apply auto[1]
apply (subst real_sqrt_abs2[symmetric], subst real_sqrt_divide[symmetric])
apply (subst cmod_def, subst real_sqrt_less_iff)
apply (rule mult_right_less_imp_less[of _ "sqrt 3 /3"])
by (auto simp: cmod_def power2_eq_square algebra_simps)
thus "x ∈ upper_circle_01"
by (auto simp: upper_circle_01_def)
qed
show "upper_circle_01 ⊆ {x. 1 / x ∈ (λx. x + 1) ` {x. sqrt 3 * Re x > Im x}}"
proof
fix x assume "x ∈ upper_circle_01"
hence "cmod (x - (1/2 + sqrt(3)/6 * 𝗂)) < sqrt 3 / 3"
by (force simp: upper_circle_01_def)
hence "sqrt ((Re x - 1/2)^2 + (Im x - sqrt(3)/6)^2) < sqrt (1/3)"
by (auto simp: cmod_def sqrt_divide_self_eq real_sqrt_inverse[symmetric])
hence 1: "- Im x * sqrt 3 + (Im x * (Im x * 3) + Re x * (Re x * 3)) < Re x * 3"
by (auto simp: power2_eq_square algebra_simps)
have 2: "- Im x + (Im x * (Im x * sqrt 3) + Re x * (Re x * sqrt 3)) < Re x * sqrt 3"
apply (rule mult_right_less_imp_less[of _ "sqrt 3"])
apply (subst mult.assoc[of _ "sqrt 3"], subst real_sqrt_mult_self)
using 1 by (auto simp: algebra_simps)
have "sqrt 3 + (-Im x) / (Im x * Im x + Re x * Re x) <
Re x * sqrt 3 / (Im x * Im x + Re x * Re x)"
apply (rule mult_right_less_imp_less[of _ "(Im x * Im x + Re x * Re x)"])
apply (rule subst, rule arg_cong2[of _ _ _ _ "(<)"])
prefer 3 apply (rule 2)
apply (subst distrib_right)
using 2 apply auto
by (auto simp: algebra_simps)
hence "0 < - Im (1/x-1) + sqrt 3 * Re (1/x-1)"
by (auto simp: power2_eq_square algebra_simps Re_divide Im_divide)
hence "sqrt 3 * Re (1/x-1) > Im (1/x-1)"
by argo
hence "(1/x-1) ∈ {x. sqrt 3 * Re x > Im x}" by fast
moreover have "1/x = (1/x-1) + 1" by simp
ultimately show "x ∈ {x. 1 / x ∈ (λx. x + 1) ` {x. sqrt 3 * Re x > Im x}}"
by blast
qed
qed
definition lower_circle_01 :: "complex set" where
"lower_circle_01 = {x. cmod (x - (1/2 - sqrt(3)/6 * 𝗂)) < sqrt 3 / 3}"
lemma cnj_upper_circle_01: "cnj ` upper_circle_01 = lower_circle_01"
proof
show "cnj ` upper_circle_01 ⊆ lower_circle_01"
proof
fix x assume "x ∈ cnj ` upper_circle_01"
then obtain y where "y ∈ upper_circle_01" and "x = cnj y" by blast
thus "x ∈ lower_circle_01"
apply (subst lower_circle_01_def, subst complex_mod_cnj[symmetric])
by (auto simp add: upper_circle_01_def del: complex_mod_cnj)
qed
show "lower_circle_01 ⊆ cnj ` upper_circle_01"
proof
fix x assume "x ∈ lower_circle_01"
hence "cnj x ∈ upper_circle_01" and "x = cnj (cnj x)"
apply (subst upper_circle_01_def, subst complex_mod_cnj[symmetric])
by (auto simp add: lower_circle_01_def del: complex_mod_cnj)
thus "x ∈ cnj ` upper_circle_01"
by blast
qed
qed
lemma lower_circle_map:
"{x::complex. 1 / x ∈ (λx. x + 1) ` {x. Im x > - sqrt 3 * Re x}} = lower_circle_01"
proof (subst cnj_upper_circle_01[symmetric], subst upper_circle_map[symmetric])
have "{x. 1 / x ∈ (λx. x + 1) ` {x. - sqrt 3 * Re x < Im x}} =
{x. 1 / x ∈ (λx. x + 1) ` {x. sqrt 3 * Re (cnj x) > Im (cnj x)}}"
by auto
also have "... = {x. 1 / x ∈ (λx. x + 1) ` cnj ` {x. sqrt 3 * Re x > Im x}}"
apply (subst(2) bij_image_Collect_eq)
apply (metis bijI' complex_cnj_cnj)
by (auto simp: inj_def inj_imp_inv_eq[of _ cnj])
also have "... = {x. 1 / x ∈ cnj ` (λx. x + 1) ` {x. sqrt 3 * Re x > Im x}}"
by (auto simp: image_image)
also have "... = {x. cnj (1 / x) ∈ (λx. x + 1) ` {x. sqrt 3 * Re x > Im x}}"
by (metis (no_types, lifting) complex_cnj_cnj image_iff)
also have "... = cnj ` {x. 1 / x ∈ (λx. x + 1) ` {x. sqrt 3 * Re x > Im x}}"
apply (subst(2) bij_image_Collect_eq)
apply (metis bijI' complex_cnj_cnj)
by (auto simp: inj_def inj_imp_inv_eq[of _ cnj])
finally show "{x. 1 / x ∈ (λx. x + 1) ` {x. - sqrt 3 * Re x < Im x}} =
cnj ` {x. 1 / x ∈ (λx. x + 1) ` {x. Im x < sqrt 3 * Re x}}" .
qed
lemma two_circles_01:
fixes P::"real poly"
assumes hP: "degree P ≤ p" and hP0: "P ≠ 0" and hp0: "p ≠ 0"
and h: "proots_count (map_poly of_real P)
(upper_circle_01 ∪ lower_circle_01) = 1"
shows "Bernstein_changes_01 p P = 1"
proof (subst Bernstein_changes_01_eq_changes[OF hP])
let ?Q = "reciprocal_poly p P ∘⇩p [:1, 1:]"
have hQ0: "?Q ≠ 0" using hP0
by (simp add: pcompose_eq_0 hP reciprocal_0_iff)
from h obtain x' where hroot': "poly (map_poly of_real P) x' = 0"
and hx':"x' ∈ upper_circle_01 ∪ lower_circle_01"
using proots_count_pos by (metis less_numeral_extra(1))
obtain x where hxx': "x' = complex_of_real x"
proof (cases "Im x' = 0")
assume "Im x' = 0" and h: "⋀x. x' = complex_of_real x ⟹ thesis"
then show thesis using h[of "Re x'"] by (simp add: complex_is_Real_iff)
next
assume hx'': "Im x' ≠ 0"
have 1: "card {x', cnj x'} = 2"
proof (subst card_2_iff)
from hx'' have "x' ≠ cnj x'" and "{x', cnj x'} = {x', cnj x'}"
by (metis cnj.simps(2) neg_equal_zero, argo)
thus "∃x y. {x', cnj x'} = {x, y} ∧ x ≠ y" by blast
qed
moreover have "{x', cnj x'} ⊆ upper_circle_01 ∪ lower_circle_01" using hx'
apply (rule UnE)
by (auto simp: cnj_upper_circle_01[symmetric])
moreover have "⋀x. x ∈ {x', cnj x'} ⟹ poly (map_poly of_real P) x = 0"
using hroot' poly_map_poly_of_real_cnj by auto
ultimately have
"proots_count (map_poly of_real P) (upper_circle_01 ∪ lower_circle_01) ≥ 2"
apply (subst 1[symmetric])
apply (rule proots_count_of_root_set)
using assms(2) of_real_poly_eq_0_iff by (blast, blast, blast)
thus thesis using assms(4) by linarith
qed
hence hroot: "poly P x = 0"
using hroot' by (metis of_real_0 of_real_eq_iff of_real_poly_map_poly)
have that: "3 * sqrt (x * x + 1 / 3 - x) < sqrt 3" using hx'
apply (rule UnE)
by (auto simp: cmod_def power2_eq_square algebra_simps upper_circle_01_def
lower_circle_01_def hxx')
have hx: "0 < x ∧ x < 1"
proof -
have "3 * sqrt (x * x + 1 / 3 - x) = sqrt (9 * (x * x + 1 / 3 - x))"
by (subst real_sqrt_mult, simp)
hence "9 * (x * x + 1 / 3 - x) < 3" using that real_sqrt_less_iff by metis
hence "x*x - x < 0" by auto
thus "0 < x ∧ x < 1"
using mult_eq_0_iff mult_less_cancel_right_disj by fastforce
qed
let ?y = "1/x - 1"
from hroot hx assms have "poly ?Q ?y = 0"
by (auto simp: poly_pcompose poly_reciprocal)
hence "[:-?y, 1:] dvd ?Q" using poly_eq_0_iff_dvd by blast
then obtain R where hR: "?Q = R * [:-?y, 1:]" by auto
hence hR0: "R ≠ 0" using hQ0 by force
interpret map_poly_idom_hom complex_of_real ..
have "normal_poly (smult (inverse (lead_coeff (reciprocal_poly p P ∘⇩p [:1, 1:]))) R)"
proof (rule normal_poly_of_roots)
show "⋀z. poly (map_poly complex_of_real
(smult (inverse (lead_coeff (reciprocal_poly p P ∘⇩p [:1, 1:]))) R)) z = 0 ⟹
Re z ≤ 0 ∧ (cmod z)⇧2 ≤ 4 * (Re z)⇧2"
proof -
fix z
assume
"poly (map_poly complex_of_real
(smult (inverse (lead_coeff (reciprocal_poly p P ∘⇩p [:1, 1:]))) R)) z = 0"
hence hroot2: "poly (map_poly complex_of_real R) z = 0"
by (auto simp: map_poly_smult hQ0)
show "Re z ≤ 0 ∧ (cmod z)⇧2 ≤ 4 * (Re z)⇧2"
proof (rule ccontr)
assume "¬ (Re z ≤ 0 ∧ (cmod z)⇧2 ≤ 4 * (Re z)⇧2)"
hence 1: "0 < Re z ∨ 4 * (Re z)⇧2 < (cmod z)⇧2" by linarith
hence hz: "z ≠ -1" by force
have "Im z > - sqrt 3 * Re z ∨ sqrt 3 * Re z > Im z"
proof (cases "Im z ≥ sqrt 3 * Re z", cases "- sqrt 3 * Re z ≥ Im z")
assume 2: "sqrt 3 * Re z ≤ Im z" "Im z ≤ - sqrt 3 * Re z"
hence "sqrt 3 * Re z ≤ sqrt 3 * - Re z" by force
hence "Re z ≤ - Re z"
apply (rule mult_left_le_imp_le)
by fastforce
hence "Re z ≤ 0" by simp
moreover have "(Im z)^2 ≤ (-sqrt 3 * Re z)^2"
apply (subst power2_eq_square, subst power2_eq_square)
apply (rule square_bounded_le)
using 2 by auto
ultimately have False using 1
by (auto simp: power2_eq_square cmod_def algebra_simps)
thus ?thesis by fast
qed auto
hence "z ∈ {z. - sqrt 3 * Re z < Im z} ∪ {z. Im z < sqrt 3 * Re z}"
by blast
hence 1: "inverse (1 + z) ∈ upper_circle_01 ∪ lower_circle_01"
by (force simp: inverse_eq_divide upper_circle_map[symmetric]
lower_circle_map[symmetric])
have hRdeg': "degree R < p"
apply (rule less_le_trans[of "degree R" "degree ?Q"])
apply (subst hR, subst degree_mult_eq[OF hR0], fastforce, fastforce)
by (auto simp: degree_pcompose degree_reciprocal hP)
hence hRdeg: "degree R ≤ p" by fastforce
have 2: "map_poly complex_of_real (reciprocal_poly p (R ∘⇩p [:-1, 1:])) ≠ 0"
apply (subst of_real_poly_eq_0_iff, subst reciprocal_0_iff)
apply (force simp: hRdeg degree_pcompose)
using hR0 pcompose_eq_0
by (metis degree_eq_zeroE gr_zeroI pCons_eq_iff pCons_one zero_neq_one)
have 3:
"poly (map_poly complex_of_real (reciprocal_poly p (R ∘⇩p [:-1, 1:])))
(inverse (1 + z)) = 0"
apply (subst map_poly_reciprocal)
using hRdeg apply (force simp: degree_pcompose)
apply auto[1]
apply (subst poly_reciprocal)
using hRdeg apply (force simp: degree_map_poly degree_pcompose)
using hz apply (metis inverse_nonzero_iff_nonzero neg_eq_iff_add_eq_0)
by (auto simp: of_real_hom.map_poly_pcompose poly_pcompose hroot2)
have "proots_count (map_poly of_real (reciprocal_poly p (R ∘⇩p [:-1, 1:])))
(upper_circle_01 ∪ lower_circle_01) > 0"
by (rule proots_count_of_root[OF 2 1 3])
moreover have "proots_count
(map_poly complex_of_real
(reciprocal_poly p ([:1 - 1 / x, 1:] ∘⇩p [:- 1, 1:])))
(upper_circle_01 ∪ lower_circle_01) > 0"
apply (subst map_poly_reciprocal)
using hp0 less_eq_Suc_le apply (simp add: degree_pcompose)
apply simp
apply (subst proots_count_reciprocal)
using hp0 less_eq_Suc_le
apply (simp add: degree_pcompose degree_map_poly)
apply (auto simp: pcompose_pCons)[1]
apply (auto simp: cmod_def power2_eq_square real_sqrt_divide
real_div_sqrt upper_circle_01_def lower_circle_01_def)[1]
apply (subst of_real_hom.map_poly_pcompose)
apply (subst proots_pcompose, fastforce, force)
apply (subst lower_circle_map[symmetric])
apply (subst upper_circle_map[symmetric])
apply (rule proots_count_of_root[of _ "of_real (1/x - 1)"])
apply simp
apply (auto simp: bij_image_Collect_eq bij_def inj_def image_iff
inverse_eq_divide inj_imp_inv_eq[of _ "λ x. x+1"] hx simp del: surj_plus_right)[1]
by force
ultimately have "proots_count
(map_poly complex_of_real (reciprocal_poly p (R ∘⇩p [:- 1, 1:])))
(upper_circle_01 ∪ lower_circle_01) +
proots_count
(map_poly complex_of_real
(reciprocal_poly p ([:1 - 1 / x, 1:] ∘⇩p [:- 1, 1:])))
(upper_circle_01 ∪ lower_circle_01) > 1"
by fastforce
also have "... = proots_count (map_poly complex_of_real
(monom 1 p * reciprocal_poly p (?Q ∘⇩p [:- 1, 1:])))
(upper_circle_01 ∪ lower_circle_01)"
apply (subst eq_commute, subst hR, subst pcompose_mult)
apply (subst reciprocal_mult, subst degree_mult_eq)
using hR0 apply (fastforce simp: pcompose_eq_0)
apply (fastforce simp: pcompose_pCons)
using hRdeg' apply (simp add: degree_pcompose)
using hRdeg apply (simp add: degree_pcompose)
using hp0 apply (auto simp: degree_pcompose)[1]
apply (subst hom_mult)
apply (subst proots_count_times)
using hp0 hRdeg' hR0
apply (fastforce simp: reciprocal_0_iff degree_pcompose pcompose_eq_0
pcompose_pCons)
by simp
also have "... = proots_count
(map_poly complex_of_real
(reciprocal_poly p (reciprocal_poly p P ∘⇩p [:1, 1:] ∘⇩p [:- 1, 1:])))
(upper_circle_01 ∪ lower_circle_01)"
apply (subst hom_mult)
apply (subst proots_count_times)
using hp0 hP hP0
apply (auto simp: map_poly_reciprocal degree_pcompose
degree_reciprocal of_real_hom.map_poly_pcompose
reciprocal_0_iff degree_map_poly pcompose_eq_0)[1]
apply (subst map_poly_monom, fastforce)
apply (subst of_real_1, subst proots_count_monom)
apply (auto simp: cmod_def power2_eq_square real_sqrt_divide
real_div_sqrt upper_circle_01_def lower_circle_01_def)[1]
by presburger
also have "... = 1"
by (auto simp: pcompose_assoc["symmetric"] pcompose_pCons
reciprocal_reciprocal hP h)
finally show False by blast
qed
qed
show "lead_coeff
(smult (inverse (lead_coeff (reciprocal_poly p P ∘⇩p [:1, 1:]))) R) = 1"
by (auto simp: hR degree_add_eq_right hR0 coeff_eq_0)
qed
hence "changes (coeffs (smult (inverse (lead_coeff ?Q)) ?Q)) = 1"
apply (subst hR, subst mult_smult_left[symmetric], rule normal_changes)
by (auto simp: hx)
moreover have "changes (coeffs (reciprocal_poly p P ∘⇩p [:1, 1:])) =
changes (coeffs (smult (inverse (lead_coeff (reciprocal_poly p P ∘⇩p [:1, 1:])))
(reciprocal_poly p P ∘⇩p [:1, 1:])))"
by (auto simp: pcompose_eq_0 reciprocal_0_iff hP hP0 coeffs_smult
changes_scale_const[symmetric])
ultimately show "changes (coeffs (reciprocal_poly p P ∘⇩p [:1, 1:])) = 1" by argo
qed
definition upper_circle :: "real ⇒ real ⇒ complex set" where
"upper_circle l r = {x::complex.
cmod ((x-of_real l)/(of_real (r-l)) - (1/2 + of_real (sqrt(3))/6 * 𝗂)) < sqrt 3 / 3}"
lemma upper_circle_rescale: assumes "l < r"
shows "upper_circle l r = (λ x . (x*(r - l) + l)) ` upper_circle_01"
proof
show "upper_circle l r ⊆
(λx. x * (of_real r - of_real l) + of_real l) ` upper_circle_01"
apply (rule subsetI)
apply (rule image_eqI[of _ _ "(_ - of_real l)/(of_real r - of_real l)"])
using assms apply (auto simp: divide_simps)[1]
by (auto simp: upper_circle_01_def upper_circle_def)
show "(λx. x * (of_real r - of_real l) + of_real l) ` upper_circle_01 ⊆
upper_circle l r"
apply (rule subsetI, subst(asm) image_iff)
using assms by (auto simp: upper_circle_01_def upper_circle_def)
qed
definition lower_circle :: "real ⇒ real ⇒ complex set" where
"lower_circle l r = {x::complex.
cmod ((x-of_real l)/(of_real (r-l)) - (1/2 - of_real (sqrt(3))/6 * 𝗂)) < sqrt 3 / 3}"
lemma lower_circle_rescale:
assumes "l < r"
shows "lower_circle l r = (λ x . (x*(r - l) + l)) ` lower_circle_01"
proof
show "lower_circle l r ⊆ (λx. x * (of_real r - of_real l) + of_real l) `
lower_circle_01"
apply (rule subsetI)
apply (rule image_eqI[of _ _ "(_ - of_real l)/(of_real r - of_real l)"])
using assms apply (auto simp: divide_simps)[1]
by (auto simp: lower_circle_01_def lower_circle_def)
show "(λx. x * (of_real r - of_real l) + of_real l) ` lower_circle_01 ⊆
lower_circle l r"
apply (rule subsetI, subst(asm) image_iff)
using assms by (auto simp: lower_circle_01_def lower_circle_def)
qed
lemma two_circles:
fixes P::"real poly" and l r::real
assumes hlr: "l < r"
and hP: "degree P ≤ p"
and hP0: "P ≠ 0"
and hp0: "p ≠ 0"
and h: "proots_count (map_poly of_real P)
(upper_circle l r ∪ lower_circle l r) = 1"
shows "Bernstein_changes p l r P = 1"
proof -
have 1: "Bernstein_changes p l r P =
Bernstein_changes_01 p (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:])"
using assms by (force simp: Bernstein_changes_eq_rescale)
have "proots_count (map_poly complex_of_real (P ∘⇩p [:l, 1:] ∘⇩p [:0, r - l:]))
(upper_circle_01 ∪ lower_circle_01) = 1"
using assms
by (auto simp: upper_circle_rescale lower_circle_rescale proots_pcompose image_Un
of_real_hom.map_poly_pcompose pcompose_eq_0 image_image algebra_simps)
thus ?thesis
apply (subst 1)
apply (rule two_circles_01)
using hP apply (force simp: degree_pcompose)
using hP0 hlr apply (fastforce simp: pcompose_eq_0)
using hp0 apply blast
by blast
qed
subsection ‹The theorem of three circles›
theorem three_circles:
fixes P::"real poly" and l r::real
assumes "l < r"
and hP: "degree P ≤ p"
and hP0: "P ≠ 0"
and hp0: "p ≠ 0"
shows "proots_count (map_poly of_real P) (circle_diam l r) = 0 ⟹
Bernstein_changes p l r P = 0"
and "proots_count (map_poly of_real P)
(upper_circle l r ∪ lower_circle l r) = 1 ⟹
Bernstein_changes p l r P = 1"
apply (rule one_circle)
using assms apply auto[4]
apply (rule two_circles)
using assms by auto
end