Generalizing the sum of consecutive cubes $\sum_{k=1}^n k^3 = \Big(\sum_{k=1}^n k\Big)^2$ to other odd powers

We have,

$$\sum_{k=1}^n k^3 = \Big(\sum_{k=1}^n k\Big)^2$$

$$2\sum_{k=1}^n k^5 = -\Big(\sum_{k=1}^n k\Big)^2+3\Big(\sum_{k=1}^n k^2\Big)^2$$

$$2\sum_{k=1}^n k^7 = \Big(\sum_{k=1}^n k\Big)^2-3\Big(\sum_{k=1}^n k^2\Big)^2+4\Big(\sum_{k=1}^n k^3\Big)^2$$

and so on (apparently).

Is it true that the sum of consecutive odd $m$ powers, for $m>1$, can be expressed as sums of squares of sums* in a manner similar to the above? What is the general formula?

*(Edited re Lord Soth's and anon's comment.)

• See Faulhaber's formula. As an historical aside, many famous mathematicians and scientists (including Richard Feynman) have been known to pose this problem to themselves in their youth and independently discover the Bernoulli numbers. – David H Nov 3 '13 at 5:54
• Ah, thanks. Didn't even know it had a name. The section on Faulhaber polynomials is close, as it discusses odd $m$ powers, but it does not express them as sums of squares. Perhaps a modification of the formula there can provide the answer. – Tito Piezas III Nov 3 '13 at 6:02
• What do you mean "in a similar manner?" For example, for the $7$th power, you have expressed in terms of the squares of sums of first, second and third powers. Perhaps you meant "squares of sums" instead of "sums of squares." - And - I guess - you want to express the sum of the $2k+1$th power in terms of only the squares of sums of $1,2,\ldots,k$th powers. – Lord Soth Nov 3 '13 at 6:03
• @LordSoth If you want to be technical, it's a sum of squares of sums. Which, in particular, is a sum of squares. Like so: $$\sum_{k=1}^n k^{2m+1}=\sum_{t=1}^m c_t\left(\sum_{\ell=1}^n\ell^t\right)^2.$$ – anon Nov 3 '13 at 6:05
• @anon: Exactly. :) – Tito Piezas III Nov 3 '13 at 6:05

This is a partial answer, it just establishes the existence.

We have $$s_m(n) = \sum_{k=1}^n k^m = \frac{1}{m+1} \left(\operatorname{B}_{m+1}(n+1)-\operatorname{B}_{m+1}(1)\right)$$ where $\operatorname{B}_m(x)$ denotes the monic Bernoulli polynomial of degree $m$, which has the following useful properties: \begin{align} \int_x^{x+1}\operatorname{B}_m(t)\,\mathrm{d}t &= x^m \quad\text{(from which everything else follows)} \\\operatorname{B}_{m+1}'(x) &= (m+1)\operatorname{B}_m(x) \\\operatorname{B}_m\left(x+\frac{1}{2}\right) &\begin{cases} \text{is even in x} & \text{for even m} \\ \text{is odd in x} & \text{for odd m} \end{cases} \\\operatorname{B}_m(0) = \operatorname{B}_m(1) &= 0 \quad\text{for odd m\geq3} \end{align}

Therefore, \begin{align} s_m(n) &\text{has degree m+1 in n} \\s_m(0) &= 0 \\s_m'(0) &= \operatorname{B}_m(1) = 0\quad\text{for odd m\geq3} \\&\quad\text{(This makes n=0 a double zero of s_m(n) for odd m\geq3)} \\s_m\left(x-\frac{1}{2}\right) &\begin{cases} \text{is even in x} & \text{for odd m} \\ \text{is odd in x} & \text{for even m\geq2} \end{cases} \end{align}

Consider the vector space $V_m$ of univariate polynomials $\in\mathbb{Q}[x]$ with degree not exceeding $2m+2$, that are even in $x$ and have a double zero at $x=\frac{1}{2}$. Thus $V_m$ has dimension $m$ and is clearly spanned by $$\left\{s_j^2\left(x-\frac{1}{2}\right)\mid j=1,\ldots,m\right\}$$ For $m>0$, we find that $s_{2m+1}(x-\frac{1}{2})$ has all the properties required for membership in $V_m$. Substituting $x-\frac{1}{2}=n$, we conclude that there exists a representation $$s_{2m+1}(n) = \sum_{j=1}^m a_{m,j}\,s_j^2(n) \quad\text{for m>0 with a_{m,j}\in\mathbb{Q}}$$ of $s_{2m+1}(n)$ as a linear combination of squares of sums.

• Out of curiosity, do you suppose you can find the $a_m$ in terms of the Bernoulli numbers? :) – Tito Piezas III Nov 8 '13 at 3:58
• @TitoPiezasIII: I am looking into it, but I do not know yet. – ccorn Nov 8 '13 at 10:53

First step was to build the matrix $A$ of coefficients of the squares of the Faulhaber-polynomials. Next step, to build the matrix $B$ of the coefficients of the un-squared Faulhaber-polynomials, but only of the odd orders.

Then the matrix $C$ of coefficients for the composition, such that $C \cdot A = B$ occurs by inversion of $A$ such that $C = B \cdot A^{-1}$. This can be done with any truncation of size. Important: Note that the first two columns of $A$ and $B$ are completely zero, so we must reduce that matrices to the remaining columns.

The inverse of the coefficients-matrix $C$ (which itself can be seen in the inner braces of Hecke's answer) might also be interesting; at least it has a slightly simpler structure: $$C^{-1}= \tiny \begin{bmatrix} 2 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1/3 & 2/3 & . & . & . & . & . & . & . & . & . & . & . \\ . & . & 1/2 & 1/2 & . & . & . & . & . & . & . & . & . & . \\ . & . & -1/15 & 2/3 & 2/5 & . & . & . & . & . & . & . & . & . \\ . & . & . & -1/6 & 5/6 & 1/3 & . & . & . & . & . & . & . & . \\ . & . & . & 1/21 & -1/3 & 1 & 2/7 & . & . & . & . & . & . & . \\ . & . & . & . & 1/6 & -7/12 & 7/6 & 1/4 & . & . & . & . & . & . \\ . & . & . & . & -1/15 & 4/9 & -14/15 & 4/3 & 2/9 & . & . & . & . & . \\ . & . & . & . & . & -3/10 & 1 & -7/5 & 3/2 & 1/5 & . & . & . & . \\ . & . & . & . & . & 5/33 & -1 & 2 & -2 & 5/3 & 2/11 & . & . & . \\ . & . & . & . & . & . & 5/6 & -11/4 & 11/3 & -11/4 & 11/6 & 1/6 & . & . \\ . & . & . & . & . & . & -691/1365 & 10/3 & -33/5 & 44/7 & -11/3 & 2 & 2/13 & . \\ . & . & . & . & . & . & . & -691/210 & 65/6 & -143/10 & 143/14 & -143/30 & 13/6 & 1/7 \end{bmatrix}$$ ... and the numerator of $691$ in rows $13,14,15,\ldots$ surely rings some bell concerning Bernoulli-numbers
Note: this matrix gives the compositions of the squared sums-of-powers by the higher odd orders of the original sums-of-powers, so helps for the inverse of your problem: $$\begin{array}{rrrrr} S(0,x)^2 &=& 2 S(1,x) & - 1 S(0,x) &\\ \hline \\ S(1,x)^2 &=& 1 S(3,x) \\ S(2,x)^2 &=& 1/3 \; S(3,x) &+2/3 \; S(5,x) \\ S(3,x)^2 &=&& 1/2 \; S(5,x) &+1/2 \; S(7,x) \\ S(4,x)^2 &=&& -1/15 \; S(5,x)&+2/3 \; S(7,x)&+2/5 \; S(9,x) \\ \vdots \end{array}$$ The first row does not fit into the pattern. The coefficients are compositions of binomials and Bernoulli-numbers, which might best be written as below according to the analysis of the matrix $E$

[update] The following matrix displays coefficients after the Bernoulli-numbers and the column-index are removed. $$D = \Tiny \begin{bmatrix} 2 & . & . & . & . & . & . & . & . & . & . & . \\ . & 2 & . & . & . & . & . & . & . & . & . & . \\ . & 4 & 2 & . & . & . & . & . & . & . & . & . \\ . & . & 9 & 2 & . & . & . & . & . & . & . & . \\ . & . & 6 & 16 & 2 & . & . & . & . & . & . & . \\ . & . & . & 20 & 25 & 2 & . & . & . & . & . & . \\ . & . & . & 8 & 50 & 36 & 2 & . & . & . & . & . \\ . & . & . & . & 35 & 105 & 49 & 2 & . & . & . & . \\ . & . & . & . & 10 & 112 & 196 & 64 & 2 & . & . & . \\ . & . & . & . & . & 54 & 294 & 336 & 81 & 2 & . & . \\ . & . & . & . & . & 12 & 210 & 672 & 540 & 100 & 2 & . \\ . & . & . & . & . & . & 77 & 660 & 1386 & 825 & 121 & 2 \end{bmatrix}$$ We have $$C^{-1}_{r,c} = D_{r,c}\cdot {B_{2(r-c)} \over c }$$ where $B_k$ are the Bernoulli-numbers and the column- and rowindices $c,r$ begin at $1$.

Perhaps an even better display is this: $$E= \Tiny \begin{bmatrix} 1 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 3 & 1 & . & . & . & . & . & . & . & . & . & . & . \\ . & . & 6 & 1 & . & . & . & . & . & . & . & . & . & . \\ . & . & 5 & 10 & 1 & . & . & . & . & . & . & . & . & . \\ . & . & . & 15 & 15 & 1 & . & . & . & . & . & . & . & . \\ . & . & . & 7 & 35 & 21 & 1 & . & . & . & . & . & . & . \\ . & . & . & . & 28 & 70 & 28 & 1 & . & . & . & . & . & . \\ . & . & . & . & 9 & 84 & 126 & 36 & 1 & . & . & . & . & . \\ . & . & . & . & . & 45 & 210 & 210 & 45 & 1 & . & . & . & . \\ . & . & . & . & . & 11 & 165 & 462 & 330 & 55 & 1 & . & . & . \\ . & . & . & . & . & . & 66 & 495 & 924 & 495 & 66 & 1 & . & . \\ . & . & . & . & . & . & 13 & 286 & 1287 & 1716 & 715 & 78 & 1 & . \\ . & . & . & . & . & . & . & 91 & 1001 & 3003 & 3003 & 1001 & 91 & 1 \end{bmatrix}$$ where $C^{-1}_{r,c} = E_{r,c} \cdot 2 B_{2(r-c)} / r$ and it seems $E_{r,c} = \binom{r}{2(r-c)}$ where $c \le r \lt 2c$ so it seems $$C^{-1}_{r,c} = B_{2(r-c)} \cdot \binom{r}{2(r-c)}\cdot \frac 2{2c-r+1} \qquad \text{ where } c \le r \lt 2c \tag 1$$ [added]
If I use this description of the composition of $C^{-1}$ then I can symbolically invert this again and get (the top-left range here) $$C= \Tiny \begin{bmatrix} 1/2 & . & . & . & . & . \\ 0 & 1 & . & . & . & . \\ 0 & -3 \; B_2 & 3/2 & . & . & . \\ 0 & 18 \; B_2^2 & -9 \; B_2 & 2 & . & . \\ 0 & -180 \; B_2^3+15 \; B_4 \; B_2 & 90 \; B_2^2-15/2 \; B_4 & -20 \; B_2 & 5/2 & . \\ 0 & 2700 \; B_2^4-495 \; B_4 \; B_2^2 & -1350 \; B_2^3+495/2 \; B_4 \; B_2 & 300 \; B_2^2-30 \; B_4 & -75/2 \; B_2 & 3 \end{bmatrix} \tag 2$$ and where the first row/column should be deleted since they are out of the pattern. For example, using the third and the fourth row in this matrix we get $$\begin{array} {rrlll} S(5,x) &=& -3 \; B_2 \; S(1,x)^2 &+ 3/2 \; S(2,x)^2 \\ S(7,x) &=& 18 \; B_2^2 \; S(1,x)^2 &- 9 \; B_2 \; S(2,x)^2 &+ 2 \; S(3,x)^2 \end{array} \tag 3$$

• Intriguing. I hope someone can build on your observation. I'm really curious to know how the $a_i$ look like in a formula that uses the Bernoulli numbers. – Tito Piezas III Nov 11 '13 at 2:29
• You can find a formula for representing polynomials as linear combinations of Bernoulli polynomials by interpreting the integral equation for the Bernoulli polynomials as a linear integral transform that maps $\operatorname{B}_m(x)$ to $x^m$, then apply the formula for the coefficients of a Maclaurin series. Use this on $\left(\operatorname{B}_m(x)-\operatorname{B}_m\right)^2$, apply the General Leibniz rule for derivatives of products and use the fact that $\operatorname{B}_k(0)\neq\operatorname{B}_k(1)$ only for $k=1$. This reproduces @GottfriedHelms's result on the coefficients of $C^{-1}$. – ccorn Nov 12 '13 at 3:31

Not an answer, but only some results I found.

Denote $s_m=\sum_{k=1}^nk^m.$ Suppose that $$s_{2m+1}=\sum_{i=1}^m a_is_i^2,a_i\in\mathbb Q.$$

Use the method of undetermined coefficients, we can get a list of $\{m,\{a_i\}\}$: $$\begin{array}{l} \{1,\{1\}\} \\ \left\{2,\left\{-\frac{1}{2},\frac{3}{2}\right\}\right\} \\ \left\{3,\left\{\frac{1}{2},-\frac{3}{2},2\right\}\right\} \\ \left\{4,\left\{-\frac{11}{12},\frac{11}{4},-\frac{10}{3},\frac{5}{2}\right\}\right\} \\ \left\{5,\left\{\frac{61}{24},-\frac{61}{8},\frac{28}{3},-\frac{25}{4},3\right\}\right\} \\ \left\{6,\left\{-\frac{1447}{144},\frac{1447}{48},-\frac{332}{9},\frac{595}{24},-\frac{21}{2},\frac{7}{2}\right\}\right\} \\ \left\{7,\left\{\frac{5771}{108},-\frac{5771}{36},\frac{5296}{27},-\frac{2375}{18},56,-\frac{49}{3},4\right\}\right\} \\ \left\{8,\left\{-\frac{53017}{144},\frac{53017}{48},-\frac{60817}{45},\frac{21817}{24},-\frac{3861}{10},\frac{1127}{10},-24,\frac{9}{2}\right\}\right\} \\ \left\{9,\left\{\frac{2755645}{864},-\frac{2755645}{288},\frac{632213}{54},-\frac{1133965}{144},\frac{13379}{4},-\frac{11725}{12},208,-\frac{135}{4},5\right\}\right\}\end{array}$$

Through some observation, we can find that $a_i a_{i+1}<0,a_m=\dfrac{m+1}2,a_2=-3a_1,a_{m-1}=\dfrac{1}{24} m^2 (m+1).$

• Last night I solved it up to 9th powers ($m=4$ in you table) and, since I was sleepy, assumed the pattern of sums of squares of sums would hold true for all odd powers $> 1$. If there is a general formula, it may involve the Bernoulli numbers (as other posters have noted). – Tito Piezas III Nov 3 '13 at 18:32
• On a hunch, I noticed that $-1/2+3/2 = 1$ and $1/2-3/2+2 = 1$. Turns out ALL have the sum of their coefficients as $\sum a_i = 1$. (P.S. The same thing happens with the formulas for sums of consecutive powers.) – Tito Piezas III Nov 4 '13 at 4:02
• @Tito Piezas III If we let $n=1$ in the formulas then we get $\sum a_i=1.$ And this is how I get the $a_i$ in the list (let $n=1,2,\cdots m$). – lsr314 Nov 4 '13 at 4:08
• Yes, I realized that. :) – Tito Piezas III Nov 4 '13 at 4:10
• @TitoPiezasIII Some experimentation and algebra confirms that $a_m=\frac{m+1}{2}, a_{m-1}=-\frac{1}{24}m^2(m+1)$ (Hecke missed a minus sign) Also, for $m \geq 4$, we have $$a_{m-2}=\frac{1}{720}(m+1)(3m-1)m(m-1)^2$$ and for $m \geq 6$, we have $$a_{m-3}=-\frac{1}{120960}(m-1)(m+1)(3m-2)(17m-19)m(m-2)^2$$ Would do more, but need to rush off for lectures. – Ivan Loh Nov 4 '13 at 9:50

Too long for a comment on @Gottfried's answer. (Note: This was written before the answer was edited to include notes about the matrix $E$.)

In matrix $D$, the main diagonal consists of $2$s, and the next consists of squares: \begin{align} D_{r,r\phantom{-1}} &= 2 & r \geq 1\\ D_{r,r-1} &= (r-1)^2 = \binom {r-1}{1}\frac{r-1}{1} & r \geq 3 \end{align}

The third diagonal consists of consecutive 4-dimensional pyramidal numbers (OEIS #A002415); the next, consecutive 6-dimensional square numbers (OEIS #A040977); the next, consecutive 8-dimensional square numbers (OEIS #A053347).

If I've done my index arithmetic correctly (please double-check), we have \begin{align} D_{r,r-2} &= \binom{r-1}{3}\frac{r-2}{2} & r \geq 5 \\ D_{r,r-3} &= \binom{r-1}{5}\frac{r-3}{3} & r \geq 7 \\ D_{r,r-4} &= \binom{r-1}{7}\frac{r-4}{4} & r \geq 9 \end{align} where row and column indices start at $1$. The two entries of the final visible column satisfy $$D_{r,r-5} = \binom{r-1}{9}\frac{r-5}{5} \qquad r \geq 11$$

So, for $r\neq c$, the non-zero entries of $D$ appear to be $$D_{r,c} = \binom{r-1}{2(r-c)-1}\frac{c}{r-c}$$

This is not an answer but does not fit in a comment-box.

Another way to describe the relation between the sums and the sums of squares of various sums-of-like-powers makes use of the Eulerian-numbers, which converts the expressions for the sums-of-like-powers into polynomials. The sums of like powers can be expressed as $$\small \begin{array} {rclllll} s_0(n) = & 1 \cdot \binom{n}{1} \\ s_1(n) = & & 1 \cdot \binom{n+1}{2} \\ s_2(n) = & & 1 \cdot \binom{n+1}{3}& +1 \cdot \binom{n+2}{3} \\ s_3(n) = & & 1 \cdot \binom{n+1}{4}& +4 \cdot \binom{n+2}{4} & +1 \cdot \binom{n+3}{4} \\ s_4(n) = & & 1 \cdot \binom{n+1}{5}& +11 \cdot \binom{n+2}{5} & +11 \cdot \binom{n+3}{5} & +1 \cdot \binom{n+4}{5} \\ s_5(n) = & & 1 \cdot \binom{n+1}{6}& +26 \cdot \binom{n+2}{6}& +66 \cdot \binom{n+3}{6} & +26 \cdot \binom{n+4}{6} & +1 \cdot \binom{n+5}{6} \\ \vdots & &\vdots \end{array}$$ This looks a bit more regular than the expressions by Bernoulli-polynomials and possibly the relations between the sums-expressions of odd orders $s_{2k+1}(n)$ can be written using these polynomials nicely and with more ease than with the Bernoulli-polynomials, don't know.

The first of your equality would then be expressed as $$\begin{array} {rcl} \sum_{k=1}^n k^3 &=&\left(\sum_{k=1}^n k \right)^2 \\ s_3(n) &=& s_1(n)^2 \\ 1 \cdot \binom{n+1}{4} +4 \cdot \binom{n+2}{4} +1 \cdot \binom{n+3}{4} &=& \left( 1 \cdot \binom{n+1}{2}\right)^2 \\ \end{array}$$ and which must then be expanded to reproduce the equality (but I've not the time at the moment to do this in the required length...)