Introduction and main results

In this manuscript we derive formulas for multifold sums of powers by utilizing Newton’s interpolation formula. Furthermore, we provide formulas for multifold sums of powers in terms of Stirling numbers of the second kind and Eulerian numbers.

The idea to derive sums of powers using difference operator and Newton’s series is quite generic, thus, formulas for sums of powers using backward and central differences can be found in the works (Kolosov 2026, 2025b).

Allow us to start from the definition of multifold sums of powers. We utilize the recurrence proposed by Donald Knuth in his article Johann Faulhaber and sums of powers, see (Knuth, Donald E. 1993) \[\begin{align*} \Sigma^{0}\,{n}^{m} &= n^m \\ \Sigma^{1}\,{n}^{m} &= \Sigma^{0}\,{1}^{m} + \Sigma^{0}\,{2}^{m} + \cdots + \Sigma^{0}\,{n}^{m} \\ \Sigma^{r+1}\,{n}^{m} &= \Sigma^{r}\,{1}^{m} + \Sigma^{r}\,{2}^{m} + \cdots + \Sigma^{r}\,{n}^{m} \end{align*}\] Throughout the paper, we utilize the Newton’s interpolation formula as stated below

Which is indeed true, because \[\begin{align*} n^3 &= 0 \binom{n}{0} + 1 \binom{n}{1} + 6 \binom{n}{2} + 6 \binom{n}{3} \\ n^3 &= 1 \binom{n-1}{0} + 7 \binom{n-1}{1} + 12 \binom{n-1}{2} + 6 \binom{n-1}{3} \\ n^3 &= 8 \binom{n-2}{0} + 19 \binom{n-2}{1} + 18 \binom{n-2}{2} + 6 \binom{n-2}{3} \end{align*}\]

Thus, for an arbitrary integer \(t\), the ordinary sum of powers is \[\begin{align*} \Sigma^{1}\,{n}^{m} = \sum_{k=1}^{n} \sum_{j=0}^{m} \binom{-t+k}{j} \Delta^{j} t^m = \sum_{j=0}^{m} \Delta^{j} t^m \sum_{k=1}^{n} \binom{-t+k}{j} \end{align*}\]

Hence, \[\begin{align*} \sum_{k=1}^{n} \binom{-t+k}{j} &= (-1)^j \binom{j+t}{j+1} + \binom{n-t+1}{j+1} - \binom{-t}{j} \\ &= (-1)^j \binom{j+t}{j+1} - (-1)^j \binom{j+t-1}{j} + \binom{n-t+1}{j+1} \\ &= (-1)^j \binom{j+t-1}{j+1} + \binom{n-t+1}{j+1} \end{align*}\] because \(\binom{-t}{j} = (-1)^j \binom{j+t-1}{j}\) and by recurrence \(\binom{t+1}{k} = \binom{t}{k} + \binom{t}{k+1}\).

This implies formula for sums of powers

The special cases for \(t=0\) and \(t=1\) are widely known and appear in literature quite frequently. For \(t=0\) and \(m=3\) we have the famous identity \[\begin{align*} \Sigma^{1}\,{n}^{3} = 0 \binom{n+1}{1} + 1 \binom{n+1}{2} + 6 \binom{n+1}{3} + 6 \binom{n+1}{4} \end{align*}\] which was discussed in (Graham, Ronald L. and Knuth, Donald E. and Patashnik, Oren 1994, 190) and in (Pfaff 2007).

The coefficients \(0, 1, 6, 6, 0, 1, 14, 36, 24, \dots\) is the sequence A131689 in the OEIS (Sloane, Neil J.A. and others 2003).

The special cases for \(t=1\) and \(m=2,3,4,5\) were discussed in (Cereceda, José L. 2022). For instance, \[\begin{align*} \Sigma^{1}\,{n}^{3} &= 1 \binom{n}{1} + 7 \binom{n}{2} + 12 \binom{n}{3} + 6 \binom{n}{4} \\ \Sigma^{1}\,{n}^{4} &= 1 \binom{n}{1} +15 \binom{n}{2} +50 \binom{n}{3} +60 \binom{n}{4} +24 \binom{n}{5} \end{align*}\]

The coefficients \(1,7,12,6,1,15, \dots\) is the sequence A028246 in the OEIS (Sloane, Neil J.A. and others 2003).

Interestingly enough that the paper (Cereceda, José L. 2022) gives the formula for sums of powers \[\begin{align*} \Sigma^{1}\,{n}^{k} = \sum_{j=0}^{k} j! \left[ \binom{n+1-r}{j+1} + (-1)^j \binom{r+j-1}{j+1} \right] \genfrac{\{}{\}}{0pt}{}{k}{j}_{r} \end{align*}\] where \(\genfrac{\{}{\}}{0pt}{}{k}{j}_{r}\) are generalized Stirling numbers of the second kind. The formula above is identical to the proposition [prop:ordinary-sums-of-powers-via-newtons-series], which implies that finite differences can be expressed in terms of generalized Stirling numbers of the second kind, that is \(\Delta^{j} t^m = j! \genfrac{\{}{\}}{0pt}{}{m}{j}_{t}\).

By considering the special cases of the proposition [prop:ordinary-sums-of-powers-via-newtons-series] for \(t=4\), we observe rather unexpected formulas for sums of powers, namely \[\begin{align*} \Sigma^{1}\,{n}^{0} &= 1 \left( \binom{n-3}{1} + \binom{3}{1} \right) \\ \Sigma^{1}\,{n}^{1} &= 4 \left( \binom{n-3}{1} + \binom{3}{1} \right) + 1 \left( \binom{n-3}{2} - \binom{4}{2} \right) \\ \Sigma^{1}\,{n}^{2} &= 16 \left( \binom{n-3}{1} + \binom{3}{1} \right) + 9 \left( \binom{n-3}{2} - \binom{4}{2} \right) + 2 \left( \binom{n-2}{3} + \binom{5}{3} \right) \\ \Sigma^{1}\,{n}^{3} &= 64 \left( \binom{n-3}{1} + \binom{3}{1} \right) + 61 \left( \binom{n-3}{2} - \binom{4}{2} \right) + 30 \left( \binom{n-3}{3} + \binom{5}{3} \right) \\ &+ 6 \left( \binom{n-3}{4} - \binom{6}{4} \right) \end{align*}\] The coefficients \(1,4,1,16,9,\dots\) is the sequence A391633 in the OEIS (Sloane, Neil J.A. and others 2003). In general, \[\begin{align*} \Sigma^{1}\,{n}^{m} = \sum_{j=0}^{m} \Delta^{j} 4^m \left[ \binom{n-3}{j+1} + (-1)^j \binom{j+3}{j+1} \right]. \end{align*}\]

To derive formula for double sum of powers, we apply summation operator over the ordinary sum of powers [prop:ordinary-sums-of-powers-via-newtons-series] again, thus \[\begin{align*} \Sigma^{2}\,{n}^{m} = \sum_{j=0}^{m} \Delta^{j} t^{m} \left[ (-1)^{j} \sum_{k=1}^{n} \binom{j+t-1}{j+1} + \sum_{k=1}^{n} \binom{k-t+1}{j+1} \right] \end{align*}\] which yields \[\begin{align*} \Sigma^{2}\,{n}^{m} = \sum_{j=0}^{m} \Delta^{j} t^{m} \left[ (-1)^{j} \binom{j+t-1}{j+1} n + \sum_{k=1}^{n} \binom{k-t+1}{j+1} \right] \end{align*}\] Thus,

For example, by setting \(t=5\), the double sums of powers are \[\begin{align*} \Sigma^{2}\,{n}^{0} &= 1 \left( \binom{n-3}{2} + \binom{4}{1} n - \binom{4}{2} \right) \\ \Sigma^{2}\,{n}^{1} &= 5 \left( \binom{n-3}{2} + \binom{4}{1} n - \binom{4}{2} \right) + 1 \left( \binom{n-3}{3} - \binom{5}{2} n + \binom{5}{3} \right) \\ \Sigma^{2}\,{n}^{2} &= 25 \left( \binom{n-3}{2} + \binom{4}{1} n - \binom{4}{2} \right) + 11 \left( \binom{n-3}{3} - \binom{5}{2} n + \binom{5}{3} \right) \\ &+ 2 \left( \binom{n-3}{4} + \binom{6}{3} n - \binom{6}{4} \right) \\ \Sigma^{2}\,{n}^{3} &= 125 \left( \binom{n-3}{2} + \binom{4}{1} n - \binom{4}{2} \right) + 91 \left( \binom{n-3}{3} - \binom{5}{2} n + \binom{5}{3} \right) \\ &+ 36 \left( \binom{n-3}{4} + \binom{6}{3} n - \binom{6}{4} \right) + 6 \left( \binom{n-3}{5} - \binom{7}{4} n + \binom{7}{5} \right) \end{align*}\] The coefficients \(1,5,1,25,11,2,\dots\) is the sequence A391635 in the OEIS (Sloane, Neil J.A. and others 2003).

In general, \[\begin{align*} \Sigma^{2}\,{n}^{m} = \sum_{j=0}^{m} \Delta^{j} 5^{m} \left[ \binom{n-3}{j+2} + (-1)^{j} \binom{j+4}{j+1} n^1 + (-1)^{j+1} \binom{j+4}{j+2} n^0 \right] \end{align*}\]

Similarly, we obtain the formula for the triple sums of powers

For example, by setting \(t=4\), the triple sums of powers are \[\begin{align*} \Sigma^{3}\,{n}^{0} &= 1 \left( \binom{n-1}{3} + \binom{3}{1} \Sigma^{2}\,{n}^{0} - \binom{3}{2} \Sigma^{1}\,{n}^{0} + \binom{3}{3} \Sigma^{0}\,{n}^{0} \right) \\ \Sigma^{3}\,{n}^{1} &= 4 \left( \binom{n-1}{3} + \binom{3}{1} \Sigma^{2}\,{n}^{0} - \binom{3}{2} \Sigma^{1}\,{n}^{0} + \binom{3}{3} \Sigma^{0}\,{n}^{0} \right) \\ &+ 1 \left( \binom{n-1}{4} - \binom{4}{2} \Sigma^{2}\,{n}^{0} + \binom{4}{3} \Sigma^{1}\,{n}^{0} - \binom{4}{4} \Sigma^{0}\,{n}^{0} \right) \\ \Sigma^{3}\,{n}^{2} &= 16 \left( \binom{n-1}{3} + \binom{3}{1} \Sigma^{2}\,{n}^{0} - \binom{3}{2} \Sigma^{1}\,{n}^{0} + \binom{3}{3} \Sigma^{0}\,{n}^{0} \right) \\ &+ 9 \left( \binom{n-1}{4} - \binom{4}{2} \Sigma^{2}\,{n}^{0} + \binom{4}{3} \Sigma^{1}\,{n}^{0} - \binom{4}{4} \Sigma^{0}\,{n}^{0} \right) \\ &+ 2 \left( \binom{n-1}{5} + \binom{5}{3} \Sigma^{2}\,{n}^{0} - \binom{5}{4} \Sigma^{1}\,{n}^{0} + \binom{5}{5} \Sigma^{0}\,{n}^{0} \right) \end{align*}\] In general, \[\begin{align*} \Sigma^{3}\,{n}^{m} &= \sum_{j=0}^{m} \Delta^{j} 4^m \Bigg[ (-1)^{j} \binom{j+3}{j+1} \Sigma^{2}\,{n}^{0} + (-1)^{j+1} \binom{j+3}{j+2} \Sigma^{1}\,{n}^{0} + \\ &+ (-1)^{j+2} \binom{j+3}{j+3} \Sigma^{0}\,{n}^{0} + \binom{n-1}{j+3} \Bigg] \end{align*}\]

Continuing similarly, we derive formula for multifold sums of powers, which is

In its explicit form

We may observe that

Which yields the following binomial variations of the multifold sums of powers [theorem:multifold-sums-of-powers-via-newtons-series]

Finite differences of powers are closely related to Stirling numbers of the second kind

Which implies variations of the formulas for sums of powers

By setting \(t=0\) into [prop:ordinary-sums-of-powers-via-stirling-numbers] yields well-known identity for sums of powers in terms of Stirling numbers

In general,

The proposition above can be presented in a pure binomial form as well, by means of the identity [prop:multifold-sum-of-zero-powers]: \(\Sigma^{r}\,{n}^{0} = \binom{r+n-1}{r}\).

In addition, we are able to express multifold sums of powers via Eulerian numbers, by expressing the forward finite difference via the Worpitzky identity (Worpitzky 1883)

where \(\genfrac{\langle}{\rangle}{0pt}{}{n}{k}\) are Eulerian numbers. Thus,

Therefore,

Which implies variations of the formulas for sums of powers

By setting \(t=0\) into [prop:sums-of-powers-via-eulerian-numbers] yields a well-known identity for sums of powers in terms of Eulerian numbers

Backward difference form

The formula for multifold sums of powers via Newton’s series [theorem:multifold-sums-of-powers-via-newtons-series] can be expressed in terms of backward differences easily, because

Thus,

Central difference form

The formula for multifold sums of powers via Newton’s series [theorem:multifold-sums-of-powers-via-newtons-series] can be expressed in terms of central differences easily, because

Thus,

Future research

In this manuscript we focus on the idea to combine Newton’s interpolation formula for forward differences [prop:newton-series-for-monomial-reversed] and the hockey-stick family identities (e.g [prop:segmented-hockey-stick-identity]) to express the sums of powers seamlessly.

This particular idea is great, however it can be generalized even further. Thus, the main idea is to utilize an interpolation formula for power \(n^m\) in terms of abstract difference operator \(D(n^m)\) and binomial coefficients \(\binom{f(n)}{k}\).

We are not limited to linear difference operators \(\Delta, \nabla, \delta\), there is a complete theory on dynamic differential equations on time-scales, founded by Bohner and Peterson (Bohner, Martin and Peterson, Allan 2001). For example, Jackson (Jackson, Frederick H. 1909) derivative is given by

\[\begin{align*} D_q(f(x)) = \frac{f(qx) - f(x)}{qx-x} \end{align*}\]

Thus, finding interpolation formula for \(f(x)=x^n\) in terms of \(D_q(f(x))\) implies new formulas for sums of powers.

In general, let be an abstract interpolation formula for powers \[\begin{align*} n^m = \sum_{k} \binom{f(n)}{k} D(n^m, k) \end{align*}\] Thus, the formula of sums of powers involves the abstract difference operator \(D\) evaluated at some point \(k\) and the hockey-stick family identity over the binomial coefficients \(\binom{n}{k}\) \[\begin{align*} \Sigma^{1}\,{n}^{m} = \sum_{k} D(n^m, k) \sum_{j \leq n} \binom{f(j)}{k} \end{align*}\] Similarly, for multifold sums of powers \[\begin{align*} \Sigma^{r}\,{n}^{m} = \sum_{k} D(n^m, k) \binom{f(n+r)}{k} \end{align*}\] Many interpolation approaches involve rising factorials \(x_{(n)}\), falling factorials \((x)_n\), or regular factorials \(n!\), and thus can be expressed in terms of binomial coefficients, because \[\begin{align*} \frac{(x)_n}{n!} = \binom{x}{n}; \quad \frac{x_{(n)}}{n!} = \binom{x+n-1}{n}. \end{align*}\]

In particular, Donald Knuth provides the formula for multifold sums of odd powers (Knuth, Donald E. 1993) based on operator of central finite differences evaluated in zero, and Newton’s interpolation formula for central differences (Kolosov 2025b)

where \(T(n,k)\) are the central factorial numbers of the second kind, see (Steffensen, Johan Frederik 1927, sec. 58) and (Carlitz and Riordan 1963, formula (10a)) \[\begin{align*} T(n, k) = \frac{1}{k!} \delta^k 0^n = \frac{1}{k!} \sum_{j=0}^{k} (-1)^j \binom{k}{j} \left( \frac{k}{2} - j \right)^n \end{align*}\] Central factorial numbers of the second kind \(T(n,k)\) were defined by Riordan in (Riordan 1968, ch. 6.5, formula (24)) by polynomial identity

The sequence A008957 in the OEIS (Sloane, Neil J.A. and others 2003) provides non-zero central factorial numbers of the second kind \(T(2n,2k)\).

As future research direction, the Knuth’s formula [prop:knuth-sums-of-odd-powers] utilizes the operator of central finite differences evaluated in zero. Thus, it is worth to investigate the existence of sums of odd powers involving the central differences evaluated at an arbitrary integer point \(t\), similar to multifold sums of powers via Newton’s series for central differences.

Conclusions

In this manuscript we have derived formulas for multifold sums of powers by utilizing Newton’s interpolation formula. Furthermore, we provide formulas for multifold sums of powers in terms of Stirling numbers of the second kind, and Eulerian numbers. In addition, in [sec:future-research] we discuss the future research directions that may lead to a complete framework for sums of powers, by means of combining interpolation approaches, binomial coefficients and variations of the hockey-stick identity for binomial coefficients. The most important results of this manuscript are validated using Mathematica programs, see [sec:mathematica-programs].

Acknowledgements

The author is grateful to Markus Scheuer for his valuable contribution of the list of references (Markus Scheuer 2020) about the fact that central factorial numbers of the second kind arise from central differences of nothing.

Metadata

• Version: 1.5.0+main.a49dfab

• License: This work is licensed under a CC BY 4.0 License.

• Sources: github.com/kolosovpetro/NewtonsInterpolationFormulaAndSumsOfPowers

• ORCID: 0000-0002-6544-8880

• Email: kolosovp94@gmail.com

• MSC2010: 05A19, 05A10, 41A15, 11B83

• Keywords: Sums of powers, Newton’s interpolation formula, Finite differences, Binomial coefficients, Newton’s series, Faulhaber’s formula, Bernoulli numbers, Bernoulli polynomials, Interpolation, Combinatorics, Central factorial numbers, Stirling numbers, Eulerian numbers, Worpitzky identity, OEIS

Proof of Segmented hockey stick identity

First we split the sum \(\sum_{k=0}^{n} \binom{-t+k}{j}\) into two sub-sums so that we discuss them separately \[\begin{align*} \sum_{k=0}^{n} \binom{-t+k}{j} &= \sum_{k=0}^{t-1} \binom{-t+k}{j} + \sum_{k=t}^{n} \binom{-t+k}{j} \end{align*}\] We assume that the two sums above run over the partition \(\{0,1,2,\cdots,t,\cdots, n\}\), with \(t<n\). Considering the sum \(\sum_{k=0}^{t-1} \binom{-t+k}{j}\) we notice that \[\begin{align*} \sum_{k=0}^{t-1} \binom{-t+k}{j} = \binom{-t}{j} &+ \binom{-t +1}{j} + \binom{-t+2}{j} + \cdots + \\ &+ \binom{-t+t-2}{j} + \binom{-t+t-1}{j} \end{align*}\] Thus \[\begin{align*} \sum_{k=0}^{t-1} \binom{-t+k}{j} &= \sum_{k=1}^{t} \binom{-k}{j} = \sum_{k=0}^{t-1} \binom{-k-1}{j} \end{align*}\] By using the identity \(\binom{-k}{j} = (-1)^j \binom{j+k-1}{j}\), we obtain \[\begin{align*} \binom{-k-1}{j} &= \binom{-(k+1)}{j} = (-1)^j \binom{j+k}{j} \end{align*}\] Thus \[\begin{align*} \sum_{k=0}^{t-1} \binom{-t+k}{j} = (-1)^j \sum_{k=0}^{t-1} \binom{j+k}{j} = (-1)^{j} \binom{j+t}{j+1} \end{align*}\] by the hockey-stick identity \(\sum_{k=0}^{t} \binom{j+k}{j} = \binom{j+t+1}{j+1}\).

Considering the sum \(\sum_{k=t}^{n} \binom{-t+k}{j}\) we notice that \[\begin{align*} \sum_{k=t}^{n} \binom{-t+k}{j} = \sum_{k=0}^{n-t} \binom{k}{j} \end{align*}\] Thus \[\begin{align*} \sum_{k=t}^{n} \binom{-t+k}{j} = \sum_{k=0}^{n-t} \binom{k}{j} = \binom{n-t+1}{j+1} \end{align*}\] again by the hockey-stick identity \(\sum_{k=0}^{t} \binom{j+k}{j} = \binom{j+t+1}{j+1}\). Combining the two parts, we obtain \[\begin{align*} \sum_{k=0}^{n} \binom{-t+k}{j} &= (-1)^j \binom{j+t}{j+1} + \binom{n-t+1}{j+1} \end{align*}\] This completes the proof.

Mathematica programs

Use the Mathematica package (Kolosov 2025a) to validate the results