Theory Real_Roots

(*  
    Author:      René Thiemann 
                 Akihisa Yamada
    License:     BSD
*)
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