# Polynomials calculating sums of powers of arithmetic progressions

The polynomials calculating sums of powers of arithmetic progressions are polynomials in a variable that depend both on the particular arithmetic progression constituting the basis of the summed powers and on the constant exponent, non-negative integer, chosen. Their degree always exceeds the constant exponent by one unit and have the property that when the polynomial variable coincides with the number of summed addends, the result of the polynomial function also coincides with that of the sum.

The problem therefore consists in finding ${\displaystyle S_{h,d}^{m}(n)}$i.e. polynomials as a function of ${\displaystyle n}$ calculating sums of ${\displaystyle n}$ addends:

${\displaystyle \sum _{k=0}^{n-1}(h+kd)^{m}=h^{m}+(h+d)^{m}+\cdots +(h+(n-1)d)^{m},}$

with ${\displaystyle m}$[1] and ${\displaystyle n}$ integers positive, ${\displaystyle h}$ first term of an arithmetic progression and ${\displaystyle d\neq 0}$ the common difference. The two parameters can be not only integers but also rational, real and even complex.

## History

### Ancient period

The history of the problem begins in antiquity and coincides with that of some of its special cases. The case ${\displaystyle m=1,}$ coincides with that of the calculation of the arithmetic series, the sum of the first ${\displaystyle n}$ values of an arithmetic progression. This problem is quite simple but the case already known by the Pythagorean school for its connection with triangular numbers is historically interesting:

${\displaystyle 1+2+\dots +n={\frac {1}{2}}n^{2}+{\frac {1}{2}}n,}$   Polynomial ${\displaystyle S_{1,1}^{1}(n)}$ calculating the sum of the first ${\displaystyle n}$ natural numbers.

For ${\displaystyle m>1,}$ the first cases encountered in the history of mathematics are:

${\displaystyle 1+3+\dots +2n-1=n^{2},}$   Polynomial ${\displaystyle S_{1,2}^{1}(n)}$ calculating the sum of the first ${\displaystyle n}$ successive odds forming a square. A property probably well known by the Pythagoreans themselves who, in constructing their figured numbers, had to add each time a gnomon consisting of an odd number of points to obtain the next perfect square.
${\displaystyle 1^{2}+2^{2}+\ldots +n^{2}={\frac {1}{3}}n^{3}+{\frac {1}{2}}n^{2}+{\frac {1}{6}}n,}$   Polynomial ${\displaystyle S_{1,1}^{2}(n)}$ calculating the sum of the squares of the successive integers. Property that we find demonstrated in Spirals, a work of Archimedes;[2]
${\displaystyle 1^{3}+2^{3}+\ldots +n^{3}={\frac {1}{4}}n^{4}+{\frac {1}{2}}n^{3}+{\frac {1}{4}}n^{2},}$   Polynomial ${\displaystyle S_{1,1}^{3}(n)}$ calculating the sum of the cubes of the successive integers. Corollary of a theorem of Nicomachus of Gerasa...[2]

L'insieme ${\displaystyle S_{1,1}^{m}(n)}$ of the cases, to which the two preceding polynomials belong, constitutes the classical problem of powers of successive integers.

### Middle period

Over time, many other mathematicians became interested in the problem and made various contributions to its solution. These include Aryabhata, Al-Karaji, Ibn al-Haytham, Thomas Harriot, Johann Faulhaber, Pierre de Fermat and Blaise Pascal who recursively solved the problem of the sum of powers of successive integers by considering an identity that allowed to obtain a polynomial of degree ${\displaystyle m+1}$ already knowing the previous ones.[2]

In 1713 the family of Jacob Bernoulli posthumously publishes his Artis Conjectandi [3] where the first 10 polynomials of this infinite series appear together with a general formula dependent on particular numbers that were soon named after him. The formula was instead attributed to Johann Faulhaber[4] for his worthy contributions recognized by Bernoulli himself. It was also immediately clear that the polynomials ${\displaystyle S_{0,1}^{m}(n)}$ calcolating the sum of ${\displaystyle n}$ powers of successive integers starting from zero were very similar to those starting from one. This is because it is evident that ${\displaystyle S_{1,1}^{m}(n)-S_{0,1}^{m}(n)=n^{m}}$and that therefore polynomials of degree ${\displaystyle m+1}$ of the form ${\displaystyle {\frac {1}{m+1}}n^{m+1}+{\frac {1}{2}}n^{m}+\cdots }$[3] subtracted the monomial difference ${\displaystyle n^{m}}$ they become ${\displaystyle {\frac {1}{m+1}}n^{m+1}-{\frac {1}{2}}n^{m}+\cdots }$.

However, a proof of Faulhaber's formula was missing, which was given more than a century later by Carl G. Jacobi[5] who benefited from the progress of mathematical analysis using the development in infinite series of an exponential function generating Bernoulli numbers.

### Modern period

In 1982 A.W.F. Edwards publishes an article [6] in which he shows that Pascal's identity can be expressed by means of triangular matrices containing the Pascal's triangle deprived of 'last element of each line:

${\displaystyle {\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\n^{5}\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0&0\\1&2&0&0&0\\1&3&3&0&0\\1&4&6&4&0\\1&5&10&10&5\end{pmatrix}}{\begin{pmatrix}n\\\sum _{k=0}^{n-1}k^{1}\\\sum _{k=0}^{n-1}k^{2}\\\sum _{k=0}^{n-1}k^{3}\\\sum _{k=0}^{n-1}k^{4}\\\end{pmatrix}}}$[7][8]

The example is limited by the choice of a fifth order matrix but is easily extendable to higher orders. The equation can be written as: ${\displaystyle {\vec {N}}=A{\vec {S}}}$ and multiplying the two sides of the equation to the left by ${\displaystyle A^{-1}}$ , inverse of the matrix A, we obtain ${\displaystyle A^{-1}{\vec {N}}={\vec {S}}}$ which allows to arrive directly at the polynomial coefficients without directly using the Bernoulli numbers. Other authors after Edwards dealing with various aspects of the power sum problem take the matrix path [9] and studying aspects of the problem in their articles useful tools such as the Vandermonde vector.[10] Other researchers continue to explore through the traditional analytic route [11] and generalize the problem of the sum of successive integers to any geometric progression[12][13]

The coefficients of the polynomials ${\displaystyle S_{h,d}^{m}}$ are found through recursive formulas and in other ways that are interesting for number theory as the expression of the result of the sum as a function of Bernoulli polynomials or the formulas involving the Stirling numbers and the r-Whitney numbers of the first and second kind [14] Finally, Edwards' matrix approach was also generalized to any arithmetic progressions [15]

## Solution by matrix method

The general problem has recently been solved [15] through the use of binomial matrices easily constructible knowing the binomial coefficients and the Pascal's triangle. It is shown that, having chosen the parameters ${\displaystyle h}$ and ${\displaystyle d}$ which determine the arithmetic progression and a positive integer ${\displaystyle m,}$ we find ${\displaystyle m+1}$ polynomials corresponding to the following sums of powers:

${\displaystyle S_{h,d}^{r-1}=c_{1}n+c_{2}n^{2}+\ldots +c_{r}n^{r},\qquad {\text{with }}r=1,2,\ldots ,m+1,}$
with the polynomial coefficients elements of the row ${\displaystyle r}$ of the triangular matrix ${\displaystyle G(h,d)=T(h,d)A^{-1}}$ of order ${\displaystyle m+1}$.

Here is the solving formula in the particular case ${\displaystyle m=3}$ which gives the polynomials of a given arithmetic progression with exponents from 0 to 3:

${\displaystyle {\begin{pmatrix}S_{h,d}^{0}(n)\\S_{h,d}^{1}(n)\\S_{h,d}^{2}(n)\\S_{h,d}^{3}(n)\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\h&d&0&0\\h^{2}&2hd&d^{2}&0\\h^{3}&3h^{2}d&3hd^{2}&d^{3}\end{pmatrix}}{\begin{pmatrix}1&0&0&0\\1&2&0&0\\1&3&3&0\\1&4&6&4\\\end{pmatrix}}^{-1}{\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\\end{pmatrix}}}$

The equation that can be easily extended to different values of m (non-negative integers) is summarized and generalized as follows:

${\displaystyle {\vec {S}}_{h,d}(n)=T(h,d)A^{-1}n{\vec {V}}(n)}$ or also by placing with ${\displaystyle G(h,d)=T(h,d)A^{-1}}$

${\displaystyle {\vec {S}}_{h,d}(n)=G(h,d)n{\vec {V}}(n)}$ [15]

Here is the rigorous definition of the matrices and the Vandermonde vector:

{\displaystyle {\begin{aligned}{[A]_{r,c}}&={\begin{cases}0,&{\text{if }}c>r,\\{\dbinom {r}{c-1}},&{\text{if }}c\leq r,\end{cases}}\\[1ex][T(h,d)]_{r,c}&={\begin{cases}0,&{\text{if }}c>r,\\{\dbinom {r-1}{c-1}}h^{r-c}d^{c-1}&{\text{if }}c\leq r.\end{cases}}\\[1ex][{\vec {V}}(n)]_{r}&={\begin{cases}1,&{\text{if }}r=1,\\n^{r-1}&{\text{if }}r>1.\end{cases}}\end{aligned}}}
for ${\displaystyle m=3}$ it results therefore
${\displaystyle n{\vec {V}}(n)=n{\begin{pmatrix}1\\n\\n^{2}\\n^{3}\\\end{pmatrix}}={\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\\end{pmatrix}}\quad {\text{and also:}}\quad {\vec {S}}_{h,d}(n)={\begin{pmatrix}S_{h,d}^{0}(n)\\S_{h,d}^{1}(n)\\S_{h,d}^{2}(n)\\S_{h,d}^{3}(n)\\\end{pmatrix}}=\sum _{k=0}^{n-1}{\vec {V}}(h+kd)}$

Matrix A is that of Edwards [8] already seen, a lower triangular matrix that reproduces, in the non-null elements, the triangle of Pascal deprived of the last element of each row. The elements of ${\displaystyle T(h,d)}$ on the other hand are the monomials of the power development ${\displaystyle (h+d)^{r-1},}$ for ${\displaystyle r=1,\ldots ,m+1.}$.

${\displaystyle T(0,1)}$ is the neutral element of the row by column product so that the general equation in this case becomes:

${\displaystyle {\vec {S}}_{0,1}(n)=A^{-1}n{\vec {V}}(n)}$ that is the one discovered by Edwards [8]

To arrive from this particular case to prove the general one, it is sufficient to multiply on the left the two members of the equation by the matrix ${\displaystyle T(h,d)}$ after having ascertained the following identity ${\displaystyle T(h,d){\vec {V}}(n)=V(h+dn)}$ [15]

### Sum of powers of successive odd numbers

We use the previous formula to solve the problem of adding powers of successive odds:.[16] The odds correspond to the arithmetic progression with the first element ${\displaystyle h=1}$ and as reason ${\displaystyle d=2.}$ We set m = 4 to find the first five polynomials calculating sums of powers of odd. Calculated ${\displaystyle T(1,2)}$ we obtain:

${\displaystyle T(1,2)={\begin{pmatrix}1&0&0&0&0\\1&2&0&0&0\\1&4&4&0&0\\1&6&12&8&0\\1&8&24&32&16\\\end{pmatrix}}\quad A={\begin{pmatrix}1&0&0&0&0\\1&2&0&0&0\\1&3&3&0&0\\1&4&6&4&0\\1&5&10&10&5\\\end{pmatrix}}\quad A^{-1}={\begin{pmatrix}1&0&0&0&0\\-{\frac {1}{2}}&{\frac {1}{2}}&0&0&0\\{\frac {1}{6}}&-{\frac {1}{2}}&{\frac {1}{3}}&0&0\\0&{\frac {1}{4}}&-{\frac {1}{2}}&{\frac {1}{4}}&0\\-{\frac {1}{30}}&0&{\frac {1}{3}}&-{\frac {1}{2}}&{\frac {1}{5}}\\\end{pmatrix}}}$

We have therefore

${\displaystyle G(1,2)=T(1,2)A^{-1}={\begin{pmatrix}1&0&0&0&0\\0&1&0&0&0\\-{\frac {1}{3}}&0&{\frac {4}{3}}&0&0\\0&-1&0&2&0\\{\frac {7}{15}}&0&-{\frac {8}{3}}&0&{\frac {16}{5}}\end{pmatrix}}.}$

At this point the general equation ${\displaystyle {\vec {S}}_{1,2}(n)=G(1,2)n{\vec {V}}(n)}$ for ${\displaystyle m=4}$ and the damage done product:

${\displaystyle {\begin{pmatrix}\sum _{k=0}^{n-1}(1+2k)^{0}\\\sum _{k=0}^{n-1}(1+2k)^{1}\\\sum _{k=0}^{n-1}(1+2k)^{2}\\\sum _{k=0}^{n-1}(1+2k)^{3}\\\sum _{k=0}^{n-1}(1+2k)^{4}\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0&0\\0&1&0&0&0\\-{\frac {1}{3}}&0&{\frac {4}{3}}&0&0\\0&-1&0&2&0\\{\frac {7}{15}}&0&-{\frac {8}{3}}&0&{\frac {16}{5}}\end{pmatrix}}{\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\n^{5}\\\end{pmatrix}}={\begin{pmatrix}n\\n^{2}\\-{\frac {1}{3}}n+{\frac {4}{3}}n^{3}\\-n^{2}+2n^{4}\\{\frac {7}{15}}n-{\frac {8}{3}}n^{3}+{\frac {16}{5}}n^{5}\\\end{pmatrix}}.}$
using the last line (${\displaystyle r=5}$) we get then
${\displaystyle {S}_{1,2}^{4}={\frac {7}{15}}n-{\frac {8}{3}}n^{3}+{\frac {16}{5}}n^{5}}$
and using the other rows:
${\displaystyle {S}_{1,2}^{3}=-n^{2}+2n^{4};\quad {S}_{1,2}^{2}=-{\frac {1}{3}}n+{\frac {4}{3}}n^{3};\quad {S}_{1,2}^{1}=n^{2};\quad {S}_{1,2}^{0}=n.}$

### Sum of successive integers starting with 1

Chosen ${\displaystyle m=3}$ and calculated ${\displaystyle A^{-1}}$ and T(1,1) which corresponds to Pascal's triangle:

${\displaystyle {\begin{pmatrix}S_{1,1}^{0}(n)\\S_{1,1}^{1}(n)\\S_{1,1}^{2}(n)\\S_{1,1}^{3}(n)\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\+{\frac {1}{2}}&{\frac {1}{2}}&0&0\\{\frac {1}{6}}&+{\frac {1}{2}}&{\frac {1}{3}}&0\\0&{\frac {1}{4}}&+{\frac {1}{2}}&{\frac {1}{4}}\end{pmatrix}}{\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\\end{pmatrix}}={\begin{pmatrix}n\\+{\frac {1}{2}}n+{\frac {1}{2}}n^{2}\\{\frac {1}{6}}n+{\frac {1}{2}}n^{2}+{\frac {1}{3}}n^{3}\\{\frac {1}{4}}n^{2}+{\frac {1}{2}}n^{3}+{\frac {1}{4}}n^{4}\end{pmatrix}}}$

### Sum of successive integers starting with 0

Chosen ${\displaystyle m=3}$ and calculated ${\displaystyle A^{-1}}$ and T(0,1) unit matrix:

${\displaystyle {\begin{pmatrix}S_{0,1}^{0}(n)\\S_{0,1}^{1}(n)\\S_{0,1}^{2}(n)\\S_{0,1}^{3}(n)\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\-{\frac {1}{2}}&{\frac {1}{2}}&0&0\\{\frac {1}{6}}&-{\frac {1}{2}}&{\frac {1}{3}}&0\\0&{\frac {1}{4}}&-{\frac {1}{2}}&{\frac {1}{4}}\end{pmatrix}}{\begin{pmatrix}n\\n^{2}\\n^{3}\\n^{4}\\\end{pmatrix}}={\begin{pmatrix}n\\-{\frac {1}{2}}n+{\frac {1}{2}}n^{2}\\{\frac {1}{6}}n-{\frac {1}{2}}n^{2}+{\frac {1}{3}}n^{3}\\{\frac {1}{4}}n^{2}-{\frac {1}{2}}n^{3}+{\frac {1}{4}}n^{4}\end{pmatrix}}}$

### Progression -1,3,7,11,15 ...

Chosen again ${\displaystyle m=3}$, calculated ${\displaystyle T(-1,4)}$, exploited the result of the previous paragraph and the associative property:

${\displaystyle {\begin{pmatrix}S_{-1,4}^{0}(n)\\S_{-1,4}^{1}(n)\\S_{-1,4}^{2}(n)\\S_{-1,4}^{3}(n)\\\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\-1&4&0&0\\1&-8&16&0\\-1&12&-48&64\end{pmatrix}}{\begin{pmatrix}n\\-{\frac {1}{2}}n+{\frac {1}{2}}n^{2}\\{\frac {1}{6}}n-{\frac {1}{2}}n^{2}+{\frac {1}{3}}n^{3}\\{\frac {1}{4}}n^{2}-{\frac {1}{2}}n^{3}+{\frac {1}{4}}n^{4}\end{pmatrix}}={\begin{pmatrix}n\\-3n+2n^{2}\\{\frac {23}{3}}n-12n^{2}+{\frac {16}{3}}n^{3}\\-15n+46n^{2}-48n^{3}+64n^{4}\end{pmatrix}}}$

## Generalization of Faulhaber's formula

The matrix ${\displaystyle G(h,d)}$ can be expressed as a function of the Bernoulli polynomials in the following way[17][18]

${\displaystyle [G(h,d)]_{r,c}={\begin{cases}0,&{\text{if }}c>r,\\{\frac {d^{r-1}}{r}}{\binom {r}{c}}B_{r-c}\left({\frac {h}{d}}\right),&{\text{if }}c\leq r,\end{cases}}}$

which for ${\displaystyle m=5}$ becomes

${\displaystyle G(h,d)={\begin{pmatrix}{\frac {1}{1}}{\binom {1}{1}}B_{0}({\frac {h}{d}})&0&0&0&0&0\\{\frac {d}{2}}{\binom {2}{1}}B_{1}({\frac {h}{d}})&{\frac {d}{3}}{\binom {2}{2}}B_{0}({\frac {h}{d}})&0&0&0&0\\{\frac {d^{2}}{3}}{\binom {3}{1}}B_{2}({\frac {h}{d}})&{\frac {d^{2}}{3}}{\binom {3}{2}}B_{1}({\frac {h}{d}})&{\frac {d^{2}}{3}}{\binom {3}{3}}B_{0}({\frac {h}{d}})&0&0&0\\{\frac {d^{3}}{4}}{\binom {4}{1}}B_{3}({\frac {h}{d}})&{\frac {d^{3}}{4}}{\binom {4}{2}}B_{2}({\frac {h}{d}})&{\frac {d^{3}}{4}}{\binom {4}{3}}B_{1}({\frac {h}{d}})&{\frac {d^{3}}{4}}{\binom {4}{4}}B_{0}({\frac {h}{d}})&0&0\\{\frac {d^{4}}{5}}{\binom {5}{1}}B_{4}({\frac {h}{d}})&{\frac {d^{4}}{5}}{\binom {5}{2}}B_{3}({\frac {h}{d}})&{\frac {d^{4}}{5}}{\binom {5}{3}}B_{2}({\frac {h}{d}})&{\frac {d^{4}}{5}}{\binom {5}{4}}B_{1}({\frac {h}{d}})&{\frac {d^{4}}{5}}{\binom {5}{5}}B_{0}({\frac {h}{d}})&0\\{\frac {d^{5}}{6}}{\binom {6}{1}}B_{5}({\frac {h}{d}})&{\frac {d^{5}}{6}}{\binom {6}{2}}B_{4}({\frac {h}{d}})&{\frac {d^{5}}{6}}{\binom {6}{3}}B_{3}({\frac {h}{d}})&{\frac {d^{5}}{6}}{\binom {6}{4}}B_{2}({\frac {h}{d}})&{\frac {d^{5}}{6}}{\binom {6}{5}}B_{1}({\frac {h}{d}})&{\frac {d^{5}}{6}}{\binom {6}{6}}B_{0}({\frac {h}{d}})\\\end{pmatrix}}}$
from which the generalized Faulhaber formula is derived:
${\displaystyle S_{h,d}^{r-1}(n)={\frac {d^{r-1}}{r}}\sum _{k=1}^{r}{\binom {r}{k}}B_{r-k}\left({\frac {h}{d}}\right)n^{k}}$
and also the well-known special cases:
{\displaystyle {\begin{aligned}S_{0,1}^{r-1}(n)&={\frac {1}{r}}\sum _{k=1}^{r}{\binom {r}{k}}B_{r-k}(0)n^{k}\\S_{1,1}^{r-1}(n)&={\frac {1}{r}}\sum _{k=1}^{r}{\binom {r}{k}}B_{r-k}(1)n^{k}\end{aligned}}}
where the Bernoulli polynomials calculated in 0 are the Bernoulli numbers and those calculated in 1 are its variant with ${\displaystyle B_{1}}$ changed of sign.[19]

Being ${\textstyle B_{m}\left({\frac {h}{d}}+n\right)=\sum _{k=0}^{m}{\binom {m}{k}}B_{m-k}\left({\frac {h}{d}}\right)n^{k}}$ for the property of translation of Bernoulli's polynomials, the generalized Faulhaber formula can become:

${\displaystyle \sum _{k=0}^{n-1}(h+dk)^{m-1}={\frac {d^{m-1}}{m}}\left(B_{m}\left({\frac {h}{d}}+n\right)-B_{m}\left({\frac {h}{d}}\right)\right)}$
very widespread, unlike the other, in the literature.[14] Hence also the two special cases:
{\displaystyle {\begin{aligned}S_{0,1}^{r-1}&=\sum _{k=0}^{n-1}k^{r-1}={\frac {1}{r}}{\bigg (}B_{r}\left(n\right)-B_{r}\left(0\right){\bigg )}={\frac {B_{r}(n)-B_{r}}{r}}\\S_{1,1}^{r-1}&=\sum _{k=0}^{n-1}(1+k)^{r-1}=\sum _{k=1}^{n}k^{r-1}={\frac {1}{r}}{\bigg (}B_{r}\left(1+n\right)-B_{r}\left(1\right){\bigg )}={\frac {B_{r}(n+1)-B_{r}^{+}}{r}}\end{aligned}}}

## References

1. ^ ${\displaystyle m=0}$ is admissible only when ${\displaystyle 0^{0}}$ is not calculated or if ${\displaystyle 0^{0}=1}$ is set.
2. ^ a b c Beery, Janet (2009). "Sum of powers of positive integers". MMA Mathematical Association of America. doi:10.4169/loci003284.
3. ^ a b Bernoulli, Jacob (1713). "Summae potestatum". Artis conjectandi. Internet Archive. p. 97.
4. ^ Johann Faulhaber (1631). Academia Algebrae - Darinnen die miraculosische Inventiones zu den höchsten Cossen weiters continuirt und profitiert werden. (online copy at Google Books)
5. ^ Jacobi, Carl (1834). "De usu legitimo formulae summatoriae Maclaurinianae". Journal für die reine und angewandte Mathematik. Vol. 12. pp. 263–72.
6. ^ Edwards, Anthony William Fairbank (1982). "Sums of powers of integers: A little of the History". The Mathematical Gazette. 66 (435): 22–28. doi:10.2307/3617302. JSTOR 3617302. S2CID 125682077.
7. ^ The first element of the vector of the sums is ${\displaystyle n}$ and not ${\displaystyle \sum _{k=0}^{n-1}k^{0}}$ because of the first addend, the indeterminate form ${\displaystyle 0^{0}}$, which should otherwise be assigned a value of 1
8. ^ a b c Edwards, A.W.F. (1987). Pascal's Arithmetical Triangle: The Story of a Mathematical Idea. Charles Griffin & C. p. 84. ISBN 0-8018-6946-3.
9. ^ Kalman, Dan (1988). "Sums of Powers by matrix method". Semantic scholar. S2CID 2656552.
10. ^ Helmes, Gottfried (2006). "Accessing Bernoulli-Numbers by Matrix-Operations" (PDF). Uni-Kassel.de.
11. ^ Howard, F.T (1994). "Sums of powers of integers via generating functions" (PDF). CiteSeerX 10.1.1.376.4044.
12. ^ Lang, Wolfdieter (2017). "On Sums of Powers of Arithmetic Progressions, and Generalized Stirling, Eulerian and Bernoulli numbers". arXiv:1707.04451 [math.NT].
13. ^ Tan Si, Do (2017). "Obtaining Easily Sums of Powers on Arithmetic Progressions and Properties of Bernoulli Polynomials by Operator Calculus". Applied Physics Research. 9. Canadian Center of Science and Education. ISSN 1916-9639.
14. ^ a b Bazsó, András; Mező, István (2015). "On the coefficients of power sums of arithmetic progressions". Journal of Number Theory. 153: 117–123. arXiv:1501.01843. doi:10.1016/j.jnt.2015.01.019. S2CID 119140210.
15. ^ a b c d Pietrocola, Giorgio (2021). "Matrici binomiali per insiemi di polinomi calcolanti somme di potenze" [Binomial matrices for sets of polynomials calculating power sums]. Archimede (in Italian). 4: 202–216. ISSN 0390-5543.
16. ^ Guo, Songbai; Shen, Youjien (2013). "On Sums of Powers of Odd Integers". doi:10.3770/j.issn:2095-2651.2013.06.003.
17. ^ Pietrocola, Giorgio (2019). "Didattica delle matrici applicata al classico problema della somma di potenze di interi successivi" [Matrix didactics applied to the classic problem of the sum of powers of successive integers] (PDF) (in Italian). 6° Simposio Mat^Nat APAV Fascino e bellezza della matematica. p. 22.
18. ^ Pietrocola, Giorgio (2022). "Matrici binomiali per insiemi di polinomi calcolanti somme di potenze con basi in p. a." [Binomial matrices for sets of polynomials calculating power sums with bases in arithmetic progression]. MatematicaMente (in Italian). No. 298, 299. ISSN 2037-6367.
19. ^ "A164555, Sequence of numerators in the variant of Bernoulli numbers with ${\displaystyle B_{1}={\frac {1}{2}}}$". OEIS. OEIS Encyclopedia of sequences of integers.