# 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 $S_{h,d}^{m}(n)$i.e. polynomials as a function of $n$ calculating sums of $n$ addends:
    $\sum _{k=0}^{n-1}(h+kd)^{m}=h^{m}+(h+d)^{m}+\cdots +(h+(n-1)d)^{m},$
with $m$ and $n$ integers positive, $h$ first term of an arithmetic progression and $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 $m=1,$ coincides with that of the calculation of the arithmetic series, the sum of the first $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:
    $1+2+\dots +n={\frac {1}{2}}n^{2}+{\frac {1}{2}}n,$ Polynomial $S_{1,1}^{1}(n)$ calculating the sum of the first $n$ natural numbers.
For $m>1,$ the first cases encountered in the history of mathematics are:
    $1+3+\dots +2n-1=n^{2},$ Polynomial $S_{1,2}^{1}(n)$ calculating the sum of the first $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.
    $1^{2}+2^{2}+\ldots +n^{2}={\frac {1}{3}}n^{3}+{\frac {1}{2}}n^{2}+{\frac {1}{6}}n,$ Polynomial $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;
    $1^{3}+2^{3}+\ldots +n^{3}={\frac {1}{4}}n^{4}+{\frac {1}{2}}n^{3}+{\frac {1}{4}}n^{2},$ Polynomial $S_{1,1}^{3}(n)$ calculating the sum of the cubes of the successive integers. Corollary of a theorem of Nicomachus of Gerasa...
L'insieme $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 $m+1$ already knowing the previous ones.
In 1713 the family of Jacob Bernoulli posthumously publishes his Artis Conjectandi 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 for his worthy contributions recognized by Bernoulli himself. It was also immediately clear that the polynomials $S_{0,1}^{m}(n)$ calcolating the sum of $n$ powers of successive integers starting from zero were very similar to those starting from one. This is because it is evident that $S_{1,1}^{m}(n)-S_{0,1}^{m}(n)=n^{m}$and that therefore polynomials of degree $m+1$ of the form ${\frac {1}{m+1}}n^{m+1}+{\frac {1}{2}}n^{m}+\cdots$ subtracted the monomial difference $n^{m}$ they become ${\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 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 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:
    ${\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}}$
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: ${\vec {N}}=A{\vec {S}}$ and multiplying the two sides of the equation to the left by $A^{-1}$ , inverse of the matrix A, we obtain $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 and studying aspects of the problem in their articles useful tools such as the Vandermonde vector. Other researchers continue to explore through the traditional analytic route and generalize the problem of the sum of successive integers to any geometric progression
The coefficients of the polynomials $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 Finally, Edwards' matrix approach was also generalized to any arithmetic progressions

## Solution by matrix method

The general problem has recently been solved 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 $h$ and $d$ which determine the arithmetic progression and a positive integer $m,$ we find $m+1$ polynomials corresponding to the following sums of powers:
$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 $r$ of the triangular matrix $G(h,d)=T(h,d)A^{-1}$ of order $m+1$.
Here is the solving formula in the particular case $m=3$ which gives the polynomials of a given arithmetic progression with exponents from 0 to 3:
    ${\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:
${\vec {S}}_{h,d}(n)=T(h,d)A^{-1}n{\vec {V}}(n)$ or also by placing with $G(h,d)=T(h,d)A^{-1}$
${\vec {S}}_{h,d}(n)=G(h,d)n{\vec {V}}(n)$
Here is the rigorous definition of the matrices and the Vandermonde vector:
${\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 $m=3$ it results therefore
$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 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 $T(h,d)$ on the other hand are the monomials of the power development $(h+d)^{r-1},$ for $r=1,\ldots ,m+1.$.
$T(0,1)$ is the neutral element of the row by column product so that the general equation in this case becomes:
${\vec {S}}_{0,1}(n)=A^{-1}n{\vec {V}}(n)$ that is the one discovered by Edwards
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 $T(h,d)$ after having ascertained the following identity $T(h,d){\vec {V}}(n)=V(h+dn)$

### Sum of powers of successive odd numbers

We use the previous formula to solve the problem of adding powers of successive odds:. The odds correspond to the arithmetic progression with the first element $h=1$ and as reason $d=2.$ We set m = 4 to find the first five polynomials calculating sums of powers of odd. Calculated $T(1,2)$ we obtain:
$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
$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 ${\vec {S}}_{1,2}(n)=G(1,2)n{\vec {V}}(n)$ for $m=4$ and the damage done product:
${\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 ($r=5$) we get then
${S}_{1,2}^{4}={\frac {7}{15}}n-{\frac {8}{3}}n^{3}+{\frac {16}{5}}n^{5}$
and using the other rows:
${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 $m=3$ and calculated $A^{-1}$ and T(1,1) which corresponds to Pascal's triangle:
${\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 $m=3$ and calculated $A^{-1}$ and T(0,1) unit matrix:
${\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 $m=3$, calculated $T(-1,4)$, exploited the result of the previous paragraph and the associative property:
${\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 $G(h,d)$ can be expressed as a function of the Bernoulli polynomials in the following way
$[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 $m=5$ becomes
$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:
$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:
${\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 $B_{1}$ changed of sign.
Being B$_m$( h/d+n)=∑$_{k=0}$$^m$(m/k)B$_{m−k}$(h/d)n$^k${\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:
$\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. Hence also the two special cases:
${\begin{aligned}S_{0,1}^{r-1}&=\sum _{k=0}^{n-1}k^{r-1}={\frac {1}{r}}{\biggB_{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}}{\biggB_{r}\left(1+n\right)-B_{r}\left(1\right){\bigg )}={\frac {B_{r}(n+1)-B_{r}^{+}}{r}}\end{aligned}}$