# Theory Power_Sum_Polynomials.Power_Sum_Polynomials_Library

```(*
File:     Power_Sum_Polynomials_Library.thy
Author:   Manuel Eberl, TU München
*)
section ‹Auxiliary material›
theory Power_Sum_Polynomials_Library
imports
"Polynomial_Factorization.Fundamental_Theorem_Algebra_Factorized"
"Symmetric_Polynomials.Symmetric_Polynomials"
"HOL-Computational_Algebra.Computational_Algebra"
begin

unbundle multiset.lifting

subsection ‹Miscellaneous›

lemma atLeastAtMost_nat_numeral:
"atLeastAtMost m (numeral k :: nat) =
(if m ≤ numeral k then insert (numeral k) (atLeastAtMost m (pred_numeral k))
else {})"

lemma sum_in_Rats [intro]: "(⋀x. x ∈ A ⟹ f x ∈ ℚ) ⟹ sum f A ∈ ℚ"
by (induction A rule: infinite_finite_induct) auto

(* TODO Move *)
lemma (in monoid_mult) prod_list_distinct_conv_prod_set:
"distinct xs ⟹ prod_list (map f xs) = prod f (set xs)"
by (induct xs) simp_all

lemma (in monoid_mult) interv_prod_list_conv_prod_set_nat:
"prod_list (map f [m..<n]) = prod f (set [m..<n])"

lemma (in monoid_mult) prod_list_prod_nth:
"prod_list xs = (∏i = 0..< length xs. xs ! i)"
using interv_prod_list_conv_prod_set_nat [of "(!) xs" 0 "length xs"] by (simp add: map_nth)

lemma gcd_poly_code_aux_reduce:
"gcd_poly_code_aux p 0 = normalize p"
"q ≠ 0 ⟹ gcd_poly_code_aux p q = gcd_poly_code_aux q (primitive_part (pseudo_mod p q))"
by (subst gcd_poly_code_aux.simps; simp)+

lemma coprimeI_primes:
fixes a b :: "'a :: factorial_semiring"
assumes "a ≠ 0 ∨ b ≠ 0"
assumes "⋀p. prime p ⟹ p dvd a ⟹ p dvd b ⟹ False"
shows   "coprime a b"
proof (rule coprimeI)
fix d assume d: "d dvd a" "d dvd b"
with assms(1) have [simp]: "d ≠ 0" by auto
show "is_unit d"
proof (rule ccontr)
assume "¬is_unit d"
then obtain p where p: "prime p" "p dvd d"
using prime_divisor_exists[of d] by auto
from assms(2)[of p] and p and d show False
using dvd_trans by auto
qed
qed

lemma coprime_pderiv_imp_squarefree:
assumes "coprime p (pderiv p)"
shows   "squarefree p"
proof (rule squarefreeI)
fix d assume d: "d ^ 2 dvd p"
then obtain q where q: "p = d ^ 2 * q"
by (elim dvdE)
hence "d dvd p" "d dvd pderiv p"
by (auto simp: pderiv_mult pderiv_power_Suc numeral_2_eq_2)
with assms show "is_unit d"
using not_coprimeI by blast
qed

lemma squarefree_field_poly_iff:
fixes p :: "'a :: {field_char_0,euclidean_ring_gcd,semiring_gcd_mult_normalize} poly"
assumes [simp]: "p ≠ 0"
shows "squarefree p ⟷ coprime p (pderiv p)"
proof
assume "squarefree p"
show "coprime p (pderiv p)"
proof (rule coprimeI_primes)
fix d assume d: "d dvd p" "d dvd pderiv p" "prime d"
from d(1) obtain q where q: "p = d * q"
by (elim dvdE)
from d(2) and q have "d dvd q * pderiv d"
with ‹prime d› have "d dvd q ∨ d dvd pderiv d"
using prime_dvd_mult_iff by blast
thus False
proof
assume "d dvd q"
hence "d ^ 2 dvd p"
by (auto simp: q power2_eq_square)
with ‹squarefree p› show False
using d(3) not_prime_unit squarefreeD by blast
next
assume "d dvd pderiv d"
hence "Polynomial.degree d = 0" by simp
moreover have "d ≠ 0" using d by auto
ultimately show False
using d(3) is_unit_iff_degree not_prime_unit by blast
qed
qed auto
qed (use coprime_pderiv_imp_squarefree[of p] in auto)

lemma coprime_pderiv_imp_rsquarefree:
assumes "coprime (p :: 'a :: field_char_0 poly) (pderiv p)"
shows   "rsquarefree p"
unfolding rsquarefree_roots
proof safe
fix x assume "poly p x = 0" "poly (pderiv p) x = 0"
hence "[:-x, 1:] dvd p" "[:-x, 1:] dvd pderiv p"
by (auto simp: poly_eq_0_iff_dvd)
with assms have "is_unit [:-x, 1:]"
using not_coprimeI by blast
thus False by auto
qed

lemma poly_of_nat [simp]: "poly (of_nat n) x = of_nat n"
by (induction n) auto

lemma poly_of_int [simp]: "poly (of_int n) x = of_int n"
by (cases n) auto

lemma order_eq_0_iff: "p ≠ 0 ⟹ order x p = 0 ⟷ poly p x ≠ 0"
by (auto simp: order_root)

lemma order_pos_iff: "p ≠ 0 ⟹ order x p > 0 ⟷ poly p x = 0"
by (auto simp: order_root)

lemma order_prod:
assumes "⋀x. x ∈ A ⟹ f x ≠ 0"
shows   "order x (∏y∈A. f y) = (∑y∈A. order x (f y))"
using assms by (induction A rule: infinite_finite_induct) (auto simp: order_mult)

lemma order_prod_mset:
assumes "0 ∉# A"
shows   "order x (prod_mset A) = sum_mset (image_mset (order x) A)"
using assms by (induction A) (auto simp: order_mult)

lemma order_prod_list:
assumes "0 ∉ set xs"
shows   "order x (prod_list xs) = sum_list (map (order x) xs)"
using assms by (induction xs) (auto simp: order_mult)

lemma order_power: "p ≠ 0 ⟹ order x (p ^ n) = n * order x p"
by (induction n) (auto simp: order_mult)

lemma smult_0_right [simp]: "MPoly_Type.smult p 0 = 0"
by (transfer, transfer) auto

lemma mult_smult_right [simp]:
fixes c :: "'a :: comm_semiring_0"
shows "p * MPoly_Type.smult c q = MPoly_Type.smult c (p * q)"

lemma mapping_single_eq_iff [simp]:
"Poly_Mapping.single a b = Poly_Mapping.single c d ⟷ b = 0 ∧ d = 0 ∨ a = c ∧ b = d"
by transfer (unfold fun_eq_iff when_def, metis)

lemma monom_of_set_plus_monom_of_set:
assumes "A ∩ B = {}" "finite A" "finite B"
shows   "monom_of_set A + monom_of_set B = monom_of_set (A ∪ B)"
using assms by transfer (auto simp: fun_eq_iff)

lemma mpoly_monom_0_eq_Const: "monom 0 c = Const c"
by (intro mpoly_eqI) (auto simp: coeff_monom when_def mpoly_coeff_Const)

lemma mpoly_Const_0 [simp]: "Const 0 = 0"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const mpoly_coeff_0)

lemma mpoly_Const_1 [simp]: "Const 1 = 1"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const mpoly_coeff_1)

lemma mpoly_Const_uminus: "Const (-a) = -Const a"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const)

lemma mpoly_Const_add: "Const (a + b) = Const a + Const b"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const)

lemma mpoly_Const_mult: "Const (a * b) = Const a * Const b"
unfolding mpoly_monom_0_eq_Const [symmetric] mult_monom by simp

lemma mpoly_Const_power: "Const (a ^ n) = Const a ^ n"
by (induction n) (auto simp: mpoly_Const_mult)

lemma of_nat_mpoly_eq: "of_nat n = Const (of_nat n)"
proof (induction n)
case 0
have "0 = (Const 0 :: 'a mpoly)"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const)
thus ?case
by simp
next
case (Suc n)
have "1 + Const (of_nat n) = Const (1 + of_nat n)"
by (intro mpoly_eqI) (auto simp: mpoly_coeff_Const mpoly_coeff_1)
thus ?case
using Suc by auto
qed

lemma insertion_of_nat [simp]: "insertion f (of_nat n) = of_nat n"

lemma insertion_monom_of_set [simp]:
"insertion f (monom (monom_of_set X) c) = c * (∏i∈X. f i)"
proof (cases "finite X")
case [simp]: True
have "insertion f (monom (monom_of_set X) c) = c * (∏a. f a ^ (if a ∈ X then 1 else 0))"
by (auto simp: lookup_monom_of_set)
also have "(∏a. f a ^ (if a ∈ X then 1 else 0)) = (∏i∈X. f i ^ (if i ∈ X then 1 else 0))"
by (intro Prod_any.expand_superset) auto
also have "… = (∏i∈X. f i)"
by (intro prod.cong) auto
finally show ?thesis .
qed (auto simp: lookup_monom_of_set)

(* TODO: Move! Version in AFP is too weak! *)
lemma symmetric_mpoly_symmetric_sum:
assumes "⋀π. π permutes A ⟹ g π permutes X"
assumes "⋀x π. x ∈ X ⟹ π permutes A ⟹ mpoly_map_vars π (f x) = f (g π x)"
shows "symmetric_mpoly A (∑x∈X. f x)"
unfolding symmetric_mpoly_def
proof safe
fix π assume π: "π permutes A"
have "mpoly_map_vars π (sum f X) = (∑x∈X. mpoly_map_vars π (f x))"
by simp
also have "… = (∑x∈X. f (g π x))"
by (intro sum.cong assms π refl)
also have "… = (∑x∈g π`X. f x)"
using assms(1)[OF π] by (subst sum.reindex) (auto simp: permutes_inj_on)
also have "g π ` X = X"
using assms(1)[OF π] by (simp add: permutes_image)
finally show "mpoly_map_vars π (sum f X) = sum f X" .
qed

lemma sym_mpoly_0 [simp]:
assumes "finite A"
shows   "sym_mpoly A 0 = 1"
using assms by (transfer, transfer) (auto simp: fun_eq_iff when_def)

lemma sym_mpoly_eq_0 [simp]:
assumes "k > card A"
shows   "sym_mpoly A k = 0"
proof (transfer fixing: A k, transfer fixing: A k, intro ext)
fix mon
have "¬(finite A ∧ (∃Y⊆A. card Y = k ∧ mon = monom_of_set Y))"
proof safe
fix Y assume Y: "finite A" "Y ⊆ A" "k = card Y" "mon = monom_of_set Y"
hence "card Y ≤ card A" by (intro card_mono) auto
with Y and assms show False by simp
qed
thus "(if finite A ∧ (∃Y⊆A. card Y = k ∧ mon = monom_of_set Y) then 1 else 0) = 0"
by auto
qed

lemma coeff_sym_mpoly_monom_of_set_eq_0:
assumes "finite X" "Y ⊆ X" "card Y ≠ k"
shows   "MPoly_Type.coeff (sym_mpoly X k) (monom_of_set Y) = 0"
using assms finite_subset[of _ X] by (auto simp: coeff_sym_mpoly)

lemma coeff_sym_mpoly_monom_of_set_eq_0':
assumes "finite X" "¬Y ⊆ X" "finite Y"
shows   "MPoly_Type.coeff (sym_mpoly X k) (monom_of_set Y) = 0"
using assms finite_subset[of _ X] by (auto simp: coeff_sym_mpoly)

subsection ‹The set of roots of a univariate polynomial›

lift_definition poly_roots :: "'a :: idom poly ⇒ 'a multiset" is
"λp x. if p = 0 then 0 else order x p"
proof -
fix p :: "'a poly"
show "finite {x. 0 < (if p = 0 then 0 else order x p)}"
by (cases "p = 0") (auto simp: order_pos_iff poly_roots_finite)
qed

lemma poly_roots_0 [simp]: "poly_roots 0 = {#}"
by transfer auto

lemma poly_roots_1 [simp]: "poly_roots 1 = {#}"
by transfer auto

lemma count_poly_roots [simp]:
assumes "p ≠ 0"
shows   "count (poly_roots p) x = order x p"
using assms by transfer auto

lemma in_poly_roots_iff [simp]: "p ≠ 0 ⟹ x ∈# poly_roots p ⟷ poly p x = 0"
by (subst count_greater_zero_iff [symmetric], subst count_poly_roots) (auto simp: order_pos_iff)

lemma set_mset_poly_roots: "p ≠ 0 ⟹ set_mset (poly_roots p) = {x. poly p x = 0}"
using in_poly_roots_iff[of p] by blast

lemma count_poly_roots': "count (poly_roots p) x = (if p = 0 then 0 else order x p)"
by transfer' auto

lemma poly_roots_const [simp]: "poly_roots [:c:] = {#}"
by (intro multiset_eqI) (auto simp: count_poly_roots' order_eq_0_iff)

lemma poly_roots_linear [simp]: "poly_roots [:-x, 1:] = {#x#}"
by (intro multiset_eqI) (auto simp: count_poly_roots' order_eq_0_iff)

lemma poly_roots_monom [simp]: "c ≠ 0 ⟹ poly_roots (Polynomial.monom c n) = replicate_mset n 0"
by (intro multiset_eqI) (auto simp: count_poly_roots' order_eq_0_iff poly_monom)

lemma poly_roots_smult [simp]: "c ≠ 0 ⟹ poly_roots (Polynomial.smult c p) = poly_roots p"
by (intro multiset_eqI) (auto simp: count_poly_roots' order_smult)

lemma poly_roots_mult: "p ≠ 0 ⟹ q ≠ 0 ⟹ poly_roots (p * q) = poly_roots p + poly_roots q"
by (intro multiset_eqI) (auto simp: count_poly_roots' order_mult)

lemma poly_roots_prod:
assumes "⋀x. x ∈ A ⟹ f x ≠ 0"
shows   "poly_roots (prod f A) = (∑x∈A. poly_roots (f x))"
using assms by (induction A rule: infinite_finite_induct) (auto simp: poly_roots_mult)

lemma poly_roots_prod_mset:
assumes "0 ∉# A"
shows   "poly_roots (prod_mset A) = sum_mset (image_mset poly_roots A)"
using assms by (induction A) (auto simp: poly_roots_mult)

lemma poly_roots_prod_list:
assumes "0 ∉ set xs"
shows   "poly_roots (prod_list xs) = sum_list (map poly_roots xs)"
using assms by (induction xs) (auto simp: poly_roots_mult)

lemma poly_roots_power: "p ≠ 0 ⟹ poly_roots (p ^ n) = repeat_mset n (poly_roots p)"
by (induction n) (auto simp: poly_roots_mult)

lemma rsquarefree_poly_roots_eq:
assumes "rsquarefree p"
shows   "poly_roots p = mset_set {x. poly p x = 0}"
proof (rule multiset_eqI)
fix x :: 'a
from assms show "count (poly_roots p) x = count (mset_set {x. poly p x = 0}) x"
by (cases "poly p x = 0") (auto simp: poly_roots_finite order_eq_0_iff rsquarefree_def)
qed

lemma rsquarefree_imp_distinct_roots:
assumes "rsquarefree p" and "mset xs = poly_roots p"
shows   "distinct xs"
proof (cases "p = 0")
case [simp]: False
have *: "mset xs = mset_set {x. poly p x = 0}"
using assms by (simp add: rsquarefree_poly_roots_eq)
hence "set_mset (mset xs) = set_mset (mset_set {x. poly p x = 0})"
by (simp only: )
hence [simp]: "set xs = {x. poly p x = 0}"
from * show ?thesis
by (subst distinct_count_atmost_1) (auto simp: poly_roots_finite)
qed (use assms in auto)

lemma poly_roots_factorization:
fixes p c A
assumes [simp]: "c ≠ 0"
defines "p ≡ Polynomial.smult c (prod_mset (image_mset (λx. [:-x, 1:]) A))"
shows   "poly_roots p = A"
proof -
have "poly_roots p = poly_roots (∏x∈#A. [:-x, 1:])"
by (auto simp: p_def)
also have "… = A"
by (subst poly_roots_prod_mset) (auto simp: image_mset.compositionality o_def)
finally show ?thesis .
qed

lemma fundamental_theorem_algebra_factorized':
fixes p :: "complex poly"
shows "p = Polynomial.smult (Polynomial.lead_coeff p)
(prod_mset (image_mset (λx. [:-x, 1:]) (poly_roots p)))"
proof (cases "p = 0")
case [simp]: False
obtain xs where
xs: "Polynomial.smult (Polynomial.lead_coeff p) (∏x←xs. [:-x, 1:]) = p"
"length xs = Polynomial.degree p"
using fundamental_theorem_algebra_factorized[of p] by auto
define A where "A = mset xs"

note xs(1)
also have "(∏x←xs. [:-x, 1:]) = prod_mset (image_mset (λx. [:-x, 1:]) A)"
unfolding A_def by (induction xs) auto
finally have *: "Polynomial.smult (Polynomial.lead_coeff p) (∏x∈#A. [:- x, 1:]) = p" .
also have "A = poly_roots p"
by (subst * [symmetric]) auto
finally show ?thesis ..
qed auto

lemma poly_roots_eq_imp_eq:
fixes p q :: "complex poly"
assumes "poly_roots p = poly_roots q"
shows   "p = q"
proof (cases "p = 0 ∨ q = 0")
case False
hence [simp]: "p ≠ 0" "q ≠ 0"
by auto
have "p = Polynomial.smult (Polynomial.lead_coeff p)
(prod_mset (image_mset (λx. [:-x, 1:]) (poly_roots p)))"
by (rule fundamental_theorem_algebra_factorized')
also have "… = Polynomial.smult (Polynomial.lead_coeff q)
(prod_mset (image_mset (λx. [:-x, 1:]) (poly_roots q)))"
also have "… = q"
by (rule fundamental_theorem_algebra_factorized' [symmetric])
finally show ?thesis .
qed (use assms in auto)

lemma Sum_any_zeroI': "(⋀x. P x ⟹ f x = 0) ⟹ Sum_any (λx. f x when P x) = 0"
by (auto simp: Sum_any.expand_set)

(* TODO: This was not really needed here, but it is important nonetheless.
It should go in the Symmetric_Polynomials entry. *)
lemma sym_mpoly_insert:
assumes "finite X" "x ∉ X"
shows   "(sym_mpoly (insert x X) (Suc k) :: 'a :: semiring_1 mpoly) =
monom (monom_of_set {x}) 1 * sym_mpoly X k + sym_mpoly X (Suc k)" (is "?lhs = ?A + ?B")
proof (rule mpoly_eqI)
fix mon
show "coeff ?lhs mon = coeff (?A + ?B) mon"
proof (cases "∀i. lookup mon i ≤ 1 ∧ (i ∉ insert x X ⟶ lookup mon i = 0)")
case False
then obtain i where i: "lookup mon i > 1 ∨ i ∉ insert x X ∧ lookup mon i > 0"
by (auto simp: not_le)

have "coeff ?A mon = prod_fun (coeff (monom (monom_of_set {x}) 1))
(coeff (sym_mpoly X k)) mon"
also have "… = (∑l. ∑q. coeff (monom (monom_of_set {x}) 1) l * coeff (sym_mpoly X k) q
when mon = l + q)"
unfolding prod_fun_def
by (intro Sum_any.cong, subst Sum_any_right_distrib, force)
(auto simp: Sum_any_right_distrib when_def intro!: Sum_any.cong)
also have "… = 0"
proof (rule Sum_any_zeroI, rule Sum_any_zeroI')
fix ma mb assume *: "mon = ma + mb"
show "coeff (monom (monom_of_set {x}) (1::'a)) ma * coeff (sym_mpoly X k) mb = 0"
proof (cases "i = x")
case [simp]: True
show ?thesis
proof (cases "lookup mb i > 0")
case True
hence "coeff (sym_mpoly X k) mb = 0" using ‹x ∉ X›
by (auto simp: coeff_sym_mpoly lookup_monom_of_set split: if_splits)
thus ?thesis
using mult_not_zero by blast
next
case False
hence "coeff (monom (monom_of_set {x}) 1) ma = 0"
using i by (auto simp: coeff_monom when_def * lookup_add)
thus ?thesis
using mult_not_zero by blast
qed
next
case [simp]: False
show ?thesis
proof (cases "lookup ma i > 0")
case False
hence "lookup mb i = lookup mon i"
using * by (auto simp: lookup_add)
hence "coeff (sym_mpoly X k) mb = 0" using i
by (auto simp: coeff_sym_mpoly lookup_monom_of_set split: if_splits)
thus ?thesis
using mult_not_zero by blast
next
case True
hence "coeff (monom (monom_of_set {x}) 1) ma = 0"
using i by (auto simp: coeff_monom when_def * lookup_add)
thus ?thesis
using mult_not_zero by blast
qed
qed
qed
finally have "coeff ?A mon = 0" .
moreover from False have "coeff ?lhs mon = 0"
by (subst coeff_sym_mpoly) (auto simp: lookup_monom_of_set split: if_splits)
moreover from False have "coeff (sym_mpoly X (Suc k)) mon = 0"
by (subst coeff_sym_mpoly) (auto simp: lookup_monom_of_set split: if_splits)
ultimately show ?thesis
by auto
next
case True
define A where "A = keys mon"
have A: "A ⊆ insert x X"
using True by (auto simp: A_def)
have [simp]: "mon = monom_of_set A"
unfolding A_def using True by transfer (force simp: fun_eq_iff le_Suc_eq)
have "finite A"
using finite_subset A assms by blast
show ?thesis
proof (cases "x ∈ A")
case False
have "coeff ?A mon = prod_fun (coeff (monom (monom_of_set {x}) 1))
(coeff (sym_mpoly X k)) (monom_of_set A)"
also have "… = (∑l. ∑q. coeff (monom (monom_of_set {x}) 1) l * coeff (sym_mpoly X k) q
when monom_of_set A = l + q)"
unfolding prod_fun_def
by (intro Sum_any.cong, subst Sum_any_right_distrib, force)
(auto simp: Sum_any_right_distrib when_def intro!: Sum_any.cong)
also have "… = 0"
proof (rule Sum_any_zeroI, rule Sum_any_zeroI')
fix ma mb assume *: "monom_of_set A = ma + mb"
hence "keys ma ⊆ A"
using ‹finite A› by transfer (auto simp: fun_eq_iff split: if_splits)
thus "coeff (monom (monom_of_set {x}) (1::'a)) ma * coeff (sym_mpoly X k) mb = 0"
using ‹x ∉ A› by (auto simp: coeff_monom when_def)
qed
finally show ?thesis
using False A assms finite_subset[of _ "insert x X"] finite_subset[of _ X]
by (auto simp: coeff_sym_mpoly)
next
case True
have "mon = monom_of_set {x} + monom_of_set (A - {x})"
using ‹x ∈ A› ‹finite A› by (auto simp: monom_of_set_plus_monom_of_set)
also have "coeff ?A … = coeff (sym_mpoly X k) (monom_of_set (A - {x}))"
by (subst coeff_monom_mult) auto
also have "… = (if card A = Suc k then 1 else 0)"
proof (cases "card A = Suc k")
case True
thus ?thesis
using assms ‹finite A› ‹x ∈ A› A
by (subst coeff_sym_mpoly_monom_of_set) auto
next
case False
thus ?thesis
using assms ‹x ∈ A› A ‹finite A› card_Suc_Diff1[of A x]
by (subst coeff_sym_mpoly_monom_of_set_eq_0) auto
qed
moreover have "coeff ?B (monom_of_set A) = 0"
using assms ‹x ∈ A› ‹finite A›
by (subst coeff_sym_mpoly_monom_of_set_eq_0') auto
moreover have "coeff ?lhs (monom_of_set A) = (if card A = Suc k then 1 else 0)"
using assms A ‹finite A› finite_subset[of _ "insert x X"] by (auto simp: coeff_sym_mpoly)
ultimately show ?thesis by simp
qed
qed
qed

lifting_update multiset.lifting
lifting_forget multiset.lifting

end```