Theory Factor_Real_Poly

(* Factorization of polynomials with real algebraic coefficients *)

subsection ‹Real Algebraic Coefficients›

text ‹We basically perform a factorization via complex algebraic numbers,
  take all real roots, and 
  then merge each pair of conjugate roots into a quadratic factor.›

theory Factor_Real_Poly
  imports
    Factor_Complex_Poly
begin

hide_const (open) Coset.order

fun delete_cnj :: "complex  nat  (complex × nat) list  (complex × nat) list" where
  "delete_cnj x i ((y,j) # yjs) = (if x = y then if j = i then yjs else if j > i then
    ((y,j - i) # yjs) else delete_cnj x (i - j) yjs else (y,j) # delete_cnj x i yjs)"
| "delete_cnj _ _ [] = []"

lemma delete_cnj_length[termination_simp]: "length (delete_cnj x i yjs)  length yjs"
  by (induct x i yjs rule: delete_cnj.induct, auto)

fun complex_roots_to_real_factorization :: "(complex × nat) list  (real poly × nat)list" where
  "complex_roots_to_real_factorization [] = []"
| "complex_roots_to_real_factorization ((x,i) # xs) = (if x   then 
    ([:-(Re x),1:],i) # complex_roots_to_real_factorization xs else 
    let xx = cnj x; ys = delete_cnj xx i xs; p = map_poly Re ([:-x,1:] * [:-xx,1:])
    in (p,i) # complex_roots_to_real_factorization ys)"

definition factor_real_poly :: "real poly  real × (real poly × nat) list" where
  "factor_real_poly p  case factor_complex_main (map_poly of_real p) of 
    (c,ris)  (Re c, complex_roots_to_real_factorization ris) "


lemma monic_imp_nonzero: "monic x  x  0" for x :: "'a :: semiring_1 poly" by auto

lemma delete_cnj_0: assumes "0  snd ` set xis" 
  shows "0  snd ` set (delete_cnj x si xis)" 
  using assms by (induct x si xis rule: delete_cnj.induct, auto)

lemma delete_cnj: assumes 
  "order x ((x, i)xis. [:- x, 1:] ^ i)  si" "si  0"
  shows "((x, i)xis. [:- x, 1:] ^ i) =
    [:- x, 1:] ^ si * ((x, i)delete_cnj x si xis. [:- x, 1:] ^ i)"
using assms
proof (induct x si xis rule: delete_cnj.induct)
  case (2 x si)
  hence "order x 1  1" by auto
  hence "[:-x,1:]^1 dvd 1" unfolding order_divides by simp
  from power_le_dvd[OF this, of 1] si  0 have "[:- x, 1:] dvd 1" by simp
  from divides_degree[OF this] 
  show ?case by auto
next
  case (1 x i y j yjs)
  note IH = 1(1-2)
  let ?yj = "[:-y,1:]^j"
  let ?yjs = "((x,i)yjs. [:- x, 1:] ^ i)"
  let ?x = "[: - x, 1 :]"
  let ?xi = "?x ^ i"
  have "monic ((x,i)(y, j) # yjs. [:- x, 1:] ^ i)"
    by (intro monic_prod_list, auto intro: monic_power)
  then have "monic (?yj * ?yjs)" by simp
  from monic_imp_nonzero[OF this] have yy0: "?yj * ?yjs  0" by auto
  have id: "((x,i)(y, j) # yjs. [:- x, 1:] ^ i) = ?yj * ?yjs" by simp
  from 1(3-) have ord: "i  order x (?yj * ?yjs)" and i: "i  0" unfolding id by auto
  from ord[unfolded order_mult[OF yy0]] have ord: "i  order x ?yj + order x ?yjs" .
  from this[unfolded order_linear_power] 
  have ord: "i  (if y = x then j else 0) + order x ?yjs" by simp
  show ?case 
  proof (cases "x = y")
    case False
    from ord False have "i  order x ?yjs" by simp
    note IH = IH(2)[OF False this i]
    from False have del: "delete_cnj x i ((y, j) # yjs) = (y,j) # delete_cnj x i yjs" by simp    
    show ?thesis unfolding del id IH
      by (simp add: ac_simps)
  next
    case True note xy = this
    note IH = IH(1)[OF True]
    show ?thesis 
    proof (cases "j  i")
      case False
      from ord have ord: "i - j  order x ?yjs" unfolding xy by simp
      have "?xi = ?x ^ (j + (i - j))" using False by simp
      also have " = ?x ^ j * ?x ^ (i - j)"
        unfolding power_add by simp
      finally have xi: "?xi = ?x ^ j * ?x ^ (i - j)" .
      from False have "j  i" "¬ i < j" "i - j  0" by auto
      note IH = IH[OF this(1,2) ord this(3)]
      from xy False have del: "delete_cnj x i ((y, j) # yjs) = delete_cnj x (i - j) yjs" by auto
      show ?thesis unfolding del id unfolding IH xi unfolding xy by simp
    next
      case True
      hence "j = i  i < j" by auto
      thus ?thesis
      proof
        assume i: "j = i"
        from xy i have del: "delete_cnj x i ((y, j) # yjs) = yjs" by simp
        show ?thesis unfolding id del unfolding xy i by simp
      next
        assume ij: "i < j"
        with xy i have del: "delete_cnj x i ((y, j) # yjs) = (y, j - i) # yjs" by simp
        from ij have idd: "j = i + (j - i)" by simp
        show ?thesis 
          apply (unfold id del)
          apply (subst idd)  
          apply (unfold power_add xy)
          by simp
      qed
    qed
  qed
qed


theorem factor_real_poly: assumes fp: "factor_real_poly p = (c,qis)"
  shows "p = smult c ((q, i)qis. q ^ i)" 
    "(q,j)  set qis  irreducible q  j  0  monic q  degree q  {1,2}"
proof -
  interpret map_poly_inj_idom_hom of_real..
  have "(p = smult c ((q, i)qis. q ^ i))  ((q,j)  set qis  irreducible q  j  0  monic q  degree q  {1,2})"
  proof (cases "p = 0")
    case True
    have yun: "yun_rf (yun_polys (0 :: complex poly)) = (0,[])" 
      by (transfer, auto simp: yun_factorization_def)
    have "factor_real_poly p = (0,[])" unfolding True
      by (simp add: factor_real_poly_def factor_complex_main_def yun)
    with fp have id: "c = 0" "qis = []" by auto
    thus ?thesis unfolding True by simp
  next
    case False note p0 = this
    let ?c = complex_of_real
    let ?rp = "map_poly Re"
    let ?cp = "map_poly ?c"
    let ?p = "?cp p"
    from fp[unfolded factor_real_poly_def]
      obtain d xis where fp: "factor_complex_main ?p = (d,xis)"
      and c: "c = Re d" and qis: "qis = complex_roots_to_real_factorization xis" 
        by (cases "factor_complex_main ?p", auto)
      from factor_complex_main[OF fp] have p: "?p = smult d ((x, i)xis. [:- x, 1:] ^ i)" 
      (is "_ = smult d ?q") and 0: "0  snd ` set xis" .
    from arg_cong[OF this(1), of "λ p. coeff p (degree p)"]
    have "coeff ?p (degree ?p) = coeff (smult d ?q) (degree (smult d ?q))" .
    also have "coeff ?p (degree ?p) = ?c (coeff p (degree p))" by simp
    also have "coeff (smult d ?q) (degree (smult d ?q)) = d * coeff ?q (degree ?q)"
      by simp
    also have "monic ?q" by (rule monic_prod_list, auto intro: monic_power)
    finally have d: "d = ?c (coeff p (degree p))" by auto
    from arg_cong[OF this, of Re, folded c] have c: "c = coeff p (degree p)" by auto
    have "set (coeffs ?p)  " by auto
    with p have q': "set (coeffs (smult d ?q))  " by auto
    from d p0 have d0: "d  0" by auto
    have "smult d ?q = [:d:] * ?q" by auto
    from real_poly_factor[OF q'[unfolded this]] d0 d
    have q: "set (coeffs ?q)  " by auto 
    have "p = ?rp ?p"
      by (rule sym, subst map_poly_map_poly, force, rule map_poly_idI, auto)
    also have " = ?rp (smult d ?q)" unfolding p ..
    also have "?q = ?cp (?rp ?q)"
      by (rule sym, rule map_poly_of_real_Re, insert q, auto)
    also have "d = ?c c" unfolding d c ..
    also have "smult (?c c) (?cp (?rp ?q)) = ?cp (smult c (?rp ?q))" by (simp add: hom_distribs)
    also have "?rp  = smult c (?rp ?q)" 
      by (subst map_poly_map_poly, force, rule map_poly_idI, auto)
    finally have p: "p = smult c (?rp ?q)" .
    let ?fact = complex_roots_to_real_factorization
    have "?rp ?q = ((q, i)qis. q ^ i) 
      ((q, j)  set qis  irreducible q  j  0  monic q  degree q  {1, 2})"
      using q 0 unfolding qis
    proof (induct xis rule: complex_roots_to_real_factorization.induct)
      case 1
      show ?case by simp
    next
      case (2 x i xis)
      note IH = 2(1-2)
      note prems = 2(3)
      from 2(4) have i: "i  0" and 0: "0  snd ` set xis" by auto
      let ?xi = "[:- x, 1:] ^ i"
      let ?xis = "((x, i)xis. [:- x, 1:] ^ i)"
      have id: "((x, i)((x,i) # xis). [:- x, 1:] ^ i) = ?xi * ?xis"
        by simp
      show ?case
      proof (cases "x  ")
        case True
        have xi: "set (coeffs ?xi)  "
          by (rule real_poly_power, insert True, auto)
        have xis: "set (coeffs ?xis)  " 
          by (rule real_poly_factor[OF prems[unfolded id] xi], rule linear_power_nonzero)
        note IH = IH(1)[OF True xis 0]
        have "?rp (?xi * ?xis) = ?rp ?xi * ?rp ?xis"
          by (rule map_poly_Re_mult[OF xi xis])
        also have "?rp ?xi = (?rp [: -x,1 :])^ i"
          by (rule map_poly_Re_power, insert True, auto)
        also have "?rp [: -x,1 :] = [:-(Re x),1:]" by auto
        also have "?rp ?xis = ( (a,b)  ?fact xis. a ^ b)"
          using IH by auto        
        also have "[:- Re x, 1:] ^ i * ( (a,b)  ?fact xis. a ^ b) = 
          ( (a,b)  ?fact ((x,i) # xis). a ^ b)" using True by simp
        finally have idd: "?rp (?xi * ?xis) = ( (a,b)  ?fact ((x,i) # xis). a ^ b)" .
        show ?thesis unfolding id idd
        proof (intro conjI, force, intro impI)
          assume "(q, j)  set (?fact ((x, i) # xis))"
          hence "(q,j)  set (?fact xis)  (q = [:- Re x, 1:]  j = i)"
            using True by auto
          thus "irreducible q  j  0  monic q  degree q  {1, 2}"
          proof
            assume "(q,j)  set (?fact xis)"
            with IH show ?thesis by auto
          next
            assume "q = [:- Re x, 1:]  j = i"
            with linear_irreducible_field[of "[:- Re x, 1:]"] i show ?thesis by auto
          qed
        qed
      next
        case False
        define xi where "xi = [:Re x * Re x + Im x * Im x, - (2 * Re x), 1:]"
        obtain xx where xx: "xx = cnj x" by auto
        have xi: "xi = ?rp ([:-x,1:] * [:-xx,1:])" unfolding xx xi_def by auto
        have cpxi: "?cp xi = [:-x,1:] * [:-xx,1:]" unfolding xi_def
          by (cases x, auto simp: xx legacy_Complex_simps)
        obtain yis where yis: "yis = delete_cnj xx i xis" by auto
        from delete_cnj_0[OF 0] have 0: "0  snd ` set yis" unfolding yis .
        from False have fact: "?fact ((x,i) # xis) = ((xi,i) # ?fact yis)"
          unfolding xi_def xx yis by simp
        note IH = IH(2)[OF False xx yis xi _ 0]
        have "irreducible xi"
          apply (fold irreducible_connect_field)
        proof (rule irreducibledI)
          show "degree xi > 0" unfolding xi by auto
          fix q p :: "real poly" 
          assume "degree q > 0" "degree q < degree xi" and qp: "xi = q * p"
          hence dq: "degree q = 1" unfolding xi by auto
          have dxi: "degree xi = 2" "xi  0" unfolding xi by auto
          with qp have "q  0" "p  0" by auto
          hence "degree xi = degree q + degree p" unfolding qp
            by (rule degree_mult_eq)
          with dq have dp: "degree p = 1" unfolding dxi by simp
          {
            fix c :: complex
            assume rt: "poly (?cp xi) c = 0"
            hence "poly (?cp q * ?cp p) c = 0" by (simp add: qp hom_distribs)
            hence "(poly (?cp q) c = 0  poly (?cp p) c = 0)" by auto
            hence "c = roots1 (?cp q)  c = roots1 (?cp p)"
              using roots1[of "?cp q"] roots1[of "?cp p"] dp dq by auto
            hence "c  " unfolding roots1_def by auto
            hence "c  x" using False by auto
          }
          hence "poly (?cp xi) x  0" by auto
          thus False unfolding cpxi by simp
        qed
        hence xi': "irreducible xi" "monic xi" "degree xi = 2"
          unfolding xi by auto
        let ?xxi = "[:- xx, 1:] ^ i"
        let ?yis = "((x, i)yis. [:- x, 1:] ^ i)"
        let ?yi = "(?cp xi)^ i"
        have yi: "set (coeffs ?yi)  "
          by (rule real_poly_power, auto simp: xi)
        have mon: "monic ((x, i)(x, i) # xis. [:- x, 1:] ^ i)"
          by (rule monic_prod_list, auto intro: monic_power)
        from monic_imp_nonzero[OF this] have xixis: "?xi * ?xis  0" unfolding id by auto
        from False have xxx: "xx  x" unfolding xx by (cases x, auto simp: legacy_Complex_simps Reals_def)
        from prems[unfolded id] have prems: "set (coeffs (?xi * ?xis))  " .
        from id have "[:- x, 1:] ^ i dvd ?xi * ?xis" by auto
        from xixis this[unfolded order_divides] 
        have "order x (?xi * ?xis)  i" by auto
        with complex_conjugate_order[OF prems xixis, of x, folded xx]
        have "order xx (?xi * ?xis)  i" by auto
        hence "order xx ?xi + order xx ?xis  i" unfolding order_mult[OF xixis] .
        also have "order xx ?xi = 0" unfolding order_linear_power using xxx by simp
        finally have "order xx ?xis  i" by simp
        hence yis: "?xis = ?xxi * ?yis" unfolding yis using i
          by (intro delete_cnj, simp)
        hence "?xi * ?xis = (?xi * ?xxi) * ?yis" by (simp only: ac_simps)
        also have "?xi * ?xxi = ([:- x, 1:] * [:- xx, 1:])^i"
          by (metis power_mult_distrib)
        also have "[:- x, 1:] * [:- xx, 1:] = ?cp xi" unfolding cpxi ..
        finally have idd: "?xi * ?xis = (?cp xi)^i * ?yis" by simp
        from prems[unfolded idd] have R: "set (coeffs ((?cp xi)^i * ?yis))  " .
        have yis: "set (coeffs ?yis)  "
          by (rule real_poly_factor[OF R yi], auto, auto simp: xi_def)
        note IH = IH[OF yis] 
        have "?rp (?xi * ?xis) = ?rp ?yi * ?rp ?yis" unfolding idd
          by (rule map_poly_Re_mult[OF yi yis])
        also have "?rp ?yi = xi^i" by (fold hom_distribs, rule map_poly_Re_of_real)
        also have "?rp ?yis = ( (a,b)  ?fact yis. a ^ b)"
          using IH by auto
        also have "xi ^ i * ( (a,b)  ?fact yis. a ^ b) = 
          ( (a,b)  ?fact ((x,i) # xis). a ^ b)" unfolding fact by simp
        finally have idd: "?rp (?xi * ?xis) = ( (a,b)  ?fact ((x,i) # xis). a ^ b)" .
        show ?thesis unfolding id idd fact using IH xi' i by auto
      qed
    qed
    thus ?thesis unfolding p by simp
  qed
  thus "p = smult c ((q, i)qis. q ^ i)" 
    "(q,j)  set qis  irreducible q  j  0  monic q  degree q  {1,2}" by blast+
qed

end