The problem of calculating the sums of the mth powers of n first integers
is investigated from antiquity by mathematicians around the world.
We learn for examples in the thesis of Coen  as so as in Edwards  that at the beginning of the 11th century ibn al-Haytham had developed the formulae
, , (1.2)
In 15th century his successors had found
About two centuries quietly passed until the day in 1631 when Faulhaber  published at Ulm the prodigious results for sums of odd powers from until in term of powers of . One may find more details on the works of Faulhaber in the reference  .
After Faulhaber, in 1636 French mathematicians Fermat utilizing the figurate numbers and in 1656 Pascal utilizing results of arithmetic triangle, found also recurrence formulae for calculating from lower-order sums  .
Then in 1713 in his posthumous Ars conjectandi, Jacob Bernoulli  , mentioning Faulhaber, published the lists of ten first . It is plausible that from this list he observed that may be written in terms of the numbers which are the same for all m
This famous conjectured formula of Bernoulli was proven in 1755 based on the calculus of finite difference by Euler  , a researcher working with the Bernoulli brothers at Zurich. As informed by Raugh  , Euler had followed de Moivre in given the denomination Bernoulli numbers, introduced the Bernoulli polynomials in 1738 via the generating function
Returning to the Faulhaber conjecture saying that is a polynomial in for all m
we know that Jacobi  has the merit of giving the right proof for this conjecture and moreover calculating the first six Faulhaber coefficients although he did not get a formula for obtaining all of them. Another merit of Jacobi consists in pioneering the use of the derivative with respect to n when observing that
Long years passed until Edwards  showed how to obtain the Faulhaber coefficients by matrix inversion and Knuth  by identification of coefficients of odd powers of n in Bernoulli formula with those in Faulhaber conjecture. Nevertheless the methods of Edwards and Knuth are not easy to apply.
Following Coen, we know the existence of the work Bernoulli numbers: bibliography (1713-1990) of Dilcher  which contained 1956 references by 839 authors!
Concerning the more general problem of powers sums on arithmetic progressions
we remark the recent formula given by Dattoli, Cesarano, Lorenzutta  in term of Bernoulli polynomials
and the formula of Chen, Fu, Zhang  in term of sums of powers of
which are not easy to apply.
After these authors we have proposed a method leading to the formula 
where and .
The Faulhaber formula is thus obtained but we see that the method for obtaining it is cumbersome and the practice calculations of for obtaining the Faulhaber coefficients fastidious.
Rethinking the problem, we observe that an arithmetic progression is a matter of translation, that there is a somehow symmetry between and and so that finally we found a more concise method for resolving the problem and theoretically and practically that we will expose in the following paragraphs.
2. Representation of a Power Sums and Bernoulli Polynomials by Differential Operators
be the powers sums on an arithmetic progression and
the powers sums of the first integers.
We may utilize the shift or translation operator
to get the differential representation
which gives directly the relations
and the generating function
Because (2.4) is valid for all integers n it is also valid for all real and complex values so that we may write
Defining now the set of polynomials by the differential representation
we see that verify
and have the generating function
allowing the identification of them with Bernoulli polynomials defined by Euler  and giving the first of them
, . (2.13)
The Bernoulli polynomials are linked to powers sums according to formulae (2.5), (2.8), (2.9) by the relations
which lead to the followed beautiful formula where the second member does not depend on n
Besides it leads also to the formula
which jointed with (2.15) gives rise to the historic Jacobi conjectured formula 
From (2.18) we see that are identifiable with Bernoulli numbers .
3. The Powers Sums
1) Powers sums in terms of Bernoulli polynomials and powers of n
From the Equation (2.16) and the boundary condition we get immediately the solution of (2.16)
which may be put under the algorithmic form very easy to remember
or taking into account (2.11)
Putting in (3.2) and replacing with we recognize the famous Bernoulli formula 
2) The Faulhaber formula on powers sums
In instead of utilizing z and n for arguments let us utilize
Equation (2.16) becomes
Happily from the definition of by differential operators (2.9) we may write down
which jointed with (2.10) leads to the important properties
As an polynomial of order 2k having for roots may be put by identification of coefficients under the form of a homogeneous polynomial of order k in we obtain the important property:
is a homogeneous polynomial of order k in (3.14)
By the way we notify that because for it is also a homogeneous polynomial of order in as conjectured Faulhaber and proven somehow by Jacobi.
Moreover the calculations of the quoted polynomials in Z or in u may be done thank to the hereinafter Table 3.
Jointed (3.14) with the formula came from the Jacobi formula (2.18)
we may define a polynomial of order k depending on Z such that
The definition of by (3.16) differs a little in index with its definition in  .
Thank to these considerations we get the solution of (3.9) corresponding to the boundary condition under the form
or under the algorithmic form
where we see that the Faulhaber coefficients are successive derivations beginning from the first one of the power sums on integers writing under the form .
Curiously by replacing z with n and consequently Z with in the definition we get the very important and very interesting formula linking the Faulhaber powers sums on integers and on arithmetic progressions
4. Obtaining Practically Bernoulli Numbers, Bernoulli Functions and Powers Sums
4.1. Calculations of Bernoulli Numbers
For calculating we remark that from
and (2.10) we get the recursion relation
which may be written under the matrix form
This matrix equation may be resolved by doing by hand or by program linear combinations over lines from the second one in order to replace them with lines each containing only some non-zero rational numbers.
For instance for calculating successively we may utilize the matrix equation
Matrix equation for calculating .
We remark that the last line of this matrix has replaced the line Some results are
, , ,
4.2. Calculations of Bernoulli Polynomials
For calculating Bernoulli polynomials, we remark that (2.11) and the Jacobi formula (2.18) give rise to the relations
which, knowing, permits the easy calculations of all the Bernoulli polynomials and the powers sums on integers from the set of Bernoulli numbers as we may see hereafter in Table 1.
4.3. Calculations of Powers Sums as Polynomials in n
For calculating we utilize the algorithmic formula (3.2) and get the results shown in Table 2.
4.4. Calculations of Faulhaber Powers Sums
Practically for transforming into from which one obtains the Faulhaber coefficients let us remark that
so that we may replace the first term of with minus a polynomial in z. The polynomial in z so obtained must begin with a term in and not so that we may continue to replace in it with a term in
Table 1. Obtaining , from .
Table 2. Obtaining from .
minus a polynomial in z. The operations continue so on until finish.
By this operation we observe that we may omit all odd powers terms in before and during the transformation of into for . From these remarks we may establish Table 3 of transformation.
This algorithm for obtaining gives rise astonishingly to the easy calculation of the Faulhaber formula for .
As examples we have
Table 3. Change of in change of into.
5. Formula for Even Powers Sums
By derivation of with respect to z and remarking that
we obtain the formulae giving in function of
・ , ,
6. Remarks and Conclusions
The main particularity of this work consists in obtaining as the transform of by a differential operator, as so as the Bernoulli polynomial , from which we deduce the new formula and get immediately as polynomials in n. From which we get also the property saying that is a homogeneous polynomial of order k in .
The second particularity consists in performing the change of arguments from z into and n into in order to get the relation which gives rise to the Faulhaber formula of . From this formula we see that the Faulhaber coefficients are successive derivatives of the function where are powers sums on integers.
The relation leads also to the Faulhaber formula where .
The third particularity of this work is proposing a method for obtaining easily the Bernoulli numbers from a matrix equation, then of
As conclusion we think that the calculations of powers sums are greatly facilitated by the utilization of derivation operators and the translation operator , parts of Operator Calculus.
Operator calculus, which is very different from Heaviside operational calculus, is thus merited to be known and introduced into Mathematical Analysis. Moreover it just had a solid foundation and many interesting applications for instance in the domains of Special functions, Differential equations, Fourier, Laplace and other transforms, quantum mechanics  .
The author acknowledges the reviewer for comments leading to the clear proof of the important formula (3.14). The author also acknowledges Mrs. Beverly GUO for her perseverance in the refinement of the article representation.
He reiterates his gratitude toward his adorable wife for all the cares she devoted for him during the long months he performed this work.
 Edwards, A.W.F. (1982) Sums of Powers of Integers: A Little of the History, The Mathematical Gazette, 66, 22-28. https://www.jstor.org/stable/3617302
 Raugh, M. (2014) The Bernoulli Summation Formula: A Pretty Natural Derivation. A Presentation for Los Angeles City College at The High School Math Contest. http://www.mikeraugh.org/