Theory Real_Roots
section ‹Real Roots›
text ‹This theory contains an algorithm to determine the set of real roots of a
rational polynomial. For polynomials with real coefficients, we refer to
the AFP entry "Factor-Algebraic-Polynomial".›
theory Real_Roots
imports
Cauchy_Root_Bound
Real_Algebraic_Numbers
begin
hide_const (open) UnivPoly.coeff
hide_const (open) Module.smult
partial_function (tailrec) roots_of_2_main ::
"int poly ⇒ root_info ⇒ (rat ⇒ rat ⇒ nat) ⇒ (rat × rat)list ⇒ real_alg_2 list ⇒ real_alg_2 list" where
[code]: "roots_of_2_main p ri cr lrs rais = (case lrs of Nil ⇒ rais
| (l,r) # lrs ⇒ let c = cr l r in
if c = 0 then roots_of_2_main p ri cr lrs rais
else if c = 1 then roots_of_2_main p ri cr lrs (real_alg_2'' ri p l r # rais)
else let m = (l + r) / 2 in roots_of_2_main p ri cr ((m,r) # (l,m) # lrs) rais)"
definition roots_of_2_irr :: "int poly ⇒ real_alg_2 list" where
"roots_of_2_irr p = (if degree p = 1
then [Rational (Rat.Fract (- coeff p 0) (coeff p 1)) ] else
let ri = root_info p;
cr = root_info.l_r ri;
B = root_bound p
in (roots_of_2_main p ri cr [(-B,B)] []))"
fun pairwise_disjoint :: "'a set list ⇒ bool" where
"pairwise_disjoint [] = True"
| "pairwise_disjoint (x # xs) = ((x ∩ (⋃ y ∈ set xs. y) = {}) ∧ pairwise_disjoint xs)"
lemma roots_of_2_irr: assumes pc: "poly_cond p" and deg: "degree p > 0"
shows "real_of_2 ` set (roots_of_2_irr p) = {x. ipoly p x = 0}" (is ?one)
"Ball (set (roots_of_2_irr p)) invariant_2" (is ?two)
"distinct (map real_of_2 (roots_of_2_irr p))" (is ?three)
proof -
note d = roots_of_2_irr_def
from poly_condD[OF pc] have mon: "lead_coeff p > 0" and irr: "irreducible p" by auto
let ?norm = "real_alg_2'"
have "?one ∧ ?two ∧ ?three"
proof (cases "degree p = 1")
case True
define c where "c = coeff p 0"
define d where "d = coeff p 1"
from True have rr: "roots_of_2_irr p = [Rational (Rat.Fract (- c) (d))]" unfolding d d_def c_def by auto
from degree1_coeffs[OF True] have p: "p = [:c,d:]" and d: "d ≠ 0" unfolding c_def d_def by auto
have *: "real_of_int c + x * real_of_int d = 0 ⟹ x = - (real_of_int c / real_of_int d)" for x
using d by (simp add: field_simps)
show ?thesis unfolding rr using d * unfolding p using of_rat_1[of "Rat.Fract (- c) (d)"]
by (auto simp: Fract_of_int_quotient hom_distribs)
next
case False
let ?r = real_of_rat
let ?rp = "map_poly ?r"
let ?rr = "set (roots_of_2_irr p)"
define ri where "ri = root_info p"
define cr where "cr = root_info.l_r ri"
define bnds where "bnds = [(-root_bound p, root_bound p)]"
define empty where "empty = (Nil :: real_alg_2 list)"
have empty: "Ball (set empty) invariant_2 ∧ distinct (map real_of_2 empty)" unfolding empty_def by auto
from mon have p: "p ≠ 0" by auto
from root_info[OF irr deg] have ri: "root_info_cond ri p" unfolding ri_def .
from False
have rr: "roots_of_2_irr p = roots_of_2_main p ri cr bnds empty"
unfolding d ri_def cr_def Let_def bnds_def empty_def by auto
note root_bound = root_bound[OF refl deg]
from root_bound(2)
have bnds: "⋀ l r. (l,r) ∈ set bnds ⟹ l ≤ r" unfolding bnds_def by auto
have "ipoly p x = 0 ⟹ ?r (- root_bound p) ≤ x ∧ x ≤ ?r (root_bound p)" for x
using root_bound(1)[of x] by (auto simp: hom_distribs)
hence rts: "{x. ipoly p x = 0}
= real_of_2 ` set empty ∪ {x. ∃ l r. root_cond (p,l,r) x ∧ (l,r) ∈ set bnds}"
unfolding empty_def bnds_def by (force simp: root_cond_def)
define rts where "rts lr = Collect (root_cond (p,lr))" for lr
have disj: "pairwise_disjoint (real_of_2 ` set empty # map rts bnds)"
unfolding empty_def bnds_def by auto
from deg False have deg1: "degree p > 1" by auto
define delta where "delta = ipoly_root_delta p"
note delta = ipoly_root_delta[OF p, folded delta_def]
define rel' where "rel' = ({(x, y). 0 ≤ y ∧ delta_gt delta x y})^-1"
define mm where "mm = (λbnds. mset (map (λ (l,r). ?r r - ?r l) bnds))"
define rel where "rel = inv_image (mult1 rel') mm"
have wf: "wf rel" unfolding rel_def rel'_def
by (rule wf_inv_image[OF wf_mult1[OF SN_imp_wf[OF delta_gt_SN[OF delta(1)]]]])
let ?main = "roots_of_2_main p ri cr"
have "real_of_2 ` set (?main bnds empty) =
real_of_2 ` set empty ∪
{x. ∃l r. root_cond (p, l, r) x ∧ (l, r) ∈ set bnds} ∧
Ball (set (?main bnds empty)) invariant_2 ∧ distinct (map real_of_2 (?main bnds empty))" (is "?one' ∧ ?two' ∧ ?three'")
using empty bnds disj
proof (induct bnds arbitrary: empty rule: wf_induct[OF wf])
case (1 lrss rais)
note rais = 1(2)[rule_format]
note lrs = 1(3)
note disj = 1(4)
note IH = 1(1)[rule_format]
note simp = roots_of_2_main.simps[of p ri cr lrss rais]
show ?case
proof (cases lrss)
case Nil
with rais show ?thesis unfolding simp by auto
next
case (Cons lr lrs)
obtain l r where lr': "lr = (l,r)" by force
{
fix lr'
assume lt: "⋀ l' r'. (l',r') ∈ set lr' ⟹
l' ≤ r' ∧ delta_gt delta (?r r - ?r l) (?r r' - ?r l')"
have l: "mm (lr' @ lrs) = mm lrs + mm lr'" unfolding mm_def by (auto simp: ac_simps)
have r: "mm lrss = mm lrs + {# ?r r - ?r l #}" unfolding Cons lr' rel_def mm_def
by auto
have "(mm (lr' @ lrs), mm lrss) ∈ mult1 rel'" unfolding l r mult1_def
proof (rule, unfold split, intro exI conjI, unfold add_mset_add_single[symmetric], rule refl, rule refl, intro allI impI)
fix d
assume "d ∈# mm lr'"
then obtain l' r' where d: "d = ?r r' - ?r l'" and lr': "(l',r') ∈ set lr'"
unfolding mm_def in_multiset_in_set by auto
from lt[OF lr']
show "(d, ?r r - ?r l) ∈ rel'" unfolding d rel'_def
by (auto simp: of_rat_less_eq)
qed
hence "(lr' @ lrs, lrss) ∈ rel" unfolding rel_def by auto
} note rel = this
from rel[of Nil] have easy_rel: "(lrs,lrss) ∈ rel" by auto
define c where "c = cr l r"
from simp Cons lr' have simp: "?main lrss rais =
(if c = 0 then ?main lrs rais else if c = 1 then
?main lrs (real_alg_2' ri p l r # rais)
else let m = (l + r) / 2 in ?main ((m, r) # (l, m) # lrs) rais)"
unfolding c_def simp Cons lr' using real_alg_2''[OF False] by auto
note lrs = lrs[unfolded Cons lr']
from lrs have lr: "l ≤ r" by auto
from root_info_condD(1)[OF ri lr, folded cr_def]
have c: "c = card {x. root_cond (p,l,r) x}" unfolding c_def by auto
let ?rt = "λ lrs. {x. ∃l r. root_cond (p, l, r) x ∧ (l, r) ∈ set lrs}"
have rts: "?rt lrss = ?rt lrs ∪ {x. root_cond (p,l,r) x}" (is "?rt1 = ?rt2 ∪ ?rt3")
unfolding Cons lr' by auto
show ?thesis
proof (cases "c = 0")
case True
with simp have simp: "?main lrss rais = ?main lrs rais" by simp
from disj have disj: "pairwise_disjoint (real_of_2 ` set rais # map rts lrs)"
unfolding Cons by auto
from finite_ipoly_roots[OF p] True[unfolded c] have empty: "?rt3 = {}"
unfolding root_cond_def[abs_def] split by simp
with rts have rts: "?rt1 = ?rt2" by auto
show ?thesis unfolding simp rts
by (rule IH[OF easy_rel rais lrs disj], auto)
next
case False
show ?thesis
proof (cases "c = 1")
case True
let ?rai = "real_alg_2' ri p l r"
from True simp have simp: "?main lrss rais = ?main lrs (?rai # rais)" by auto
from card_1_Collect_ex1[OF c[symmetric, unfolded True]]
have ur: "unique_root (p,l,r)" .
from real_alg_2'[OF ur pc ri]
have rai: "invariant_2 ?rai" "real_of_2 ?rai = the_unique_root (p, l, r)" by auto
with rais have rais: "⋀ x. x ∈ set (?rai # rais) ⟹ invariant_2 x"
and dist: "distinct (map real_of_2 rais)" by auto
have rt3: "?rt3 = {real_of_2 ?rai}"
using ur rai by (auto intro: the_unique_root_eqI theI')
have "real_of_2 ` set (roots_of_2_main p ri cr lrs (?rai # rais)) =
real_of_2 ` set (?rai # rais) ∪ ?rt2 ∧
Ball (set (roots_of_2_main p ri cr lrs (?rai # rais))) invariant_2 ∧
distinct (map real_of_2 (roots_of_2_main p ri cr lrs (?rai # rais)))"
(is "?one ∧ ?two ∧ ?three")
proof (rule IH[OF easy_rel, of "?rai # rais", OF conjI lrs])
show "Ball (set (real_alg_2' ri p l r # rais)) invariant_2" using rais by auto
have "real_of_2 (real_alg_2' ri p l r) ∉ set (map real_of_2 rais)"
using disj rt3 unfolding Cons lr' rts_def by auto
thus "distinct (map real_of_2 (real_alg_2' ri p l r # rais))" using dist by auto
show "pairwise_disjoint (real_of_2 ` set (real_alg_2' ri p l r # rais) # map rts lrs)"
using disj rt3 unfolding Cons lr' rts_def by auto
qed auto
hence ?one ?two ?three by blast+
show ?thesis unfolding simp rts rt3
by (rule conjI[OF _ conjI[OF ‹?two› ‹?three›]], unfold ‹?one›, auto)
next
case False
let ?m = "(l+r)/2"
let ?lrs = "[(?m,r),(l,?m)] @ lrs"
from False ‹c ≠ 0› have simp: "?main lrss rais = ?main ?lrs rais"
unfolding simp by (auto simp: Let_def)
from False ‹c ≠ 0› have "c ≥ 2" by auto
from delta(2)[OF this[unfolded c]] have delta: "delta ≤ ?r (r - l) / 4" by auto
have lrs: "⋀ l r. (l,r) ∈ set ?lrs ⟹ l ≤ r"
using lr lrs by (fastforce simp: field_simps)
have "?r ?m ∈ ℚ" unfolding Rats_def by blast
with poly_cond_degree_gt_1[OF pc deg1, of "?r ?m"]
have disj1: "?r ?m ∉ rts lr" for lr unfolding rts_def root_cond_def by auto
have disj2: "rts (?m, r) ∩ rts (l, ?m) = {}" using disj1[of "(l,?m)"] disj1[of "(?m,r)"]
unfolding rts_def root_cond_def by auto
have disj3: "(rts (l,?m) ∪ rts (?m,r)) = rts (l,r)"
unfolding rts_def root_cond_def by (auto simp: hom_distribs)
have disj4: "real_of_2 ` set rais ∩ rts (l,r) = {}" using disj unfolding Cons lr' by auto
have disj: "pairwise_disjoint (real_of_2 ` set rais # map rts ([(?m, r), (l, ?m)] @ lrs))"
using disj disj2 disj3 disj4 by (auto simp: Cons lr')
have "(?lrs,lrss) ∈ rel"
proof (rule rel, intro conjI)
fix l' r'
assume mem: "(l', r') ∈ set [(?m,r),(l,?m)]"
from mem lr show "l' ≤ r'" by auto
from mem have diff: "?r r' - ?r l' = (?r r - ?r l) / 2" by auto
(metis eq_diff_eq minus_diff_eq mult_2_right of_rat_add of_rat_diff,
metis of_rat_add of_rat_mult of_rat_numeral_eq)
show "delta_gt delta (?r r - ?r l) (?r r' - ?r l')" unfolding diff
delta_gt_def by (rule order.trans[OF delta], insert lr,
auto simp: field_simps of_rat_diff of_rat_less_eq)
qed
note IH = IH[OF this, of rais, OF rais lrs disj]
have "real_of_2 ` set (?main ?lrs rais) =
real_of_2 ` set rais ∪ ?rt ?lrs ∧
Ball (set (?main ?lrs rais)) invariant_2 ∧ distinct (map real_of_2 (?main ?lrs rais))"
(is "?one ∧ ?two")
by (rule IH)
hence ?one ?two by blast+
have cong: "⋀ a b c. b = c ⟹ a ∪ b = a ∪ c" by auto
have id: "?rt ?lrs = ?rt lrs ∪ ?rt [(?m,r),(l,?m)]" by auto
show ?thesis unfolding rts simp ‹?one› id
proof (rule conjI[OF cong[OF cong] conjI])
have "⋀ x. root_cond (p,l,r) x = (root_cond (p,l,?m) x ∨ root_cond (p,?m,r) x)"
unfolding root_cond_def by (auto simp:hom_distribs)
hence id: "Collect (root_cond (p,l,r)) = {x. (root_cond (p,l,?m) x ∨ root_cond (p,?m,r) x)}"
by auto
show "?rt [(?m,r),(l,?m)] = Collect (root_cond (p,l,r))" unfolding id list.simps by blast
show "∀ a ∈ set (?main ?lrs rais). invariant_2 a" using ‹?two› by auto
show "distinct (map real_of_2 (?main ?lrs rais))" using ‹?two› by auto
qed
qed
qed
qed
qed
hence idd: "?one'" and cond: ?two' ?three' by blast+
define res where "res = roots_of_2_main p ri cr bnds empty"
have e: "set empty = {}" unfolding empty_def by auto
from idd[folded res_def] e have idd: "real_of_2 ` set res = {} ∪ {x. ∃l r. root_cond (p, l, r) x ∧ (l, r) ∈ set bnds}"
by auto
show ?thesis
unfolding rr unfolding rts id e norm_def using cond
unfolding res_def[symmetric] image_empty e idd[symmetric] by auto
qed
thus ?one ?two ?three by blast+
qed
definition roots_of_2 :: "int poly ⇒ real_alg_2 list" where
"roots_of_2 p = concat (map roots_of_2_irr
(factors_of_int_poly p))"
lemma roots_of_2:
shows "p ≠ 0 ⟹ real_of_2 ` set (roots_of_2 p) = {x. ipoly p x = 0}"
"Ball (set (roots_of_2 p)) invariant_2"
"distinct (map real_of_2 (roots_of_2 p))"
proof -
let ?rr = "roots_of_2 p"
note d = roots_of_2_def
note frp1 = factors_of_int_poly
{
fix q r
assume "q ∈ set ?rr"
then obtain s where
s: "s ∈ set (factors_of_int_poly p)" and
q: "q ∈ set (roots_of_2_irr s)"
unfolding d by auto
from frp1(1)[OF refl s] have "poly_cond s" "degree s > 0" by (auto simp: poly_cond_def)
from roots_of_2_irr[OF this] q
have "invariant_2 q" by auto
}
thus "Ball (set ?rr) invariant_2" by auto
{
assume p: "p ≠ 0"
have "real_of_2 ` set ?rr = (⋃ ((λ p. real_of_2 ` set (roots_of_2_irr p)) `
(set (factors_of_int_poly p))))"
(is "_ = ?rrr")
unfolding d set_concat set_map by auto
also have "… = {x. ipoly p x = 0}"
proof -
{
fix x
assume "x ∈ ?rrr"
then obtain q s where
s: "s ∈ set (factors_of_int_poly p)" and
q: "q ∈ set (roots_of_2_irr s)" and
x: "x = real_of_2 q" by auto
from frp1(1)[OF refl s] have s0: "s ≠ 0" and pt: "poly_cond s" "degree s > 0"
by (auto simp: poly_cond_def)
from roots_of_2_irr[OF pt] q have rt: "ipoly s x = 0" unfolding x by auto
from frp1(2)[OF refl p, of x] rt s have rt: "ipoly p x = 0" by auto
}
moreover
{
fix x :: real
assume rt: "ipoly p x = 0"
from rt frp1(2)[OF refl p, of x] obtain s where s: "s ∈ set (factors_of_int_poly p)"
and rt: "ipoly s x = 0" by auto
from frp1(1)[OF refl s] have s0: "s ≠ 0" and ty: "poly_cond s" "degree s > 0"
by (auto simp: poly_cond_def)
from roots_of_2_irr(1)[OF ty] rt obtain q where
q: "q ∈ set (roots_of_2_irr s)" and
x: "x = real_of_2 q" by blast
have "x ∈ ?rrr" unfolding x using q s by auto
}
ultimately show ?thesis by auto
qed
finally show "real_of_2 ` set ?rr = {x. ipoly p x = 0}" by auto
}
show "distinct (map real_of_2 (roots_of_2 p))"
proof (cases "p = 0")
case True
from factors_of_int_poly_const[of 0] True show ?thesis unfolding roots_of_2_def by auto
next
case p: False
note frp1 = frp1[OF refl]
let ?fp = "factors_of_int_poly p"
let ?cc = "concat (map roots_of_2_irr ?fp)"
show ?thesis unfolding roots_of_2_def distinct_conv_nth length_map
proof (intro allI impI notI)
fix i j
assume ij: "i < length ?cc" "j < length ?cc" "i ≠ j" and id: "map real_of_2 ?cc ! i = map real_of_2 ?cc ! j"
from ij id have id: "real_of_2 (?cc ! i) = real_of_2 (?cc ! j)" by auto
from nth_concat_diff[OF ij, unfolded length_map] obtain j1 k1 j2 k2 where
*: "(j1,k1) ≠ (j2,k2)"
"j1 < length ?fp" "j2 < length ?fp" and
"k1 < length (map roots_of_2_irr ?fp ! j1)"
"k2 < length (map roots_of_2_irr ?fp ! j2)"
"?cc ! i = map roots_of_2_irr ?fp ! j1 ! k1"
"?cc ! j = map roots_of_2_irr ?fp ! j2 ! k2" by blast
hence **: "k1 < length (roots_of_2_irr (?fp ! j1))"
"k2 < length (roots_of_2_irr (?fp ! j2))"
"?cc ! i = roots_of_2_irr (?fp ! j1) ! k1"
"?cc ! j = roots_of_2_irr (?fp ! j2) ! k2"
by auto
from * have mem: "?fp ! j1 ∈ set ?fp" "?fp ! j2 ∈ set ?fp" by auto
from frp1(1)[OF mem(1)] frp1(1)[OF mem(2)]
have pc1: "poly_cond (?fp ! j1)" "degree (?fp ! j1) > 0" and pc10: "?fp ! j1 ≠ 0"
and pc2: "poly_cond (?fp ! j2)" "degree (?fp ! j2) > 0"
by (auto simp: poly_cond_def)
show False
proof (cases "j1 = j2")
case True
with * have neq: "k1 ≠ k2" by auto
from **[unfolded True] id *
have "map real_of_2 (roots_of_2_irr (?fp ! j2)) ! k1 = real_of_2 (?cc ! j)"
"map real_of_2 (roots_of_2_irr (?fp ! j2)) ! k1 = real_of_2 (?cc ! j)"
by auto
hence "¬ distinct (map real_of_2 (roots_of_2_irr (?fp ! j2)))"
unfolding distinct_conv_nth using * ** True by auto
with roots_of_2_irr(3)[OF pc2] show False by auto
next
case neq: False
with frp1(4)[of p] * have neq: "?fp ! j1 ≠ ?fp ! j2" unfolding distinct_conv_nth by auto
let ?x = "real_of_2 (?cc ! i)"
define x where "x = ?x"
from ** have "x ∈ real_of_2 ` set (roots_of_2_irr (?fp ! j1))" unfolding x_def by auto
with roots_of_2_irr(1)[OF pc1] have x1: "ipoly (?fp ! j1) x = 0" by auto
from ** id have "x ∈ real_of_2 ` set (roots_of_2_irr (?fp ! j2))" unfolding x_def
by (metis image_eqI nth_mem)
with roots_of_2_irr(1)[OF pc2] have x2: "ipoly (?fp ! j2) x = 0" by auto
have "ipoly p x = 0" using x1 mem unfolding roots_of_2_def by (metis frp1(2) p)
from frp1(3)[OF p this] x1 x2 neq mem show False by blast
qed
qed
qed
qed
lift_definition (code_dt) roots_of_3 :: "int poly ⇒ real_alg_3 list" is roots_of_2
by (insert roots_of_2, auto simp: list_all_iff)
lemma roots_of_3:
shows "p ≠ 0 ⟹ real_of_3 ` set (roots_of_3 p) = {x. ipoly p x = 0}"
"distinct (map real_of_3 (roots_of_3 p))"
proof -
show "p ≠ 0 ⟹ real_of_3 ` set (roots_of_3 p) = {x. ipoly p x = 0}"
by (transfer; intro roots_of_2, auto)
show "distinct (map real_of_3 (roots_of_3 p))"
by (transfer; insert roots_of_2, auto)
qed
lift_definition roots_of_real_alg :: "int poly ⇒ real_alg list" is roots_of_3 .
lemma roots_of_real_alg:
"p ≠ 0 ⟹ real_of ` set (roots_of_real_alg p) = {x. ipoly p x = 0}"
"distinct (map real_of (roots_of_real_alg p))"
proof -
show "p ≠ 0 ⟹ real_of ` set (roots_of_real_alg p) = {x. ipoly p x = 0}"
by (transfer', insert roots_of_3, auto)
show "distinct (map real_of (roots_of_real_alg p))"
by (transfer, insert roots_of_3(2), auto)
qed
definition real_roots_of_int_poly :: "int poly ⇒ real list" where
"real_roots_of_int_poly p = map real_of (roots_of_real_alg p)"
definition real_roots_of_rat_poly :: "rat poly ⇒ real list" where
"real_roots_of_rat_poly p = map real_of (roots_of_real_alg (snd (rat_to_int_poly p)))"
abbreviation rpoly :: "rat poly ⇒ 'a :: field_char_0 ⇒ 'a"
where "rpoly f ≡ poly (map_poly of_rat f)"
lemma real_roots_of_int_poly: "p ≠ 0 ⟹ set (real_roots_of_int_poly p) = {x. ipoly p x = 0}"
"distinct (real_roots_of_int_poly p)"
unfolding real_roots_of_int_poly_def using roots_of_real_alg[of p] by auto
lemma real_roots_of_rat_poly: "p ≠ 0 ⟹ set (real_roots_of_rat_poly p) = {x. rpoly p x = 0}"
"distinct (real_roots_of_rat_poly p)"
proof -
obtain c q where cq: "rat_to_int_poly p = (c,q)" by force
from rat_to_int_poly[OF this]
have pq: "p = smult (inverse (of_int c)) (of_int_poly q)"
and c: "c ≠ 0" by auto
have id: "{x. rpoly p x = (0 :: real)} = {x. ipoly q x = 0}"
unfolding pq by (simp add: c of_rat_of_int_poly hom_distribs)
show "distinct (real_roots_of_rat_poly p)" unfolding real_roots_of_rat_poly_def cq snd_conv
using roots_of_real_alg(2)[of q] .
assume "p ≠ 0"
with pq c have q: "q ≠ 0" by auto
show "set (real_roots_of_rat_poly p) = {x. rpoly p x = 0}" unfolding id
unfolding real_roots_of_rat_poly_def cq snd_conv using roots_of_real_alg(1)[OF q]
by auto
qed
end