(* File: Power_Sum_Puzzle.thy Author: Manuel Eberl, TU München *) section ‹Power sum puzzles› theory Power_Sum_Puzzle imports Power_Sum_Polynomials "Polynomial_Factorization.Rational_Root_Test" begin subsection ‹General setting and results› text ‹ We now consider the following situation: Given unknown complex numbers $x_1,\ldots,x_n$, define $p_k = x_1^k + \ldots + x_n^k$. Also, define $e_k := e_k(x_1,\ldots,x_n)$ where $e_k(X_1,\ldots,X_n)$ is the $k$-th elementary symmetric polynomial. What is the relationship between the sequences $e_k$ and $p_k$; in particular, how can we determine one from the other? › locale power_sum_puzzle = fixes x :: "nat ⇒ complex" fixes n :: nat begin text ‹ We first introduce the notation $p_k := x_1 ^ k + \ldots + x_n ^ k$: › definition p where "p k = (∑i<n. x i ^ k)" lemma p_0 [simp]: "p 0 = of_nat n" by (simp add: p_def) lemma p_altdef: "p k = insertion x (powsum_mpoly {..<n} k)" by (simp add: p_def) text ‹ Similarly, we introduce the notation $e_k = e_k(x_1,\ldots, x_n)$ where $e_k(X_1,\ldots,X_n)$ is the $k$-th elementary symmetric polynomial (i.\,e. the sum of all monomials that can be formed by taking the product of exactly $k$ distinct variables). › definition e where "e k = (∑Y | Y ⊆ {..<n} ∧ card Y = k. prod x Y)" lemma e_altdef: "e k = insertion x (sym_mpoly {..<n} k)" by (simp add: e_def insertion_sym_mpoly) text ‹ It is clear that $e_k$ vanishes for $k > n$. › lemma e_eq_0 [simp]: "k > n ⟹ e k = 0" by (simp add: e_altdef) lemma e_0 [simp]: "e 0 = 1" by (simp add: e_altdef) text ‹ The recurrences we got from the Girard--Newton Theorem earlier now directly give us analogous recurrences for $e_k$ and $p_k$: › lemma e_recurrence: assumes k: "k > 0" shows "e k = -(∑i=1..k. (- 1) ^ i * e (k - i) * p i) / of_nat k" using assms unfolding e_altdef p_altdef by (subst sym_mpoly_recurrence) (auto simp: insertion_sum insertion_add insertion_mult insertion_power insertion_sym_mpoly) lemma p_recurrence: assumes k: "k > 0" shows "p k = -of_nat k * (-1) ^ k * e k - (∑i=1..<k. (-1) ^ i * e i * p (k - i))" using assms unfolding e_altdef p_altdef by (subst powsum_mpoly_recurrence) (auto simp: insertion_sum insertion_add insertion_mult insertion_diff insertion_power insertion_sym_mpoly) lemma p_recurrence'': assumes k: "k > n" shows "p k = -(∑i=1..n. (-1) ^ i * e i * p (k - i))" using assms unfolding e_altdef p_altdef by (subst powsum_mpoly_recurrence') (auto simp: insertion_sum insertion_add insertion_mult insertion_diff insertion_power insertion_sym_mpoly) text ‹ It is clear from this recurrence that if $p_1$ to $p_n$ are rational, then so are the $e_k$: › lemma e_in_Rats: assumes "⋀k. k ∈ {1..n} ⟹ p k ∈ ℚ" shows "e k ∈ ℚ" proof (cases "k ≤ n") case True thus ?thesis proof (induction k rule: less_induct) case (less k) show ?case proof (cases "k = 0") case False thus ?thesis using assms less by (subst e_recurrence) (auto intro!: Rats_divide) qed auto qed qed auto text ‹ Analogously, if $p_1$ to $p_n$ are rational, then so are all the other $p_k$: › lemma p_in_Rats: assumes "⋀k. k ∈ {1..n} ⟹ p k ∈ ℚ" shows "p k ∈ ℚ" proof (induction k rule: less_induct) case (less k) consider "k = 0" | "k ∈ {1..n}" | "k > n" by force thus ?case proof cases assume "k > n" thus ?thesis using less assms by (subst p_recurrence'') (auto intro!: sum_in_Rats Rats_mult e_in_Rats) qed (use assms in auto) qed text ‹ Next, we define the unique monic polynomial that has $x_1, \ldots, x_n$ as its roots (respecting multiplicity): › definition Q :: "complex poly" where "Q = (∏i<n. [:-x i, 1:])" lemma degree_Q [simp]: "Polynomial.degree Q = n" by (simp add: Q_def degree_prod_eq_sum_degree) lemma lead_coeff_Q [simp]: "Polynomial.coeff Q n = 1" using monic_prod[of "{..<n}" "λi. [:-x i, 1:]"] by (simp add: Q_def degree_prod_eq_sum_degree) text ‹ By Vieta's Theorem, we then have: \[Q(X) = \sum_{k=0}^n (-1)^{n-k} e_{n-k} X^k\] In other words: The above allows us to determine the $x_1, \ldots, x_n$ explicitly. They are, in fact, precisely the roots of the above polynomial (respecting multiplicity). Since this polynomial depends only on the $e_k$, which are in turn determined by $p_1, \ldots, p_n$, this means that these are the ∗‹only› solutions of this puzzle (up to permutation of the $x_i$). › lemma coeff_Q: "Polynomial.coeff Q k = (if k > n then 0 else (-1) ^ (n - k) * e (n - k))" proof (cases "k ≤ n") case True thus ?thesis using coeff_poly_from_roots[of "{..<n}" k x] by (auto simp: Q_def e_def) qed (auto simp: Polynomial.coeff_eq_0) lemma Q_altdef: "Q = (∑k≤n. Polynomial.monom ((-1) ^ (n - k) * e (n - k)) k)" by (subst poly_as_sum_of_monoms [symmetric]) (simp add: coeff_Q) text ‹ The following theorem again shows that $x_1, \ldots, x_n$ are precisely the roots of \<^term>‹Q›, respecting multiplicity. › theorem mset_x_eq_poly_roots_Q: "{#x i. i ∈# mset_set {..<n}#} = poly_roots Q" proof - have "poly_roots Q = (∑i<n. {#x i#})" by (simp add: Q_def poly_roots_prod) also have "… = {#x i. i ∈# mset_set {..<n}#}" by (induction n) (auto simp: lessThan_Suc) finally show ?thesis .. qed end subsection ‹Existence of solutions› text ‹ So far, we have assumed a solution to the puzzle and then shown the properties that this solution must fulfil. However, we have not yet shown that there ∗‹is› a solution. We will do that now. Let $n$ be a natural number and $f_k$ some sequence of complex numbers. We will show that there are $x_1, \ldots, x_n$ so that $x_1 ^ k + \ldots + x_n ^ k = f_k$ for any $1\leq k\leq n$. › locale power_sum_puzzle_existence = fixes f :: "nat ⇒ complex" and n :: nat begin text ‹ First, we define a sequence of numbers ‹e'› analogously to the sequence ‹e› before, except that we replace all occurrences of the power sum $p_k$ with $f_k$ (recall that in the end we want $p_k = f_k$). › fun e' :: "nat ⇒ complex" where "e' k = (if k = 0 then 1 else if k > n then 0 else -(∑i=1..k. (-1) ^ i * e' (k - i) * f i) / of_nat k)" lemmas [simp del] = e'.simps lemma e'_0 [simp]: "e' 0 = 1" by (simp add: e'.simps) lemma e'_eq_0 [simp]: "k > n ⟹ e' k = 0" by (auto simp: e'.simps) text ‹ Just as before, we can show the following recurrence for ‹f› in terms of ‹e'›: › lemma f_recurrence: assumes k: "k > 0" "k ≤ n" shows "f k = -of_nat k * (-1) ^ k * e' k - (∑i=1..<k. (- 1) ^ i * e' i * f (k - i))" proof - have "-of_nat k * e' k = (∑i=1..k. (- 1) ^ i * e' (k - i) * f i)" using assms by (subst e'.simps) (simp add: field_simps) hence "(-1)^k * (-of_nat k * e' k) = (-1)^k * (∑i=1..k. (- 1) ^ i * e' (k - i) * f i)" by simp also have "… = f k + (-1) ^ k * (∑i=1..<k. (- 1) ^ i * e' (k - i) * f i)" using assms by (subst sum.last_plus) (auto simp: minus_one_power_iff) also have "(-1) ^ k * (∑i=1..<k. (- 1) ^ i * e' (k - i) * f i) = (∑i=1..<k. (- 1) ^ (k - i) * e' (k - i) * f i)" unfolding sum_distrib_left by (intro sum.cong) (auto simp: minus_one_power_iff) also have "… = (∑i=1..<k. (- 1) ^ i * e' i * f (k - i))" by (intro sum.reindex_bij_witness[of _ "λi. k - i" "λi. k - i"]) auto finally show ?thesis by (simp add: algebra_simps) qed text ‹ We now define a polynomial whose roots will be precisely the solution $x_1, \ldots, x_n$ to our problem. › lift_definition Q' :: "complex poly" is "λk. if k > n then 0 else (-1) ^ (n - k) * e' (n - k)" using eventually_gt_at_top[of n] unfolding cofinite_eq_sequentially by eventually_elim auto lemma coeff_Q': "Polynomial.coeff Q' k = (if k > n then 0 else (-1) ^ (n - k) * e' (n - k))" by transfer auto lemma lead_coeff_Q': "Polynomial.coeff Q' n = 1" by (simp add: coeff_Q') lemma degree_Q' [simp]: "Polynomial.degree Q' = n" proof (rule antisym) show "Polynomial.degree Q' ≥ n" by (rule le_degree) (auto simp: coeff_Q') show "Polynomial.degree Q' ≤ n" by (rule degree_le) (auto simp: coeff_Q') qed text ‹ Since the complex numbers are algebraically closed, this polynomial splits into linear factors: › definition Root :: "nat ⇒ complex" where "Root = (SOME Root. Q' = (∏i<n. [:-Root i, 1:]))" lemma Root: "Q' = (∏i<n. [:-Root i, 1:])" proof - obtain rs where rs: "(∏r←rs. [:-r, 1:]) = Q'" "length rs = n" using fundamental_theorem_algebra_factorized[of Q'] lead_coeff_Q' by auto have "Q' = (∏r←rs. [:-r, 1:])" by (simp add: rs) also have "… = (∏r=0..<n. [:-rs ! r, 1:])" by (subst prod_list_prod_nth) (auto simp: rs) also have "{0..<n} = {..<n}" by auto finally have "∃Root. Q' = (∏i<n. [:-Root i, 1:])" by blast thus ?thesis unfolding Root_def by (rule someI_ex) qed text ‹ We can therefore now use the results from before for these $x_1, \ldots, x_n$. › sublocale power_sum_puzzle Root n . text ‹ Vieta's theorem gives us an expression for the coefficients of ‹Q'› in terms of $e_k(x_1,\ldots,x_n)$. This shows that our ‹e'› is indeed exactly the same as ‹e›. › lemma e'_eq_e: "e' k = e k" proof (cases "k ≤ n") case True from True have "e' k = (-1) ^ k * poly.coeff Q' (n - k)" by (simp add: coeff_Q') also have "Q' = (∏x<n. [:-Root x, 1:])" using Root by simp also have "(-1) ^ k * poly.coeff … (n - k) = e k" using True coeff_poly_from_roots[of "{..<n}" "n - k" Root] by (simp add: insertion_sym_mpoly e_altdef) finally show "e' k = e k" . qed auto text ‹ It then follows by a simple induction that $p_k = f_k$ for $1\leq k\leq n$, as intended: › lemma p_eq_f: assumes "k > 0" "k ≤ n" shows "p k = f k" using assms proof (induction k rule: less_induct) case (less k) thus "p k = f k" using p_recurrence[of k] f_recurrence[of k] less by (simp add: e'_eq_e) qed end text ‹ Here is a more condensed form of the above existence theorem: › theorem power_sum_puzzle_has_solution: fixes f :: "nat ⇒ complex" shows "∃Root. ∀k∈{1..n}. (∑i<n. Root i ^ k) = f k" proof - interpret power_sum_puzzle_existence f . from p_eq_f have "∀k∈{1..n}. (∑i<n. Root i ^ k) = f k" by (auto simp: p_def) thus ?thesis by blast qed subsection ‹A specific puzzle› text ‹ We now look at one particular instance of this puzzle, which was given as an exercise in ∗‹Abstract Algebra› by Dummit and Foote (Exercise 23 in Section 14.6)~\<^cite>‹"dummit"›. Suppose we know that $x + y + z = 1$, $x^2 + y^2 + z^2 = 2$, and $x^3 + y^3 + z^3 = 3$. Then what is $x^5+y^5+z^5$? What about any arbitrary $x^n+y^n+z^n$? › locale power_sum_puzzle_example = fixes x y z :: complex assumes xyz: "x + y + z = 1" "x^2 + y^2 + z^2 = 2" "x^3 + y^3 + z^3 = 3" begin text ‹ We reuse the results we have shown in the general case before. › definition f where "f n = [x,y,z] ! n" sublocale power_sum_puzzle f 3 . text ‹ We can simplify \<^term>‹p› a bit more now. › lemma p_altdef': "p k = x ^ k + y ^ k + z ^ k" unfolding p_def f_def by (simp add: eval_nat_numeral) lemma p_base [simp]: "p (Suc 0) = 1" "p 2 = 2" "p 3 = 3" using xyz by (simp_all add: p_altdef') text ‹ We can easily compute all the non-zero values of \<^term>‹e› recursively: › lemma e_Suc_0 [simp]: "e (Suc 0) = 1" by (subst e_recurrence; simp) lemma e_2 [simp]: "e 2 = -1/2" by (subst e_recurrence; simp add: atLeastAtMost_nat_numeral) lemma e_3 [simp]: "e 3 = 1/6" by (subst e_recurrence; simp add: atLeastAtMost_nat_numeral) text ‹ Plugging in all the values, the recurrence relation for \<^term>‹p› now looks like this: › lemma p_recurrence''': "k > 3 ⟹ p k = p (k-3) / 6 + p (k-2) / 2 + p (k-1)" using p_recurrence''[of k] by (simp add: atLeastAtMost_nat_numeral) text ‹ Also note again that all $p_k$ are rational: › lemma p_in_Rats': "p k ∈ ℚ" proof - have *: "{1..3} = {1, 2, (3::nat)}" by auto also have "∀k∈…. p k ∈ ℚ" by auto finally show ?thesis using p_in_Rats[of k] by simp qed text ‹ The above recurrence has the characteristic polynomial $X^3 - X^2 - \frac{1}{2} X - \frac{1}{6}$ (which is exactly our \<^term>‹Q›), so we know that can now specify $x$, $y$, and $z$ more precisely: They are the roots of that polynomial (in unspecified order). › lemma xyz_eq: "{#x, y, z#} = poly_roots [:-1/6, -1/2, -1, 1:]" proof - have "image_mset f (mset_set {..<3}) = poly_roots Q" using mset_x_eq_poly_roots_Q . also have "image_mset f (mset_set {..<3}) = {#x, y, z#}" by (simp add: numeral_3_eq_3 lessThan_Suc f_def Multiset.union_ac) also have "Q = [:-1/6, -1/2, -1, 1:]" by (simp add: Q_altdef atMost_nat_numeral Polynomial.monom_altdef power3_eq_cube power2_eq_square) finally show ?thesis . qed text ‹ Using the rational root test, we can easily show that $x$, $y$, and $z$ are irrational. › lemma xyz_irrational: "set_mset (poly_roots [:-1/6, -1/2, -1, 1::complex:]) ∩ ℚ = {}" proof - define p :: "rat poly" where "p = [:-1/6, -1/2, -1, 1:]" have "rational_root_test p = None" unfolding p_def by code_simp hence "¬(∃x::rat. poly p x = 0)" by (rule rational_root_test) hence "¬(∃x∈ℚ. poly (map_poly of_rat p) x = (0 :: complex))" by (auto simp: Rats_def) also have "map_poly of_rat p = [:-1/6, -1/2, -1, 1 :: complex:]" by (simp add: p_def of_rat_minus of_rat_divide) finally show ?thesis by auto qed text ‹ This polynomial is ∗‹squarefree›, so these three roots are, in fact, unique (so that there are indeed $3! = 6$ possible permutations). › lemma rsquarefree: "rsquarefree [:-1/6, -1/2, -1, 1 :: complex:]" by (rule coprime_pderiv_imp_rsquarefree) (auto simp: pderiv_pCons coprime_iff_gcd_eq_1 gcd_poly_code gcd_poly_code_def content_def primitive_part_def gcd_poly_code_aux_reduce pseudo_mod_def pseudo_divmod_def Let_def Polynomial.monom_altdef normalize_poly_def) lemma distinct_xyz: "distinct [x, y, z]" by (rule rsquarefree_imp_distinct_roots[OF rsquarefree]) (simp_all add: xyz_eq) text ‹ While these roots ∗‹can› be written more explicitly in radical form, they are not very pleasant to look at. We therefore only compute a few values of ‹p› just for fun: › lemma "p 4 = 25 / 6" and "p 5 = 6" and "p 10 = 15539 / 432" by (simp_all add: p_recurrence''') text ‹ Lastly, let us (informally) examine the asymptotics of this problem. Two of the roots have a norm of roughly $\beta \approx 0.341$, while the remaining root ‹α› is roughly 1.431. Consequently, $x^n + y^n + z^n$ is asymptotically equivalent to $\alpha^n$, with the error being bounded by $2\cdot \beta^n$ and therefore goes to 0 very quickly. For $p(10) = \frac{15539}{432} \approx 35.97$, for instance, this approximation is correct up to 6 decimals (a relative error of about 0.0001\,\%). › end text ‹ To really emphasise that the above puzzle has a solution and the locale is not `vacuous', here is an interpretation of the locale using the existence theorem from before: › notepad begin define f :: "nat ⇒ complex" where "f = (λk. [1,2,3] ! (k - 1))" obtain Root :: "nat ⇒ complex" where Root: "⋀k. k ∈ {1..3} ⟹ (∑i<3. Root i ^ k) = f k" using power_sum_puzzle_has_solution[of 3 f] by metis define x y z where "x = Root 0" "y = Root 1" "z = Root 2" have "x + y + z = 1" and "x^2 + y^2 + z^2 = 2" and "x^3 + y^3 + z^3 = 3" using Root[of 1] Root[of 2] Root[of 3] by (simp_all add: eval_nat_numeral x_y_z_def f_def) then interpret power_sum_puzzle_example x y z by unfold_locales have "p 5 = 6" by (simp add: p_recurrence''') end end