(* Author: René Thiemann Akihisa Yamada License: BSD *) section ‹Cauchy's Root Bound› text ‹This theory contains a formalization of Cauchy's root bound, i.e., given an integer polynomial it determines a bound ‹b› such that all real or complex roots of the polynomials have a norm below ‹b›.› theory Cauchy_Root_Bound imports Algebraic_Numbers_Pre_Impl begin hide_const (open) UnivPoly.coeff hide_const (open) Module.smult text ‹Division of integers, rounding to the upper value.› definition div_ceiling :: "int ⇒ int ⇒ int" where "div_ceiling x y = (let q = x div y in if q * y = x then q else q + 1)" definition root_bound :: "int poly ⇒ rat" where "root_bound p ≡ let n = degree p; m = 1 + div_ceiling (max_list_non_empty (map (λi. abs (coeff p i)) [0..<n])) (abs (lead_coeff p)) ― ‹round to the next higher number ‹2^n›, so that bisection will› ― ‹stay on integers for as long as possible› in of_int (2 ^ (log_ceiling 2 m))" lemma root_imp_deg_nonzero: assumes "p ≠ 0" "poly p x = 0" shows "degree p ≠ 0" proof assume "degree p = 0" from degree0_coeffs[OF this] assms show False by auto qed lemma cauchy_root_bound: fixes x :: "'a :: real_normed_field" assumes x: "poly p x = 0" and p: "p ≠ 0" shows "norm x ≤ 1 + max_list_non_empty (map (λ i. norm (coeff p i)) [0 ..< degree p]) / norm (lead_coeff p)" (is "_ ≤ _ + ?max / ?nlc") proof - let ?n = "degree p" let ?p = "coeff p" let ?lc = "lead_coeff p" define ml where "ml = ?max / ?nlc" from p have lc: "?lc ≠ 0" by auto hence nlc: "norm ?lc > 0" by auto from root_imp_deg_nonzero[OF p x] have *: "0 ∈ set [0 ..< degree p]" by auto have "0 ≤ norm (?p 0)" by simp also have "… ≤ ?max" by (rule max_list_non_empty, insert *, auto) finally have max0: "?max ≥ 0" . with nlc have ml0: "ml ≥ 0" unfolding ml_def by auto hence easy: "norm x ≤ 1 ⟹ ?thesis" unfolding ml_def[symmetric] by auto show ?thesis proof (cases "norm x ≤ 1") case True thus ?thesis using easy by auto next case False hence nx: "norm x > 1" by simp hence x0: "x ≠ 0" by auto hence xn0: "0 < norm x ^ ?n" by auto from x[unfolded poly_altdef] have "x ^ ?n * ?lc = x ^ ?n * ?lc - (∑i≤?n. x ^ i * ?p i)" unfolding poly_altdef by (simp add: ac_simps) also have "(∑i≤?n. x ^ i * ?p i) = x ^ ?n * ?lc + (∑i < ?n. x ^ i * ?p i)" by (subst sum.remove[of _ ?n], auto intro: sum.cong) finally have "x ^ ?n * ?lc = - (∑i < ?n. x ^ i * ?p i)" by simp with lc have "x ^ ?n = - (∑i < ?n. x ^ i * ?p i) / ?lc" by (simp add: field_simps) from arg_cong[OF this, of norm] have "norm x ^ ?n = norm ((∑i < ?n. x ^ i * ?p i) / ?lc)" unfolding norm_power by simp also have "(∑i < ?n. x ^ i * ?p i) / ?lc = (∑i < ?n. x ^ i * ?p i / ?lc)" by (rule sum_divide_distrib) also have "norm … ≤ (∑i < ?n. norm (x ^ i * (?p i / ?lc)))" by (simp add: field_simps, rule norm_sum) also have "… = (∑i < ?n. norm x ^ i * norm (?p i / ?lc))" unfolding norm_mult norm_power .. also have "… ≤ (∑i < ?n. norm x ^ i * ml)" proof (rule sum_mono) fix i assume "i ∈ {..<?n}" hence i: "i < ?n" by simp show "norm x ^ i * norm (?p i / ?lc) ≤ norm x ^ i * ml" proof (rule mult_left_mono) show "0 ≤ norm x ^ i" using nx by auto show "norm (?p i / ?lc) ≤ ml" unfolding norm_divide ml_def by (rule divide_right_mono[OF max_list_non_empty], insert nlc i, auto) qed qed also have "… = ml * (∑i < ?n. norm x ^ i)" unfolding sum_distrib_right[symmetric] by simp also have "(∑i < ?n. norm x ^ i) = (norm x ^ ?n - 1) / (norm x - 1)" by (rule geometric_sum, insert nx, auto) finally have "norm x ^ ?n ≤ ml * (norm x ^ ?n - 1) / (norm x - 1)" by simp from mult_left_mono[OF this, of "norm x - 1"] have "(norm x - 1) * (norm x ^ ?n) ≤ ml * (norm x ^ ?n - 1)" using nx by auto also have "… = (ml * (1 - 1 / (norm x ^ ?n))) * norm x ^ ?n" using nx False x0 by (simp add: field_simps) finally have "(norm x - 1) * (norm x ^ ?n) ≤ (ml * (1 - 1 / (norm x ^ ?n))) * norm x ^ ?n" . from mult_right_le_imp_le[OF this xn0] have "norm x - 1 ≤ ml * (1 - 1 / (norm x ^ ?n))" by simp hence "norm x ≤ 1 + ml - ml / (norm x ^ ?n)" by (simp add: field_simps) also have "… ≤ 1 + ml" using ml0 xn0 by auto finally show ?thesis unfolding ml_def . qed qed lemma div_le_div_ceiling: "x div y ≤ div_ceiling x y" unfolding div_ceiling_def Let_def by auto lemma div_ceiling: assumes q: "q ≠ 0" shows "(of_int x :: 'a :: floor_ceiling) / of_int q ≤ of_int (div_ceiling x q)" proof (cases "q dvd x") case True then obtain k where xqk: "x = q * k" unfolding dvd_def by auto hence id: "div_ceiling x q = k" unfolding div_ceiling_def Let_def using q by auto show ?thesis unfolding id unfolding xqk using q by simp next case False { assume "x div q * q = x" hence "x = q * (x div q)" by (simp add: ac_simps) hence "q dvd x" unfolding dvd_def by auto with False have False by simp } hence id: "div_ceiling x q = x div q + 1" unfolding div_ceiling_def Let_def using q by auto show ?thesis unfolding id by (metis floor_divide_of_int_eq le_less add1_zle_eq floor_less_iff) qed lemma max_list_non_empty_map: assumes hom: "⋀ x y. max (f x) (f y) = f (max x y)" shows "xs ≠ [] ⟹ max_list_non_empty (map f xs) = f (max_list_non_empty xs)" by (induct xs rule: max_list_non_empty.induct, auto simp: hom) lemma root_bound: assumes "root_bound p = B" and deg: "degree p > 0" shows "ipoly p (x :: 'a :: real_normed_field) = 0 ⟹ norm x ≤ of_rat B" "B ≥ 0" proof - let ?r = "of_rat :: _ ⇒ 'a" let ?i = "of_int :: _ ⇒ 'a" let ?p = "map_poly ?i p" define n where "n = degree p" let ?lc = "coeff p n" let ?list = "map (λi. abs (coeff p i)) [0..<n]" let ?list' = "(map (λi. real_of_int (abs ((coeff p i)))) [0..<n])" define m where "m = max_list_non_empty ?list" define m_up where "m_up = 1 + div_ceiling m (abs ?lc)" define C where "C = rat_of_int (2^(log_ceiling 2 m_up))" from deg have p0: "p ≠ 0" by auto from p0 have alc0: "abs ?lc ≠ 0" unfolding n_def by auto from deg have mem: "abs (coeff p 0) ∈ set ?list" unfolding n_def by auto from max_list_non_empty[OF this, folded m_def] have m0: "m ≥ 0" by auto have "div_ceiling m (abs ?lc) ≥ 0" by (rule order_trans[OF _ div_le_div_ceiling[of m "abs ?lc"]], subst pos_imp_zdiv_nonneg_iff, insert p0 m0, auto simp: n_def) hence mup: "m_up ≥ 1" unfolding m_up_def by auto have "m_up ≤ 2 ^ (log_ceiling 2 m_up)" using mup log_ceiling_sound(1) by auto hence Cmup: "C ≥ of_int m_up" unfolding C_def by linarith with mup have C: "C ≥ 1" by auto from assms(1)[unfolded root_bound_def Let_def] have B: "C = B" unfolding C_def m_up_def n_def m_def by auto note dc = div_le_div_ceiling[of m "abs ?lc"] with C show "B ≥ 0" unfolding B by auto assume "ipoly p x = 0" hence rt: "poly ?p x = 0" by simp from root_imp_deg_nonzero[OF _ this] p0 have n0: "n ≠ 0" unfolding n_def by auto from cauchy_root_bound[OF rt] p0 have "norm x ≤ 1 + max_list_non_empty ?list' / real_of_int (abs ?lc)" by (simp add: n_def) also have "?list' = map real_of_int ?list" by simp also have "max_list_non_empty … = real_of_int m" unfolding m_def by (rule max_list_non_empty_map, insert mem, auto) also have "1 + m / real_of_int (abs ?lc) ≤ real_of_int m_up" unfolding m_up_def using div_ceiling[OF alc0, of m] by auto also have "… ≤ real_of_rat C" using Cmup using of_rat_less_eq by force finally have "norm x ≤ real_of_rat C" . thus "norm x ≤ real_of_rat B" unfolding B by simp qed end