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)[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