Theory Power_Sum_Polynomials
section ‹Power sum polynomials›
theory Power_Sum_Polynomials
imports
  "Symmetric_Polynomials.Symmetric_Polynomials"
  "HOL-Computational_Algebra.Field_as_Ring"
  Power_Sum_Polynomials_Library
begin
subsection ‹Definition›
text ‹
  For $n$ indeterminates $X_1,\ldots,X_n$, we define the $k$-th power sum polynomial as
  \[p_k(X_1, \ldots, X_n) = X_1^k + \ldots + X_n^k\ .\]
›
lift_definition powsum_mpoly_aux :: "nat set ⇒ nat ⇒ (nat ⇒⇩0 nat) ⇒⇩0 'a :: {semiring_1,zero_neq_one}" is
  "λX k mon. if infinite X ∨ k = 0 ∧ mon ≠ 0 then 0
             else if k = 0 ∧ mon = 0 then of_nat (card X)
             else if finite X ∧ (∃x∈X. mon = Poly_Mapping.single x k) then 1 else 0"
  by auto
lemma lookup_powsum_mpoly_aux:
  "Poly_Mapping.lookup (powsum_mpoly_aux X k) mon =
     (if infinite X ∨ k = 0 ∧ mon ≠ 0 then 0
             else if k = 0 ∧ mon = 0 then of_nat (card X)
             else if finite X ∧ (∃x∈X. mon = Poly_Mapping.single x k) then 1 else 0)"
  by transfer' simp
lemma lookup_sym_mpoly_aux_monom_singleton [simp]:
  assumes "finite X" "x ∈ X" "k > 0"
  shows   "Poly_Mapping.lookup (powsum_mpoly_aux X k) (Poly_Mapping.single x k) = 1"
  using assms by (auto simp: lookup_powsum_mpoly_aux)
lemma lookup_sym_mpoly_aux_monom_singleton':
  assumes "finite X" "k > 0"
  shows   "Poly_Mapping.lookup (powsum_mpoly_aux X k) (Poly_Mapping.single x k) = (if x ∈ X then 1 else 0)"
  using assms by (auto simp: lookup_powsum_mpoly_aux)
lemma keys_powsum_mpoly_aux: "m ∈ keys (powsum_mpoly_aux A k) ⟹ keys m ⊆ A"
  by transfer' (auto split: if_splits simp: keys_monom_of_set)
lift_definition powsum_mpoly :: "nat set ⇒ nat ⇒ 'a :: {semiring_1,zero_neq_one} mpoly" is
  "powsum_mpoly_aux" .
lemma vars_powsum_mpoly_subset: "vars (powsum_mpoly A k) ⊆ A"
  using keys_powsum_mpoly_aux by (auto simp: vars_def powsum_mpoly.rep_eq)
lemma powsum_mpoly_infinite: "¬finite A ⟹ powsum_mpoly A k = 0"
  by (transfer, transfer) auto
lemma coeff_powsum_mpoly:
  "MPoly_Type.coeff (powsum_mpoly X k) mon =
     (if infinite X ∨ k = 0 ∧ mon ≠ 0 then 0
             else if k = 0 ∧ mon = 0 then of_nat (card X)
             else if finite X ∧ (∃x∈X. mon = Poly_Mapping.single x k) then 1 else 0)"
  by transfer' (simp add: lookup_powsum_mpoly_aux)
lemma coeff_powsum_mpoly_0_right:
  "MPoly_Type.coeff (powsum_mpoly X 0) mon = (if mon = 0 then of_nat (card X) else 0)"
  by transfer' (auto simp add: lookup_powsum_mpoly_aux)
lemma coeff_powsum_mpoly_singleton:
  assumes "finite X" "k > 0"
  shows   "MPoly_Type.coeff (powsum_mpoly X k) (Poly_Mapping.single x k) = (if x ∈ X then 1 else 0)"
  using assms by transfer' (simp add: lookup_powsum_mpoly_aux)
lemma coeff_powsum_mpoly_singleton_eq_1 [simp]:
  assumes "finite X" "x ∈ X" "k > 0"
  shows   "MPoly_Type.coeff (powsum_mpoly X k) (Poly_Mapping.single x k) = 1"
  using assms by (simp add: coeff_powsum_mpoly_singleton)
lemma coeff_powsum_mpoly_singleton_eq_0 [simp]:
  assumes "finite X" "x ∉ X" "k > 0"
  shows   "MPoly_Type.coeff (powsum_mpoly X k) (Poly_Mapping.single x k) = 0"
  using assms by (simp add: coeff_powsum_mpoly_singleton)
lemma powsum_mpoly_0 [simp]: "powsum_mpoly X 0 = of_nat (card X)"
  by (intro mpoly_eqI ext) (auto simp: coeff_powsum_mpoly_0_right of_nat_mpoly_eq mpoly_coeff_Const)
lemma powsum_mpoly_empty [simp]: "powsum_mpoly {} k = 0"
  by (intro mpoly_eqI) (auto simp: coeff_powsum_mpoly)
lemma powsum_mpoly_altdef: "powsum_mpoly X k = (∑x∈X. monom (Poly_Mapping.single x k) 1)"
proof (cases "finite X")
  case [simp]: True
  show ?thesis
  proof (cases "k = 0")
    case True
    thus ?thesis by auto
  next
    case False
    show ?thesis
    proof (intro mpoly_eqI, goal_cases)
      case (1 mon)
      show ?case using False
        by (cases "∃x∈X. mon = Poly_Mapping.single x k")
           (auto simp: coeff_powsum_mpoly coeff_monom when_def)
    qed
  qed
qed (auto simp: powsum_mpoly_infinite)
text ‹
  Power sum polynomials are symmetric:
›
lemma symmetric_powsum_mpoly [intro]:
  assumes "A ⊆ B"
  shows   "symmetric_mpoly A (powsum_mpoly B k)"
  unfolding powsum_mpoly_altdef
proof (rule symmetric_mpoly_symmetric_sum)
  fix x π
  assume "x ∈ B" "π permutes A"
  thus "mpoly_map_vars π (MPoly_Type.monom (Poly_Mapping.single x k) 1) =
        MPoly_Type.monom (Poly_Mapping.single (π x) k) 1"
    using assms by (auto simp: mpoly_map_vars_monom permutes_bij permutep_single
                               bij_imp_bij_inv permutes_inv_inv)
qed (use assms in ‹auto simp: permutes_subset›)
lemma insertion_powsum_mpoly [simp]: "insertion f (powsum_mpoly X k) = (∑i∈X. f i ^ k)"
  unfolding powsum_mpoly_altdef insertion_sum insertion_single by simp
lemma powsum_mpoly_nz:
  assumes "finite X" "X ≠ {}" "k > 0"
  shows   "(powsum_mpoly X k :: 'a :: {semiring_1, zero_neq_one} mpoly) ≠ 0"
proof -
  from assms obtain x where "x ∈ X" by auto
  hence "coeff (powsum_mpoly X k) (Poly_Mapping.single x k) = (1 :: 'a)"
    using assms by (auto simp: coeff_powsum_mpoly)
  thus ?thesis by auto
qed
lemma powsum_mpoly_eq_0_iff:
  assumes "k > 0"
  shows   "powsum_mpoly X k = 0 ⟷ infinite X ∨ X = {}"
  using assms powsum_mpoly_nz[of X k] by (auto simp: powsum_mpoly_infinite)
subsection ‹The Girard--Newton Theorem›
text ‹
  The following is a nice combinatorial proof of the Girard--Newton Theorem due to
  Doron Zeilberger~\<^cite>‹"zeilberger"›.
  The precise statement is this:
  Let $e_k$ denote the $k$-th elementary symmetric polynomial in $X_1,\ldots,X_n$.
  This is the sum of all monomials that can be formed by taking the product of $k$ 
  distinct variables.
  Next, let $p_k = X_1^k + \ldots + X_n^k$ denote that $k$-th symmetric power sum polynomial
  in $X_1,\ldots,X_n$.
  Then the following equality holds:
  \[(-1)^k k e_k + \sum_{i=0}^{k-1} (-1)^i e_i p_{k-i}\]
›
theorem Girard_Newton:
  assumes "finite X"
  shows   "(-1) ^ k * of_nat k * sym_mpoly X k +
           (∑i<k. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i)) =
             (0 :: 'a :: comm_ring_1 mpoly)"
  (is "?lhs = 0")
proof -
  write Poly_Mapping.single (‹sng›)
  define n where "n = card X"
  define 𝒜 :: "(nat set × nat) set"
    where "𝒜 = {(A, j). A ⊆ X ∧ card A ≤ k ∧ j ∈ X ∧ (card A = k ⟶ j ∈ A)}"
  define 𝒜1 :: "(nat set × nat) set"
    where "𝒜1 = {A∈Pow X. card A < k} × X"
  define 𝒜2 :: "(nat set × nat) set"
    where "𝒜2 = (SIGMA A:{A∈Pow X. card A = k}. A)"
  have 𝒜_split: "𝒜 = 𝒜1 ∪ 𝒜2" "𝒜1 ∩ 𝒜2 = {}"
    by (auto simp: 𝒜_def 𝒜1_def 𝒜2_def)
  have [intro]: "finite 𝒜1" "finite 𝒜2"
    using assms finite_subset[of _ X] by (auto simp: 𝒜1_def 𝒜2_def intro!: finite_SigmaI)
  have [intro]: "finite 𝒜"
    by (subst 𝒜_split) auto
  
  define w :: "nat set × nat ⇒ 'a mpoly"
    where "w = (λ(A, j). monom (monom_of_set A + sng j (k - card A)) ((-1) ^ card A))"
  
  have "?lhs = (∑x∈𝒜. w x)"
  proof -
    have "(∑x∈𝒜. w x) = (∑x∈𝒜1. w x) + (∑x∈𝒜2. w x)"
      by (subst 𝒜_split, subst sum.union_disjoint, use 𝒜_split(2) in auto)
    also have "(∑x∈𝒜1. w x) = (∑i<k. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
    proof -
      have "(∑x∈𝒜1. w x) = (∑A | A ⊆ X ∧ card A < k. ∑j∈X. w (A, j))"
        using assms by (subst sum.Sigma) (auto simp: 𝒜1_def)
      also have "… = (∑A | A ⊆ X ∧ card A < k. ∑j∈X.
                        monom (monom_of_set A) ((-1) ^ card A) * monom (sng j (k - card A)) 1)"
        unfolding w_def by (intro sum.cong) (auto simp: mult_monom)
      also have "… = (∑A | A ⊆ X ∧ card A < k. monom (monom_of_set A) ((-1) ^ card A) *
                        powsum_mpoly X (k - card A))"
        by (simp add: sum_distrib_left powsum_mpoly_altdef)
      also have "… = (∑(i,A) ∈ (SIGMA i:{..<k}. {A. A ⊆ X ∧ card A = i}).
                        monom (monom_of_set A) ((-1) ^ i) * powsum_mpoly X (k - i))"
        by (rule sum.reindex_bij_witness[of _ snd "λA. (card A, A)"]) auto
      also have "… = (∑i<k. ∑A | A ⊆ X ∧ card A = i.
                        monom (monom_of_set A) 1 * monom 0 ((-1) ^ i) * powsum_mpoly X (k - i))"
        using assms by (subst sum.Sigma) (auto simp: mult_monom)
      also have "… = (∑i<k. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
        by (simp add: sum_distrib_left sum_distrib_right mpoly_monom_0_eq_Const 
                      mpoly_Const_power mpoly_Const_uminus algebra_simps sym_mpoly_altdef)
      finally show ?thesis .
    qed
    also have "(∑x∈𝒜2. w x) = (-1) ^ k * of_nat k * sym_mpoly X k"
    proof -
      have "(∑x∈𝒜2. w x) = (∑(A,j)∈𝒜2. monom (monom_of_set A) ((- 1) ^ k))"
        by (intro sum.cong) (auto simp: 𝒜2_def w_def mpoly_monom_0_eq_Const intro!: sum.cong)
      also have "… = (∑A | A ⊆ X ∧ card A = k. ∑j∈A. monom (monom_of_set A) ((- 1) ^ k))"
        using assms finite_subset[of _ X] by (subst sum.Sigma) (auto simp: 𝒜2_def)
      also have "(λA. monom (monom_of_set A) ((- 1) ^ k) :: 'a mpoly) =
                   (λA. monom 0 ((-1) ^ k) * monom (monom_of_set A) 1)"
        by (auto simp: fun_eq_iff mult_monom)
      also have "monom 0 ((-1) ^ k) = (-1) ^ k"
        by (auto simp: mpoly_monom_0_eq_Const mpoly_Const_power mpoly_Const_uminus)
      also have "(∑A | A ⊆ X ∧ card A = k. ∑j∈A. (- 1) ^ k * monom (monom_of_set A) 1) =
                   ((-1) ^ k * of_nat k * sym_mpoly X k :: 'a mpoly)"
        by (auto simp: sum_distrib_left sum_distrib_right mult_ac sym_mpoly_altdef)
      finally show ?thesis .
    qed
    finally show ?thesis by (simp add: algebra_simps)
  qed
  
  also have "(∑x∈𝒜. w x) = 0"
  proof -
    
    define T :: "nat set × nat ⇒ nat set × nat"
      where "T = (λ(A, j). if j ∈ A then (A - {j}, j) else (insert j A, j))"
    have [simp]: "T (T x) = x" for x
      by (auto simp: T_def split: prod.splits)
    have [simp]: "T x ∈ 𝒜" if "x ∈ 𝒜" for x
    proof -
      have [simp]: "n ≤ n - Suc 0 ⟷ n = 0" for n
        by auto
      show ?thesis using that assms finite_subset[of _ X]
        by (auto simp: T_def 𝒜_def split: prod.splits)
    qed
    have "snd (T x) ∈ fst (T x) ⟷ snd x ∉ fst x" if "x ∈ 𝒜" for x
      by (auto simp: T_def split: prod.splits)
    hence bij: "bij_betw T {x∈𝒜. snd x ∈ fst x} {x∈𝒜. snd x ∉ fst x}"
      by (intro bij_betwI[of _ _ _ T]) auto
    
    have [simp]: "w (T x) = -w x" if "x ∈ 𝒜" for x
    proof -
      obtain A j where [simp]: "x = (A, j)" by force
      
      
      have aux: "w (T (A, j)) = - w (A, j)" if "(A, j) ∈ 𝒜" "j ∈ A" for j A
      proof -
        from that have [simp]: "j ∈ A" "A ⊆ X" and "k > 0"
          using finite_subset[OF _ assms, of A] by (auto simp: 𝒜_def intro!: Nat.gr0I)
        have [simp]: "finite A"
          using finite_subset[OF _ assms, of A] by auto
        from that have "card A ≤ k"
          by (auto simp: 𝒜_def)
        have card: "card A = Suc (card (A - {j}))"
          using card.remove[of A j] by auto
        hence card_less: "card (A - {j}) < card A" by linarith
        have "w (T (A, j)) = monom (monom_of_set (A - {j}) + sng j (k - card (A - {j})))
                         ((- 1) ^ card (A - {j}))" by (simp add: w_def T_def)
        also have "(- 1) ^ card (A - {j}) = ((- 1) ^ Suc (Suc (card (A - {j}))) :: 'a)"
          by simp
        also have "Suc (card (A - {j})) = card A"
          using card by simp
        also have "k - card (A - {j}) = Suc (k - card A)"
          using ‹k > 0› ‹card A ≤ k› card_less by (subst card) auto
        also have "monom_of_set (A - {j}) + sng j (Suc (k - card A)) =
                   monom_of_set A + sng j (k - card A)"
          by (transfer fixing: A j k) (auto simp: fun_eq_iff)
        also have "monom … ((-1)^ Suc (card A)) = -w (A, j)"
          by (simp add: w_def monom_uminus)
        finally show ?thesis .
      qed
      show ?thesis
      proof (cases "j ∈ A")
        case True
        with aux[of A j] that show ?thesis by auto
      next
        case False
        hence "snd (T x) ∈ fst (T x)"
          by (auto simp: T_def split: prod.splits)
        with aux[of "fst (T x)" "snd (T x)"] that show ?thesis by auto
      qed
    qed
    text ‹
      We can now show fairly easily that the sum is equal to zero.
    ›
    have *: "𝒜 = {x∈𝒜. snd x ∈ fst x} ∪ {x∈𝒜. snd x ∉ fst x}"
      by auto
    have "(∑x∈𝒜. w x) = (∑x | x ∈ 𝒜 ∧ snd x ∈ fst x. w x) + (∑x | x ∈ 𝒜 ∧ snd x ∉ fst x. w x)"
      using ‹finite 𝒜› by (subst *, subst sum.union_disjoint) auto
    also have "(∑x | x ∈ 𝒜 ∧ snd x ∉ fst x. w x) = (∑x | x ∈ 𝒜 ∧ snd x ∈ fst x. w (T x))"
      using sum.reindex_bij_betw[OF bij, of w] by simp
    also have "… = -(∑x | x ∈ 𝒜 ∧ snd x ∈ fst x. w x)"
      by (simp add: sum_negf)
    finally show "(∑x∈𝒜. w x) = 0"
      by simp
  qed
  finally show ?thesis .
qed
text ‹
  The following variant of the theorem holds for ‹k > n›. Note that this is now a
  linear recurrence relation with constant coefficients for $p_k$ in terms of
  $e_0, \ldots, e_n$.
›
corollary Girard_Newton':
  assumes "finite X" and "k > card X"
  shows   "(∑i≤card X. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i)) =
             (0 :: 'a :: comm_ring_1 mpoly)"
proof -
  have "(0 :: 'a mpoly) = (∑i<k. (- 1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
    using Girard_Newton[of X k] assms by simp
  also have "… = (∑i≤card X. (- 1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
    using assms by (intro sum.mono_neutral_right) auto
  finally show ?thesis ..
qed  
text ‹
  The following variant is the Newton--Girard Theorem solved for $e_k$, giving us
  an explicit way to determine $e_k$ from $e_0, \ldots, e_{k-1}$ and $p_1, \ldots, p_k$:
›
corollary sym_mpoly_recurrence:
  assumes k: "k > 0" and "finite X"
  shows   "(sym_mpoly X k :: 'a :: field_char_0 mpoly) =
             -smult (1 / of_nat k) (∑i=1..k. (-1) ^ i * sym_mpoly X (k - i) * powsum_mpoly X i)"
proof -
  define e p :: "nat ⇒ 'a mpoly" where [simp]: "e = sym_mpoly X" "p = powsum_mpoly X"
  have *: "0 = (-1) ^ k * of_nat k * e k +
              (∑i<k. (- 1) ^ i * e i * p (k - i) :: 'a mpoly)"
    using Girard_Newton[of X k] assms by simp
  have "0 = (-1) ^ k * smult (1 / of_nat k) (0 :: 'a mpoly)"
    by simp
  also have "… = smult (1 / of_nat k) (of_nat k) * e k +
                  smult (1 / of_nat k) (∑i<k. (-1)^(k+i) * e i * p (k - i))"
    unfolding smult_conv_mult
    using k by (subst *) (simp add: power_add sum_distrib_left sum_distrib_right field_simps 
                               del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)
  also have "smult (1 / of_nat k :: 'a) (of_nat k) = 1"
    using k by (simp add: of_nat_monom smult_conv_mult mult_monom del: monom_of_nat)
  also have "(∑i<k. (-1) ^ (k+i) * e i * p (k - i)) = (∑i=1..k. (-1) ^ i * e (k-i) * p i)"
    by (intro sum.reindex_bij_witness[of _ "λi. k - i" "λi. k - i"])
       (auto simp: minus_one_power_iff)
  finally show ?thesis unfolding e_p_def by algebra
qed
text ‹
  Analogously, the following is the theorem solved for $p_k$, giving us a
  way to determine $p_k$ from $e_0, \ldots, e_k$ and $p_1, \ldots, p_{k-1}$:
›
corollary powsum_mpoly_recurrence:
  assumes k: "k > 0" and X: "finite X"
  shows   "(powsum_mpoly X k :: 'a :: comm_ring_1 mpoly) =
             (-1) ^ (k + 1) * of_nat k * sym_mpoly X k -
             (∑i=1..<k. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
proof -
  define e p :: "nat ⇒ 'a mpoly" where [simp]: "e = sym_mpoly X" "p = powsum_mpoly X"
  have *: "0 = (-1) ^ k * of_nat k * e k +
                 (∑i<k. (-1) ^ i * e i * p (k - i) :: 'a mpoly)"
    using Girard_Newton[of X k] assms by simp
  also have "{..<k} = insert 0 {1..<k}"
    using assms by auto
  finally have "(-1) ^ k * of_nat k * e k + (∑i=1..<k. (-1) ^ i * e i * p (k - i)) + p k = 0"
    using assms by (simp add: algebra_simps)
  from add.inverse_unique[OF this] show ?thesis by simp
qed
text ‹
  Again, if we assume $k > n$, the above takes a much simpler form and is, in fact,
  a linear recurrence with constant coefficients:
›
lemma powsum_mpoly_recurrence':
  assumes k: "k > card X" and X: "finite X"
  shows   "(powsum_mpoly X k :: 'a :: comm_ring_1 mpoly) =
             -(∑i=1..card X. (-1) ^ i * sym_mpoly X i * powsum_mpoly X (k - i))"
proof -
  define e p :: "nat ⇒ 'a mpoly" where [simp]: "e = sym_mpoly X" "p = powsum_mpoly X"
  have "p k = (-1) ^ (k + 1) * of_nat k * e k - (∑i=1..<k. (-1) ^ i * e i * p (k - i))"
    unfolding e_p_def using assms by (intro powsum_mpoly_recurrence) auto
  also have "… = -(∑i=1..<k. (-1) ^ i * e i * p (k - i))"
    using assms by simp
  also have "(∑i=1..<k. (-1) ^ i * e i * p (k - i)) = (∑i=1..card X. (-1) ^ i * e i * p (k - i))"
    using assms by (intro sum.mono_neutral_right) auto
  finally show ?thesis by simp
qed
end