Theory Cauchy_Root_Bound
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))
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