Selection of Coherent and Concise Formulae on Bernoulli Polynomials-Numbers-Series and Power Sums-Faulhaber Problems

Do Tan Si

The HoChiMinh-City Physical Association, Ho Chi Minh City, Vietnam.

DOI: 10.4236/am.2022.1310051                

Abstract

Utilizing the translation operator to represent Bernoulli polynomials and power sums as polynomials of Sheffer-type, we obtain concisely almost all their known properties as so as many new ones, especially new recursion relations for calculating Bernoulli polynomials and numbers, new formulae for obtaining power sums of entire and complex numbers. Then by the change of arguments from z into Z = z(z-1) and n into λ which is the 1^st order power sum we obtain the Faulhaber formula for powers sums in term of polynomials in λ having coefficients depending on Z. Practically we give tables for calculating in easiest possible manners, the Bernoulli numbers, polynomials, the general powers sums.

Keywords

Bernoulli Numbers, Bernoulli Polynomials, Powers Sums, Zeta Function, Faulhaber Conjecture

1. Introduction

In many branches of mathematics the problem of Bernoulli numbers related to the millenary problem of power sums is probably the most studied since the publication of the book Ars Conjectandi by Euler in 1738 [1] as we can see on the net and, specially, in a didactical thesis of Coen [2] , the explicative work of Raugh [3] , Beardon [4] , the bibliography of thousands of articles on Bernoulli numbers realized by Dilcher, Shula, Slavutskii [5] , etc.

Concerning Bernoulli polynomials $B_m\left(z\right)$ , classically defined from a generating function, there had not so much properties, the most remarkable is its representation by a hyper-differential operator, the Hurwitz expansion of them into Fourier series, the Roman formula for $B_m\left(nz\right)$ , the Euler-McLaurin formula, etc.

As for the power sums on real and complex numbers, including the famous Faulhaber conjecture, there has no valuable formula linking them with Bernoulli polynomials until only some years ago [6] .

Regarding the situation, we would like to perform a selection of as many as possible known and new interesting properties of Bernoulli polynomials then of Bernoulli numbers in a coherent way, i.e., by only one approach, which utilizes principally operator calculus lying on the couple of operators position and derivation, similar as the couple $\stackrel{\to }{r},\nabla$ in quantum mechanics.

In Section 2, we will treat the problem of Bernoulli polynomials, from their representation by a hyper-differential operator to almost all of their algebraic properties to the fact that $B_m\left(n\right)$ is equal to the primitive of the power sums of natural integers. Afterward we show that the formula giving Bernoulli polynomials of a sum of two arguments $B_m\left(z+y\right)$ leads to two new recurrence relations for obtaining $B_m\left(z\right)$ . We also give another approach for calculating without integrations, the Fourier series of Bernoulli polynomials and the Bernoulli series of functions, the relation of $B_m\left(z\right)$ with the Euler zeta function. Afterward we show an up-to-date procedure for obtaining $B_m\left(z\right)$ from and only from $B_{m-1}\left(z\right)$ leading to the rapid establishment of Table of Bernoulli polynomials and numbers. Finally, we show a new way for obtaining Fourier series of Bernoulli polynomials, Euler zeta function, and vice-versa, the series of functions in term of a set of Bernoulli polynomials.

In Section 3, we treat the problem of Bernoulli numbers $B_m$ , from its initial definition by Jakob Bernoulli in 1713 who related them by conjecture with the power sums on natural numbers. By comparison of this relation with the preceding formula linking $B_m\left(n\right)$ with power sums, we may identify $B_m$ with $B_m\left(0\right)$ then calculate $B_m$ by a simple matrix method side-by-side with the method, more powerful, link with l, coming from the special recurrence formula coming from $B_m\left(z+y\right)$ .

In Section 4, we prove by utilizing the translation operator ${\text{e}}^{a{\partial }_z}$ , coming from the Newtonian binomial, that the power sums on complex numbers are simply related to those on natural numbers. On the other hand, we prove that they are also related very simply to Bernoulli polynomials, from that we get again the recurrence relation between Bernoulli polynomials.

Section 5 is devoted to the Faulhaber problem regarding power sums on complex numbers. Here we show that power sums on complex numbers may be calculated from sums of entire numbers somehow by writing $B_{2m}\left(z\right)$ in function of the new argument $Z=z\left(z-1\right)$ .

2. Bernoulli Polynomials

2.1. Definition and Principal Properties

In 1738, Euler introduced the Bernoulli polynomials $B_m\left(z\right)$ via the generating function [1]

$\frac{t}{{\text{e}}^t-1}{\text{e}}^{zt}={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}\frac{1}{m!}{B}_m\left(z\right)t^m}$ (2.1)

which directly gives by identification

$B_0\left(z\right)=1$ , $B_1\left(z\right)=z-\frac{1}{2}$ , $B_2\left(0\right)=\frac{1}{6}$ (2.2)

Utilizing the translation operator ${\text{e}}^{a{\partial }_z}$ coming from the Newtonian binomial

${\left(x+a\right)}^m={\displaystyle \underset{k=0}{\overset{m}{\sum }}\left(\begin{array}{c}m\\ k\end{array}\right)a^k{x}^{m-k}}={\displaystyle \underset{k=0}{\overset{m}{\sum }}\frac{a^k}{k!}{\partial }_z^k{x}^m}={\text{e}}^{a{\partial }_z}{x}^m$ (2.3)

and having the property

${\text{e}}^{a{\partial }_z}f\left(x\right)=f\left(x+a\right)$

${\text{e}}^{{\partial }_z}{\text{e}}^{tz}={\text{e}}^{t\left(z+1\right)}={\text{e}}^t{\text{e}}^{tz}$ (2.4)

$\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{\text{e}}^{tz}=\frac{t}{{\text{e}}^t-1}{\text{e}}^{tz}$ (2.5)

we directly find from (2.1) that $B_m\left(z\right)$ is the transform of $z^m$ via a differential operator

$B_m\left(z\right)=\frac{{\partial }_z}{e^{{\partial }_z}-1}{z}^m,m>0$ (2.6)

From (2.6) we get the famous known formulae

${B^{\prime }}_m\left(z\right)=m{B}_{m-1}\left(z\right)$ (2.7)

$B_m\left(z+1\right)-B_m\left(z\right)=\left({\text{e}}^{{\partial }_z}-1\right)B_m\left(z\right)={\partial }_z{z}^m=m{z}^{m-1}$ (2.8)

$B_m\left(1\right)-B_m\left(0\right)={\delta }_{m1}$ (2.9)

and the following formula which gives ${\text{e}}^{inxz}$ as series of Bernoulli polynomials.

${\displaystyle \underset{m=0}{\overset{\infty }{\sum }}\frac{t^m{B}_m\left(z\right)}{m!}}=\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{\text{e}}^{t\text{}z}=\frac{t}{{\text{e}}^t-1}{\text{e}}^{t\text{}z}$ (2.10)

From (2.8) we get the formula given by Roman [7]

$B_{m+1}\left(\frac{z}{y}+N\right)-B_{m+1}\left(\frac{z}{y}\right)=\left(m+1\right){\displaystyle \underset{n=0}{\overset{N-1}{\sum }}{\left(\frac{z}{y}+n\right)}^m}$ (2.11)

From (2.10) we get the formulae on relations of Bernoulli polynomials versus trigonometric functions, especially the Castellanos formula [8]

${\displaystyle \underset{m=2}{\overset{\infty }{\sum }}\frac{{\left(2ix\right)}^m{B}_m\left(0\right)}{m!}}=\frac{x\cos x}{\sin x}-1$ (2.12)

The formulae (2.7) and (2.8) give the important formulae

${\displaystyle {\int }_0^1{B}_{m-1}\left(z\right)\text{d}z}=\frac{1}{m}\left(B_m\left(1\right)-B_m\left(0\right)\right)={\delta }_{m1}$ (2.13)

$B_m\left(n\right)-B_m\left(0\right)=m\left(0^{m-1}+1^{m-1}+\cdots +{\left(n-1\right)}^{m-1}\right)$ (2.14)

and the Taylor expansion

$B_m\left(z\right)=B_m\left(a\right)+\cdots +\left(\begin{array}{c}m\\ k\end{array}\right){\left(z-a\right)}^k{B}_{m-k}\left(a\right)+\cdots +{\left(z-a\right)}^m{B}_0(\; a\; )$

which may be put under symbolic form

$B_m\left(z+a\right)=:{\left(B\left(a\right)+z\right)}^m$ (2.15)

where undefined symbols $B^k\left(a\right)$ are to be replaced with $B_k\left(a\right)$ .

Exploring now the inter-relations between Bernoulli polynomials.

From (2.4) and (2.7) we get the complementary of (2.15)

$\begin{array}{c}{B}_m\left(z+a\right)={\text{e}}^{a{\partial }_z}{B}_m\left(z\right)=\left(1+\cdots +\frac{a^k}{k!}{\partial }^k+\cdots +\frac{a^m}{m!}{\partial }^m\right)B_m\left(z\right)\\ =B_m\left(z\right)+\cdots +\left(\begin{array}{c}m\\ k\end{array}\right)B_{m-k}\left(z\right)a^k+\cdots +B_0\left(z\right)a^m\\ =:{\left(B\left(z\right)+a\right)}^m,0^0=1\end{array}$ (2.16)

From (2.13)

${\displaystyle {\int }_z^{z+1}{B}_m\left(y\right)\text{d}y}={\displaystyle {\int }_0^1{B}_m\left(z+y\right)\text{d}y}=z^m{\displaystyle {\int }_0^1{B}_0\left(y\right)\text{d}y}=z^m$

${\displaystyle {\int }_0^1{B}_m\left(y\right)\text{d}y}=0^m$

${\displaystyle {\int }_0^n{B}_m\left(y\right)\text{d}y}=0^m+1^m+\cdots +{\left(n-1\right)}^m$ (2.17)

i.e.,

“The sum of powers of order m of n first entire numbers from 0 to $\left(n-1\right)$ , denoted by $S_m\left(n\right)$ , is equal to the simple primitive (without constant of integration) of the Bernoulli polynomial $B_m\left(n\right)$ ” and vice-versa,

“The Bernoulli polynomial $B_m\left(n\right)$ is equal to the derivative of the power sums $S_m\left(n\right)$ ”

As for $B_m\left(-z\right)$ we see that

$\begin{array}{c}{B}_m\left(-z\right)=\frac{-{\partial }_z}{1-{\text{e}}^{-{\partial }_z}}{\left(-z\right)}^m={\left(-1\right)}^m{\text{e}}^{{\partial }_z}\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{z}^m\\ ={\left(-1\right)}^m{B}_m\left(z+1\right)=:{\left(-1\right)}^m{\left(B\left(z\right)+1\right)}^m\end{array}$ (2.18)

which leads to

$B_m\left(-z+\frac{1}{2}\right)={\left(-1\right)}^m{B}_m\left(z+\frac{1}{2}\right)$ (2.19)

i.e., to the theorem

“The graph of a Bernoulli polynomial is symmetric with respect to the axis $z=\frac{1}{2}$ if m is pair and anti-symmetric if m is impair”.

Joint (2.19) with (2.9) we get the famous property [1]

$B_{2m+1}\left(1\right)=-B_{2m+1}\left(0\right)=\frac{1}{2}{\delta }_{m0}$ (2.20)

Now, by replacing in (2.6) z with $\frac{z}{n}$ so that ${\partial }_z$ is with $n{\partial }_z$ we get

$B_m\left(\frac{z}{n}\right)=\frac{n{\partial }_z}{{\text{e}}^{n{\partial }_z}-1}{\left(\frac{z}{n}\right)}^m=\frac{n{\partial }_z}{\left({\text{e}}^{{\partial }_z}-1\right)\left(1+{\text{e}}^{{\partial }_z}+{\text{e}}^{2{\partial }_z}+\cdots +{\text{e}}^{\left(n-1\right){\partial }_z}\right)}{\left(\frac{z}{n}\right)}^m$

and the formula

${\displaystyle \underset{k=0}{\overset{n-1}{\sum }}{B}_m\left(\frac{z+k}{n}\right)}=n^{1-m}{B}_m\left(z\right)$ (2.21)

saying that

$B_m\left(z\right)$ is $n^{m-1}$ times the sum of $B_m\left(\frac{z+k}{n}\right),k<n$

For examples:

$2^{1-m}{B}_m\left(z\right)=B_m\left(\frac{z}{2}\right)+B_m\left(\frac{z+1}{2}\right)$

$B_{2m+1}\left(\frac{1}{3}\right)+B_{2m+1}\left(\frac{2}{3}\right)=0=\left(3^{-2m}-1\right)B_{2m+1}(\; 1\; )$

By replacing in (2.6) z with nz and ${\partial }_z$ with $\frac{1}{n}{\partial }_z$ we find again the formula given by Raabe [9] in 1851

$B_m\left(nz\right)=n^{m-1}\left(B_m\left(z\right)+B_m\left(z+\frac{1}{n}\right)+\cdots +B_m\left(z+\frac{n-1}{n}\right)\right)$ (2.22)

saying that

“ $B_m\left(nz\right)$ is $n^{m-1}$ times the sum of $B_m\left(z+\frac{k}{n}\right),k<n$ .”

For examples

$B_m\left(2z\right)=2^{m-1}\left(B_m\left(z\right)+B_m\left(z+\frac{1}{2}\right)\right)$

$B_m\left(\frac{1}{2}\right)=\left(2^{1-m}-1\right)B_m(\; 0\; )$

$B_1\left(3z\right)=B_1\left(z\right)+B_1\left(z+\frac{1}{3}\right)+B_1\left(z+\frac{2}{3}\right)$

$5^{-m}{B}_m\left(0\right)=\left(1+{\left(-1\right)}^m\right)\left(B_m\left(\frac{1}{5}\right)+B_m\left(\frac{2}{5}\right)\right)$

2.2. Bernoulli Polynomials of Sum of Two Arguments

From the following property of operators that we characterize fundamental [10]

$f\left({\partial }_z\right)g\left(z\right)\equiv g\left(z\right)f\left({\partial }_z\right)+\frac{1}{1!}{g}^{\prime }\left(z\right)f^{\prime }\left({\partial }_z\right)+\frac{1}{2!}{g}^{″}\left(z\right)f^{″}\left({\partial }_z\right)+\cdots$ (2.23)

we get

$\begin{array}{c}{B}_{m+1}\left(z\right)=\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}z{z}^m=z\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{z}^m+\left(\frac{1}{{\text{e}}^{{\partial }_z}-1}-\frac{{\partial }_z{\text{e}}^{{\partial }_z}}{{\left({\text{e}}^{{\partial }_z}-1\right)}^2}\right)z^m\\ =z{B}_m\left(z\right)+\frac{1}{{\text{e}}^{{\partial }_z}-1}{z}^m-\frac{{\text{e}}^{{\partial }_z}{\partial }_z}{{\left({\text{e}}^{{\partial }_z}-1\right)}^2}{z}^m\end{array}$

${\partial }_z{B}_{m+1}\left(z\right)={\partial }_{z}z{B}_m\left(z\right)+B_m\left(z\right)-\frac{{\text{e}}^{{\partial }_z}{\partial }_z^2}{{\left({\text{e}}^{{\partial }_z}-1\right)}^2}{z}^m$

$\left(m-1\right)B_m\left(z\right)=z{\partial }_z{B}_m\left(z\right)-{\text{e}}^{{\partial }_z}\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{B}_m(\; z\; )$

Now, because

${\partial }_{z+y}f\left(z+y\right)={\partial }_{z}f\left(z+y\right)={\partial }_{y}f\left(z+y\right)$ (2.24)

$\begin{array}{c}\left(m-1\right)B_m\left(z+y\right)=m\left(z+y\right)B_{m-1}\left(z+y\right)-{\text{e}}^{{\partial }_y}\frac{{\partial }_y}{{\text{e}}^{{\partial }_y}-1}{B}_m\left(z+y\right)\\ =:m\left(z+y\right)B_{m-1}\left(z+y\right)-{\left(B\left(z\right)+B\left(y+1\right)\right)}^m\end{array}$ (2.25)

The above recurrence formula is to be compare with that given by Weisstein [11] without proof where there seems has a little mistake

$\left(1-m\right)B_m\left(z+y\right)+m\left(z+y-1\right)B_{m-1}\left(z+y\right)=:{\left(B\left(z\right)+B\left(y\right)\right)}^m$

From (2.25) and knowing that $B_k\left(1\right)={\left(-1\right)}^k{B}_k\left(0\right)$ we obtain another type of recurrence formula for Bernoulli polynomials

$\left(m-1\right)B_m\left(z\right)=mz{B}_{m-1}\left(z\right)-{\left(B\left(z\right)+B\left(1\right)\right)}^m$ (2.26)

$B_m\left(z\right)=B_1\left(z\right)B_{m-1}\left(z\right)-\frac{1}{m}{\displaystyle \underset{k=2}{\overset{m}{\sum }}{\left(-1\right)}^k\left(\begin{array}{c}m\\ k\end{array}\right)B_k\left(0\right)B_{m-k}\left(z\right)}$ (2.27)

For examples, with $B_1\left(z\right)=z-\frac{1}{2}$ , $B_2\left(0\right)=\frac{1}{6}$ ,

$B_2\left(z\right)=B_1\left(z\right)B_1\left(z\right)-\frac{1}{2}{B}_2\left(0\right)B_0\left(z\right)={\left(z-\frac{1}{2}\right)}^2-\frac{1}{12}=z^2-z+\frac{1}{6}$

$\begin{array}{c}{B}_3\left(z\right)=B_1\left(z\right)B_2\left(z\right)-\frac{1}{3}3B_2\left(0\right)B_1\left(z\right)\\ =\left(z-\frac{1}{2}\right)\left(z^2-z+\frac{1}{6}-\frac{1}{6}\right)=z^3-\frac{3}{2}{z}^2+\frac{1}{2}z\end{array}$

$B_4\left(z\right)=B_1\left(z\right)B_3\left(z\right)-\frac{1}{4}\left(B_2\left(z\right)+B_4(\; 0\; )\right)$

2.3. The Fourier Series of Bernoulli Polynomials. Euler Zeta Function. Powers of pi

By successive integrations by parts and utilizing the formula (2.13) for $n,m\ge 1$ we get, knowing (2.9),

$\begin{array}{l}{\displaystyle {\int }_0^1{B}_n\left(z\right)B_m\left(z\right)\text{d}z}=\frac{1}{m+1}{\displaystyle {\int }_0^1{B}_n\left(z\right){B^{\prime }}_{m+1}\left(z\right)\text{d}z}\\ =\frac{1}{m+1}{\left(B_n\left(z\right)B_{m+1}\left(z\right)\right)|}_0^1-\frac{n}{m+1}{\displaystyle {\int }_0^1{B}_{n-1}\left(z\right)B_{m+1}\left(z\right)\text{d}z}\\ ={\left(-1\right)}^{n-1}\frac{n!m!}{\left(m+n\right)!}{\left(B_1\left(z\right)B_{m+n}\left(z\right)\right)|}_0^1\\ ={\left(-1\right)}^{n-1}\frac{n!m!}{\left(m+n\right)!}{B}_{m+n}\left(0\right)\end{array}$ (2.28)

Because of the factor ${\left(-1\right)}^{n-1}$ we may conclude that

$B_{2n+1}\left(0\right)=0$ for $n>0$ and $B_{2n+2}\left(0\right)$ has opposite sign with respect to $B_{2n}\left(0\right)$ .

The same method also gives

$\begin{array}{c}{\displaystyle {\int }_0^1{B}_m\left(z\right){\text{e}}^{-2ik\pi \text{}z}\text{d}z}=\frac{-1}{2\pi ik}{\displaystyle {\int }_0^1{B}_m\left(z\right)}{\left({\text{e}}^{-2ik\pi \text{}z}\right)}^{\prime }\text{d}z\\ =\frac{-1}{2\pi ik}{\delta }_{m1}-\frac{-m}{2\pi ik}{\displaystyle {\int }_0^1{B}_{m-1}\left(z\right){\text{e}}^{-2i\pi kz}\text{d}z}\\ =\cdots =-\frac{m!}{{\left(2\pi ik\right)}^m}\end{array}$ (2.29)

which provides us the following formula on Fourier series of $B_m\left(z\right)$ proven by Hurwitz in 1890 by another method [10]

$B_m\left(z\right)={\displaystyle \underset{k\in Z,k\ne 0}{\overset{\infty }{\sum }}\left({\displaystyle {\int }_0^1{B}_m\left(z\right){\text{e}}^{-2ik\pi \text{}z}\text{d}z}\right){\text{e}}^{i2\pi kz}}=-\frac{m!}{{\left(2i\pi \right)}^m}{\displaystyle \underset{k\in Z,k\ne 0}{\overset{\infty }{\sum }}\frac{1}{k^m}{\text{e}}^{i2\pi kz}},0\le z\le 1$ (2.30)

2.4. Bernoulli Series of Functions

Let $f\left(z\right)$ be a periodic function defined on an interval $a\le z<b$ and has the period $P=b-a$ . For expanding $f\left(z\right)$ into a Fourier series of exponentials

$f\left(z\right)={\displaystyle \underset{n\in Z}{\sum }c\left(n\right){\text{e}}^{i2\pi n\frac{z}{P}}}$ , $a\le z<b=a+P$ (2.31)

we firstly write

${\displaystyle {\int }_a^{b^{-}}{\text{e}}^{-i2\pi n_0\frac{z}{P}}}f\left(z\right)\frac{\text{d}z}{P}={\displaystyle \underset{n\in Z}{\sum }c\left(n\right){\displaystyle {\int }_a^{b^{-}}{\text{e}}^{i2\pi \left(n-n_0\right)\frac{z}{P}}}\frac{\text{d}z}{P}}$

and see that the second member is equal uniquely to $c\left(n_0\right)$ so that

$c\left(n_0\right)=\frac{1}{P}{\displaystyle {\int }_a^{b^{-}}{\text{e}}^{-i2\pi n_0\frac{z}{P}}f\left(z\right)\text{d}z}$ (2.32)

The Fourier series of a function, if it exists, is then

$f\left(z\right)=\frac{1}{P}{\displaystyle \underset{n\in Z}{\overset{\infty }{\sum }}{\text{e}}^{i2\pi n\frac{z}{P}}}{\displaystyle {\int }_a^{b^{-}}{\text{e}}^{-i2\pi n\frac{z}{P}}}f\left(z\right)\text{d}z$ (2.33)

To avoid integrations in the calculation, we may utilize the method of integrations by parts and get

$c\left(n\right)=\frac{1}{P}{\displaystyle {\int }_a^{b^{-}}f\left(z\right){\text{e}}^{-i2\pi n\frac{z}{P}}\text{d}z}$

$Pc\left(n\right)=\frac{-P}{2i\pi n}\left(f\left(b^{-}\right){\text{e}}^{-i2\pi n\frac{b^{-}}{P}}-f\left(a\right){\text{e}}^{-i2\pi n\frac{a}{P}}\right)+\frac{P}{2i\pi n}{\displaystyle {\int }_a^{b^{-}}{f}^{\prime }\left(z\right){\text{e}}^{-i2\pi n\frac{z}{P}}}\text{d}z$

$Pc\left(n\right)=-{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\left(\frac{P}{2i\pi n}\right)}^{k+1}}\left(f^{\left(k\right)}\left(b^{-}\right){\text{e}}^{-i2\pi n\frac{b^{-}}{P}}-f^{\left(k\right)}\left(a\right){\text{e}}^{-i2\pi n\frac{a}{P}}\right)-0(\; k\; )$

so that we may write down the Fourier series formula

$f\left(z\right)=\frac{1}{P}{\displaystyle {\int }_a^{b^{-}}f\left(z\right)\text{d}z}-\frac{1}{P}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{f^{\left(k\right)}\left(z\right)|}_a^{b^{-}}}{\displaystyle \underset{n\in Z,n\ne 0}{\overset{}{\sum }}{\left(\frac{P}{2i\pi n}\right)}^{k+1}{\text{e}}^{i2\pi n\frac{z-a}{P}}}$ (2.34)

In the case $0\le z<1$ , jointed the preceding formula written under the form

$f\left(z\right)={\displaystyle {\int }_0^{1-}f\left(z\right)\text{d}z}-{\displaystyle \underset{n\in Z,n\ne 0}{\overset{}{\sum }}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{f^{\left(k\right)}\left(z\right)|}_0^{1-}}}{\left(\frac{1}{2i\pi n}\right)}^{k+1}{\text{e}}^{i2\pi nz}$

with the Hurwitz formula we get the new and precious formula on expansion of derivable functions into series of Bernoulli polynomials

$f\left(z\right)={\displaystyle {\int }_0^{1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\left[f^{\left(k\right)}\left(1\right)-f^{\left(k\right)}\left(0\right)\right]}\frac{1}{\left(k+1\right)!}{B}_{k+1}\left(z\right)$ (2.35)

or

$f\left(z\right)={\displaystyle {\int }_0^{1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{N}{\sum }}{f^{\left(k\right)}\left(z\right)|}_0^1}\frac{B_{k+1}\left(z\right)}{\left(k+1\right)!}-{\displaystyle \underset{\begin{array}{l}n\in Z\\ n\ne 0\end{array}}{\overset{}{\sum }}{\displaystyle \underset{k=N+1}{\overset{\infty }{\sum }}{f^{\left(k\right)}\left(z\right)|}_0^1}{\left(\frac{1}{2i\pi n}\right)}^{k+1}{\text{e}}^{i2\pi nz}}$ (2.36)

For examples, under matrix form

$\left(\begin{array}{c}f\left(z\right)\\ 1\\ z\\ z^2\\ ⋮\\ z^m\end{array}\right)=\left(\begin{array}{cccccc}{\left[{\displaystyle \int f\left(z\right)}\right]}_0^1& {\left[f\left(z\right)\right]}_0^1& {\left[f^{\prime }\left(z\right)\right]}_0^1& {\left[f^{″}\left(z\right)\right]}_0^1& \cdots & {\left[f^{\left(m-1\right)}\left(z\right)\right]}_0^1\\ 1& & & & & \\ 1/2& 1& & & & \\ 1/3& 1& 2& & & \\ ⋮& ⋮& ⋮& ⋮& \ddots & \\ 1/\left(m+1\right)& 1& m& m\left(m-1\right)& \cdots & m!\end{array}\right)\left(\begin{array}{c}1\\ B_1\left(z\right)/1!\\ B_2\left(z\right)/2!\\ B_3\left(z\right)/3!\\ ⋮\\ B_m\left(z\right)/m!\end{array}\right)$ (2.37)

to be compared with

$\left(\begin{array}{c}f\left(z\right)\\ 1\\ z\\ z^2\\ ⋮\\ z^m\end{array}\right)=\left(\begin{array}{cccccc}{\left[{\displaystyle \int f\left(z\right)}\right]}_0^1& {\left[f\left(z\right)\right]}_0^1& {\left[f^{\prime }\left(z\right)\right]}_0^1& {\left[f^{″}\left(z\right)\right]}_0^1& \cdots & {\left[f^{\left(m-1\right)}\left(z\right)\right]}_0^1\\ 1& & & & & \\ 1/2& 1& & & & \\ 1/3& 1& 2& & & \\ ⋮& ⋮& ⋮& ⋮& \ddots & \\ 1/\left(m+1\right)& 1& m& m\left(m-1\right)& \cdots & m!\end{array}\right)\left(\begin{array}{c}1\\ {\displaystyle \underset{n\ne 0}{\sum }{\text{e}}^{2i\pi nz}/\left(2i\pi n\right)}\\ {\displaystyle \underset{n\ne 0}{\sum }{\text{e}}^{2i\pi nz}/{\left(2i\pi n\right)}^2}\\ {\displaystyle \underset{n\ne 0}{\sum }{\text{e}}^{2i\pi nz}/{\left(2i\pi n\right)}^3}\\ ⋮\\ {\displaystyle \underset{n\ne 0}{\sum }{\text{e}}^{2i\pi nz}/{\left(2i\pi n\right)}^m}\end{array}\right)$ (2.38)

Formula (2.36) leads also to

$f\left(0\right)={\displaystyle {\int }_0^{1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\left[f^{\left(k\right)}\left(1\right)-f^{\left(k\right)}\left(0\right)\right]}\frac{1}{\left(k+1\right)!}{B}_{k+1}\left(0\right)$ (2.39)

$f^{\prime }\left(z\right)={\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\left[f^{\left(k\right)}\left(1\right)-f^{\left(k\right)}\left(0\right)\right]}\frac{1}{\left(k+1\right)!}{B}_{k+1}\left(z\right)$ (2.40)

As first interesting applications

${\text{e}}^z=\left(e-1\right)+\left(e-1\right){\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{B_{k+1}\left(z\right)}{\left(k+1\right)!}}0\le z<1$

$\frac{e-2}{e-1}=-{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{B_{k+1}\left(0\right)}{\left(k+1\right)!}}$ (2.41)

By (2.36) we also obtain a precious recurrence formula of Bernoulli polynomials

$z^m={\displaystyle {\int }_0^1{z}^m\text{d}z}+{\displaystyle \underset{k=1}{\overset{m}{\sum }}\left(\begin{array}{c}m\\ k-1\end{array}\right)}\frac{B_k\left(z\right)}{k}0\le z<1$ (2.42)

i.e., under matrix form

$\left(\begin{array}{c}z\\ z^2\\ z^3\\ z^4\\ ⋮\end{array}\right)=\left(\begin{array}{c}{2}^{-1}\\ 3^{-1}\\ 4^{-1}\\ 5^{-1}\\ ⋮\end{array}\right)+\left(\begin{array}{ccccc}1& & & & \cdots \\ 1& 2& & & \cdots \\ 1& 3& 3& & \cdots \\ 1& 4& 6& 4& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right)\left(\begin{array}{c}{B}_1\left(z\right)/1\\ B_2\left(z\right)/2\\ B_3\left(z\right)/3\\ B_4\left(z\right)/4\\ ⋮\end{array}\right)$ (2.43)

which may be resolved for $B_m\left(z\right)$ and $B_m\left(0\right)$ my matrix calculus.

2.5. Obtaining $B_m\left(z\right)$ from $B_{m-1}\left(z\right)$ and Table of Bernoulli Polynomials

Integrating two times as followed the Hurwitz formula on Fourier series of Bernoulli polynomials we get

${\displaystyle {\int }_0^x{B}_m\left(z\right)\text{d}z}=-m!{\displaystyle \underset{\begin{array}{l}n\in Z\\ n\ne 0\end{array}}{\overset{\infty }{\sum }}{\left(\frac{1}{2i\pi n}\right)}^{m+1}\left({\text{e}}^{i2\pi nz}-1\right)}=\frac{1}{m+1}{B}_{m+1}\left(z\right)+m!{\displaystyle \underset{\begin{array}{l}n\in Z\\ n\ne 0\end{array}}{\overset{\infty }{\sum }}{\left(\frac{1}{2i\pi n}\right)}^{m+1}}$

${\displaystyle {\int }_0^1\text{d}z{\displaystyle {\int }_0^z{B}_m\left(x\right)\text{d}x}}=m!{\displaystyle \underset{\begin{array}{l}n\in Z\\ n\ne 0\end{array}}{\overset{\infty }{\sum }}{\left(\frac{1}{2i\pi n}\right)}^{m+1}}$

$B_{m+1}\left(z\right)=\left(m+1\right){\displaystyle {\int }_0^z{B}_m\left(x\right)\text{d}x}-\left(m+1\right){\displaystyle {\int }_0^1\text{d}z{\displaystyle {\int }_0^z{B}_m\left(x\right)\text{d}x}}$ (2.44)

i.e.,

$B_{m+1}\left(z\right)$ is equal to $\left(m+1\right)$ times the primitive of $B_m\left(z\right)$ minus the double primitive of $B_m\left(z\right)$ calculated for $z=1$ . The second term is so equal to $B_m\left(0\right)={\left(-1\right)}^m{B}_m\left(1\right)$ . (2.45)

This new algorithm for obtaining $B_{m+1}\left(z\right)$ from $B_m\left(z\right)$ and $B_m\left(0\right)$ is very easy to perform and may be utilized to establish Table of Bernoulli polynomials.

For examples:

$B_0\left(x\right)=1$

$B_1\left(x\right)=x-{\frac{x^2}{2}|}_{x=1}=x-\frac{1}{2}$

$B_2\left(x\right)=2\left(\frac{x^2}{2}-\frac{x}{2}\right)-2{\left(\frac{x^3}{6}-\frac{x^2}{4}\right)|}_{x=1}=x^2-x+\frac{1}{6}$

$B_3\left(x\right)=3\left(\frac{x^3}{3}-\frac{x^2}{2}+\frac{x}{6}\right)-3\left(\frac{1}{12}-\frac{1}{6}+\frac{1}{12}\right)=x^3-\frac{3}{2}{x}^2+\frac{x}{2}$

$B_4\left(x\right)=4\left(\frac{x^4}{4}-\frac{x^3}{2}+\frac{x^2}{4}\right)-4\left(\frac{1}{20}-\frac{1}{8}+\frac{1}{12}\right)=x^4-2x^3+x^2-\frac{1}{30}$

$B_5\left(x\right)=5\left(\frac{x^5}{5}-\frac{x^4}{2}+\frac{x^3}{3}-\frac{x}{30}\right)=x^5-\frac{5}{2}{x}^4+\frac{5}{3}{x}^3-\frac{1}{6}x$

$\begin{array}{c}{B}_6\left(x\right)=6\left(\frac{x^6}{6}-\frac{x^5}{2}+\frac{5x^4}{12}-\frac{x^2}{12}\right)+6\left(\frac{1}{6.7}-\frac{1}{2.4}+\frac{5}{12.4}-\frac{1}{12.3}\right)\\ =x^6-3x^5+\frac{5x^4}{2}-\frac{x^2}{2}+\frac{1}{42}\end{array}$

$B_7\left(x\right)=7\left(\frac{x^7}{7}-\frac{x^6}{2}+\frac{x^5}{2}-\frac{x^3}{6}+\frac{x}{42}\right)+0=x^7-\frac{7}{2}{x}^6+\frac{7x^5}{2}-\frac{7x^3}{6}+\frac{x}{6}$

$\begin{array}{c}{B}_8\left(x\right)=8\left(\frac{x^8}{8}-\frac{x^7}{2}+\frac{7x^6}{12}-\frac{7x^4}{24}+\frac{x^2}{12}\right)-8\left(\frac{1}{8.9}-\frac{1}{2.8}+\frac{7}{8.6}-\frac{7}{24.5}+\frac{1}{12.3}\right)\\ =x^8-4x^7+\frac{14}{3}{x}^6-\frac{7}{3}{x}^4+\frac{2}{3}{x}^2-\frac{1}{30}\end{array}$

$\begin{array}{c}{B}_9\left(x\right)=9\left(\frac{x^9}{9}-\frac{x^8}{2}+\frac{2x^7}{3}-\frac{7x^5}{15}+\frac{2x^3}{9}-\frac{x}{30}\right)+0\\ =x^9-\frac{9x^8}{2}+6x^7-\frac{21x^5}{5}+2x^3-\frac{3x}{10}\end{array}$

This method for establishing a table of Bernoulli polynomials is extremely easier if we utilize the list of fifty Bernoulli numbers $B_m\left(0\right)$ conscientiously established by Coen [2] . For examples

$B_{10}\left(x\right)=x^{10}-5x^9+\frac{15}{2}{x}^8-7x^6+5x^4-\frac{3}{2}{x}^2+\frac{5}{66}$

$B_{11}\left(x\right)=x^{11}-\frac{11}{2}{x}^{10}+\frac{55}{6}{x}^9-11x^7+11x^5-\frac{11}{2}{x}^3+\frac{5}{6}x$

$B_{12}\left(x\right)=x^{12}-6x^{11}+11x^{10}-\frac{33}{2}{x}^8+22x^6-\frac{33}{2}{x}^4+5x^2-\frac{691}{2730}$ (2.46)

2.6. Bernoulli Polynomials and Euler Zeta Function

From the Hurwitz formula

$\frac{1}{k!}{B}_k\left(z\right)=-\frac{1}{{\left(2i\pi \right)}^k}{\displaystyle \underset{n\in Z,n\ne 0}{\overset{\infty }{\sum }}\frac{1}{n^k}{\text{e}}^{i2\pi nz}}0\le z\le 1$

we get the Euler zeta function one may find references in Coen [2] and Raugh [3]

$\zeta \left(2m\right)={\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\frac{1}{k^{2m}}}={\left(-1\right)}^{m+1}\frac{1}{\left(2m\right)!2}{\left(2\pi \right)}^{2m}{B}_{2m}\left(0\right)$ (2.47)

as so as

${\left(2\pi \right)}^{2m}={\left(-1\right)}^{m+1}\frac{\left(2m\right)!2}{B_{2m}\left(z\right)}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\frac{1}{k^{2m}}}\cos 2\pi kz$ (2.48)

${\left(2\pi \right)}^{2m+1}={\left(-1\right)}^{m+1}\frac{\left(2m+1\right)!2}{B_{2m+1}\left(z\right)}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\frac{1}{k^{2m+1}}}\sin 2\pi kz$ (2.49)

Moreover, by taking $z=0,\frac{1}{8},\frac{1}{6},\frac{1}{4},\frac{1}{3},\frac{1}{2},\frac{2}{3},1$ in these formulae we get the known property

$B_{2m+1}\left(0\right)=B_{2m+1}\left(1\right)=-\left(2m+1\right)!{\displaystyle \underset{\begin{array}{l}k\in Z\\ k\ne 0\end{array}}{\overset{\infty }{\sum }}{\left(\frac{1}{2i\pi k}\right)}^{2m+1}}=0\text{for}m>0$ (2.50)

and the powers of pi.

For examples

${\left(2\pi \right)}^{2m}={\left(-1\right)}^{m+1}\frac{\left(2m\right)!2}{B_{2m}\left(1/2\right)}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\frac{{\left(-1\right)}^k}{k^{2m}}}$ (2.51)

${\pi }^{2m}={\left(-1\right)}^{m+1}\frac{\left(2m\right)!}{2^{2m-1}\left(1-2^{2m-1}\right)B_{2m}}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}{\left(-1\right)}^k\frac{1}{k^{2m}}}$ (2.52)

${\left(2\pi \right)}^{2m}={\left(-1\right)}^{m+1}\frac{\left(2m\right)!2}{B_{2m}\left(1/6\right)}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\cos \left(\frac{k\pi }{3}\right)\frac{1}{k^{2m}}}$ (2.53)

${\left(2\pi \right)}^{2m+1}={\left(-1\right)}^{m+1}\frac{\left(2m+1\right)!2}{B_{2m+1}\left(1/4\right)}{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\sin \left(\frac{k\pi }{2}\right)\frac{1}{k^{2m+1}}}$ (2.54)

and

$\frac{\pi }{4}=\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots$

${\left(\frac{\pi }{4}\right)}^2=1-\frac{1}{2^2}+\frac{1}{3^2}-\frac{1}{4^2}+\cdots$

${\pi }^2=36\left(\frac{1}{2}\frac{1}{1^2}-\frac{1}{2}\frac{1}{2^2}-\frac{1}{3^2}-\frac{1}{2}\frac{1}{4^2}+\frac{1}{2}\frac{1}{5^2}+\frac{1}{6^2}\right)+\left(\frac{1}{2}\frac{1}{7^2}-\frac{1}{2}\frac{1}{8^2}\right)+\cdots$

${\pi }^3=32\left(1-\frac{1}{3^3}+\frac{1}{5^3}-\frac{1}{7^3}+\cdots \right)$

etc.

3. Bernoulli Numbers

3.1. Definition and Properties

In 1713, according to Jacob Bernoulli (1655-1705), was published the list of ten first sums of powers of entire numbers [3]

${\displaystyle \sum n^m}=1^m+2^m+\cdots +n^m$ (3.1)

in terms of the numbers $B_k$ which are conjectured to be the same for all m

${\displaystyle \sum n^m}=\frac{1}{m+1}{\displaystyle \underset{k=0}{\overset{m}{\sum }}{\left(-1\right)}^k\left(\begin{array}{c}m+1\\ k\end{array}\right)B_k{n}^{m+1-k}}.$ (3.2)

Afterward, the $B_k$ were baptized Bernoulli numbers.

By comparison of the relation coming from (3.2)

$\begin{array}{c}{\partial }_n{\displaystyle \sum n^m}={\displaystyle \underset{k=0}{\overset{m}{\sum }}{\left(-1\right)}^k\frac{m!}{\left(m-k\right)!k!}{B}_k{n}^{m-k}}=:{\left(B-n\right)}^m\\ =B_0{n}^m-m{B}_1{n}^{m-1}+\frac{m\left(m-1\right)}{2}{B}_2{n}^{m-2}+\cdots +{\left(-1\right)}^m{B}_m\end{array}$ (3.3)

with the formula coming from (2.16), (2.17)

$\begin{array}{l}{\partial }_n\left(1^m+\cdots +n^{m-1}\right)+{\partial }_n{n}^m=B_m\left(n\right)+{\partial }_n{n}^m=:{\left(B\left(0\right)+n\right)}^m+m{n}^{m-1}\\ =B_0\left(0\right)n^m+m\left(1+B_1\left(0\right)\right)n^{m-1}+\left(\begin{array}{c}m\\ 2\end{array}\right)B_2\left(0\right)n^{m-2}+\cdots +B_m\left(0\right)\end{array}$ (3.4)

we get, combining with (2.20),

$B_0=B_0(\; 0\; )$

$B_1=-B_1\left(0\right)-1=-\frac{1}{2}$

$B_{2m}=B_{2m}(\; 0\; )$

$B_{2m+1}=-B_{2m+1}\left(0\right)=B_{2m+1}\left(1\right)=\frac{1}{2}{\delta }_{m0}$

$B_m=B_m\left(0\right)$ (3.5)

i.e.

“The Bernoulli numbers $B_m$ are equal to the values at origin of the Bernoulli polynomial $B_m\left(z\right)$ ”.

3.2. Obtaining Bernoulli Numbers

The above formula (3.5) and the recurrence formula for Bernoulli polynomials (2.43) corresponding to $z=0$

$0^m=\frac{1}{m+1}+{\displaystyle \underset{k=1}{\overset{m}{\sum }}\left(\begin{array}{c}m\\ k-1\end{array}\right)}\frac{B_k\left(0\right)}{k}$ (3.6)

lead to that for Bernoulli numbers

$\frac{1}{m+1}{B}_0+\left(\begin{array}{c}m\\ 0\end{array}\right)\frac{B_1}{1}+\left(\begin{array}{c}m\\ 1\end{array}\right)\frac{B_2}{2}+\left(\begin{array}{c}m\\ 2\end{array}\right)\frac{B_3}{3}+\cdots +\left(\begin{array}{c}m\\ m-1\end{array}\right)\frac{B_m}{m}=0,m>0$ (3.7)

which, knowing $B_0=B_m\left(0\right)=1$ , gives $B_1,B_2,B_4,\cdots ,B_m$ according to following Table 1.

This matrix equation may be resolved by doing linear combinations over lines from the second one in order to replace them with lines containing only some non-zero numbers.

Table 1. Matrix equation for calculating B_m.

For instance, for calculating successively $\left\{B_0,B_1,B_2,B_4,B_6,\cdots ,B_{18}\right\}$ we may utilize the matrix equation (Table 2).

We remark that the last line of this matrix has replaced $\left\{\left(\begin{array}{c}19\\ i\end{array}\right),i=0,1,2,4,6,8,10,12,14,16,18\right\}$ .

The results are

$B_0=1$ , $B_0+2B_1=0$ , $B_1+3B_2=0$ , $B_2+5B_4=0$ , $-B_2+7B_6=0$

$\frac{9}{5}{B}_2+9B_8=0$ . $5B_2+11B_{10}=0$ , $3B_1-\frac{61}{105}{B}_2+13B_{12}=0$ ,

$35B_1+15B_{14}=0$ , $240B_1+17B_3+17B_{16}=0$

$2052B_1-775B_4+19B_{18}=0=\frac{43867}{798}+B_{18}$ (3.8)

Another method, maybe more interesting, for establishing table of Bernoulli numbers is obtained from the formula (2.27). It is

$\left(-1-2m\right)B_{2m}=\left(\begin{array}{c}2m\\ 2\end{array}\right)B_{2m-2}{B}_2+\left(\begin{array}{c}2m\\ 4\end{array}\right)B_{2m-4}{B}_4+\cdots +\left(\begin{array}{c}2m\\ 2\end{array}\right)B_2{B}_{2m-2},m>1$

or, symbolically,

$\left(1-m\right)B_m=:{\left(B-B\right)}^m$ (3.9)

For examples

$\left(1-2\right)B_2=:{\left(B-B\right)}^2=2B_0{B}_2-2B_1{B}_1$

$-4B_3=:-3B_1{B}_2+3B_2{B}_1=0$

$-5B_4=\left(\begin{array}{c}4\\ 2\end{array}\right)B_2{B}_2=6B_2{B}_2\Rightarrow B_4=\frac{-1}{30}$

$-7B_6=2\left(\begin{array}{c}6\\ 4\end{array}\right)B_4{B}_2=30B_4{B}_2=-\frac{1}{6}\Rightarrow B_6=\frac{1}{42}$

$-9B_8=8\times 7\times B_6{B}_2+\left(\begin{array}{c}8\\ 4\end{array}\right)B_4{B}_4\Rightarrow B_8=-\frac{1}{30}$

$-11B_{10}=10\times 9\times B_8{B}_2+2\left(\begin{array}{c}10\\ 6\end{array}\right)B_6{B}_4\Rightarrow B_{10}=\frac{5}{66}$

Table 2. Simplified matrix equation for calculating B_m.

$-13B_{12}=12\times 11\times B_{10}{B}_2+2\left(\begin{array}{c}12\\ 8\end{array}\right)B_8{B}_4+\left(\begin{array}{c}12\\ 6\end{array}\right)B_6{B}_6\Rightarrow B_{12}=-\frac{691}{2730}$

$-15B_{14}=14\times 13\times B_{12}{B}_2+2\left(\begin{array}{c}14\\ 10\end{array}\right)B_{10}{B}_4+2\left(\begin{array}{c}14\\ 8\end{array}\right)B_8{B}_6\Rightarrow B_{14}=\frac{7}{6}$

$-17B_{16}=16\times 15\times B_{14}{B}_2+2\left(\begin{array}{c}16\\ 10\end{array}\right)B_{12}{B}_4+\left(\begin{array}{c}16\\ 8\end{array}\right)B_8{B}_8\Rightarrow B_{16}=-\frac{3617}{510}$

$-19B_{18}=18\times 17\times B_{16}{B}_2+2\left(\begin{array}{c}18\\ 14\end{array}\right)B_{14}{B}_4+2\left(\begin{array}{c}18\\ 12\end{array}\right)B_{12}{B}_6+2\left(\begin{array}{c}18\\ 10\end{array}\right)B_{10}{B}_8$

$B_{18}=\frac{43867}{798}$ , $B_{20}=-\frac{174611}{330}$ , etc. (3.10)

We see that $B_{18}$ is a sum over only four terms $B_{16}{B}_2,B_{14}{B}_4,B_{12}{B}_6,B_{10}{B}_8$ ; $B_{20}$ is over five, $B_{40}$ over ten, $B_{50}$ over twelve terms.

3.3. Obtaining Bernoulli Polynomials and Power Sums from Bernoulli Numbers

From the formula (2.15)

$B_m\left(z+a\right)=:{\left(B\left(a\right)+z\right)}^m$

we get the symbolic Lucas formula

$B_m\left(z\right)=:{\left(B+z\right)}^m$ (3.11)

for calculating Bernoulli polynomials $B_m\left(z\right)$ from the set of Bernoulli numbers.

For examples

$B_1\left(z\right)=:\left(B+z\right)=B_{0}z+B_1{z}^0=z-\frac{1}{2}$

$B_2\left(z\right)=:{\left(B+z\right)}^2=B_0{z}^2+2B_{1}z+B_2=z^2-z+\frac{1}{6}$

$B_3\left(z\right)=:{\left(B+z\right)}^3=B_0{z}^3+3B_1{z}^2+3B_{2}z+B_3=z^3-\frac{3}{2}{z}^2+\frac{1}{2}z$

As for the power sums $S_m\left(n\right)$ we begin by calculating the formula coming from (2.17)

${\partial }_n{S}_m\left(n\right)=B_m\left(n\right)=:{\left(B+n\right)}^m$ (3.12)

then take the primitives of both members.

For examples

${\partial }_n{S}_1\left(n\right)=B_1\left(n\right)\Rightarrow S_1\left(n\right)=0+1+\cdots +\left(n-1\right)={\displaystyle \int B_1}\left(n\right)={\displaystyle \int \left(n-\frac{1}{2}\right)}=\frac{n^2}{2}-\frac{n}{2}$

$S_2\left(n\right)=0^2+1^2+\cdots +{\left(n-1\right)}^2={\displaystyle \int B_2}\left(n\right)={\displaystyle \int \left(n^2-n+\frac{1}{6}\right)}=\frac{n^3}{3}-\frac{n^2}{2}+\frac{n}{6}$

$S_3\left(n\right)={\displaystyle \int B_3\left(n\right)}={\displaystyle \int \left(n^3-\frac{3}{2}{n}^2+\frac{n}{2}\right)}=\frac{n^4}{4}-\frac{n^3}{2}+\frac{n^2}{4}=\frac{n^2{\left(n-1\right)}^2}{4}$

3.4. Bernoulli Numbers and the Euler-McLaurin Formula

From the formula for expansion of derivable functions into series of Bernoulli polynomials

$f\left(z\right)={\displaystyle {\int }_0^{1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\left[f^{\left(k\right)}\left(1\right)-f^{\left(k\right)}\left(0\right)\right]}\frac{1}{\left(k+1\right)!}{B}_{k+1}\left(z\right)$ which leads, for periodic functions $B_m\left(z\right)$ identical to $B_m\left(z\right)$ in the interval $\left(0,1\right)$ , to

$f\left(z+m\right)={\displaystyle {\int }_m^{m+1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\left[f^{\left(k\right)}\left(m+1\right)-f^{\left(k\right)}\left(m\right)\right]}\frac{1}{\left(k+1\right)!}{B}_{k+1}\left(z\right)$ (3.13)

we get the formula

$\begin{array}{c}f\left(m\right)={\displaystyle {\int }_m^{m+1}f\left(z\right)\text{d}z}+{\displaystyle \underset{k=0}{\overset{r-1}{\sum }}\left[f^{\left(k\right)}\left(m+1\right)-f^{\left(k\right)}\left(m\right)\right]\frac{B_{k+1}\left(0\right)}{\left(k+1\right)!}}\\ +\frac{{\left(-1\right)}^{r+1}}{r!}{\displaystyle {\int }_m^{m+1}{f}^{\left(r\right)}\left(z\right)B_r\left(z\right)\text{d}z}\end{array}$ (3.14)

analogue to the Euler-McLaurin formula one may find in [11]

For example, with $f\left(z\right)=z^3$ , $B_1\left(0\right)=-\frac{1}{2}$ , $B_2\left(0\right)=\frac{1}{6}$ , $B_3\left(0\right)=0$ it is verified that

$2^3={\displaystyle {\int }_2^3{z}^3\text{d}z}-\frac{1}{2}\left(3^3-2^3\right)+3\left(3^2-2^2\right)\frac{1}{12}+6\left(3^1-2^1\right)\frac{0}{3}=\frac{81-16}{4}-\frac{19}{2}+\frac{15}{12}=8$

4. Obtaining Powers Sums of Real and Complex Numbers

4.1. From Power Sums of Integers

From the definition of the power sums on real and complex numbers

$S_m\left(z,n\right)=z^m+{\left(z+1\right)}^m+\cdots +{\left(z+\left(n-1\right)\right)}^m$ (4.1)

we get, by utilizing the translation operator ${\text{e}}^{{\partial }_z}$ mentioned in (2.4),

$S_m\left(z,n\right)=\left(1+{\text{e}}^{{\partial }_z}+\cdots +{\text{e}}^{\left(n-1\right){\partial }_z}\right)z^m$ (4.2)

and the formula for sums of geometric progressions, the compact formula

$S_m\left(z,n\right)=\frac{{\text{e}}^{n{\partial }_z}-1}{{\text{e}}^{{\partial }_z}-1}{z}^m$ (4.3)

From (4.3) and the fact that

${\partial }_{z+y}f\left(z+y\right)={\partial }_{z}f\left(z+y\right)={\partial }_{y}f\left(z+y\right)$ (4.4)

we get the symbolic formula

$S_m\left(z+y,n\right)=\frac{{\text{e}}^{n{\partial }_{z+y}}-1}{{\text{e}}^{{\partial }_{z+y}}-1}{\left(z+y\right)}^m=\frac{{\text{e}}^{n{\partial }_y}-1}{{\text{e}}^{{\partial }_y}-1}{\left(z+y\right)}^m=:{\left(z+S\left(y,n\right)\right)}^m$

leading to the very interesting new formula given powers sums of complex numbers from powers sums of integers

$S_m\left(z,n\right)={\left(S\left(n\right)+z\right)}^m$ (4.5)

where the undefined symbol $S^k\left(n\right)$ is to be replaced with the power sums on integers (2.17)

$S_k\left(n\right)=0^k+1^k+\cdots +{\left(n-1\right)}^k={\displaystyle \int B_k\left(n\right)}$ , $0^0=1$ (4.6)

Another way, more shortly, to obtain (4.5) is by remarking that

$\left(z+n\right)={\text{e}}^{z{\partial }_n}(\; n\; )$

so that

$S_m\left(z,n\right)={\text{e}}^{z{\partial }_n}{S}_m\left(n\right)={\displaystyle \underset{k=0}{\overset{m}{\sum }}\frac{z^k}{k!}{\partial }_n^k}{S}_m\left(n\right)={\displaystyle \underset{k=0}{\overset{m}{\sum }}\left(\begin{array}{c}m\\ k\end{array}\right)z^k}{S}_m{}_{-k}\left(n\right)={\left(S\left(n\right)+z\right)}^m$

For examples

$S_1\left(z,n\right)=S_1\left(n\right)z^0+S_0\left(n\right)z^1=\frac{n\left(n-1\right)}{2}+nz$

$S_2\left(z,n\right)=:S_2\left(n\right)+2S_1\left(n\right)z+S_0\left(n\right)z^2=\left(\frac{n^3}{3}-\frac{n^2}{2}+\frac{n}{6}\right)+\left(n^2-n\right)z+n{z}^2$

$S_3\left(z,n\right)=:{\left(S\left(n\right)+z\right)}^3=:{\displaystyle \int B_3\left(n\right)}+3z{\displaystyle \int B_2\left(n\right)}+3z^2{\displaystyle \int B_1\left(n\right)}+n{z}^3$

4.2. From Bernoulli Polynomials

Now, because n may go until infinity, ${\partial }_n$ is well defined so that

${\partial }_n{S}_m\left(z,n\right)={\partial }_n\frac{{\text{e}}^{n{\partial }_z}-1}{{\text{e}}^{{\partial }_z}-1}{z}^m=\frac{{\text{e}}^{n{\partial }_z}{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{z}^m=B_m\left(z+n\right)$ (4.7)

On the other hand, from (2.18)

${\partial }_z{S}_m\left(z,n\right)=\left({\text{e}}^{n{\partial }_z}-1\right)\frac{{\partial }_z}{{\text{e}}^{{\partial }_z}-1}{z}^m=B_m\left(z+n\right)-B_m\left(z\right)$ (4.8)

so that we obtain the following beautiful important formula

$\left({\partial }_n-{\partial }_z\right)S_m\left(z,n\right)=B_m\left(z\right)$ (4.9)

as so as the historic Jacobi conjectured formula

${\partial }_n{S}_m\left(n\right)=B_m\left(n\right)$ (4.10)

Formula (4.9) leads to the formula giving $S_m\left(z,n\right)$ directly from $B_m(\; z\; )$

$\begin{array}{c}{S}_m\left(z,n\right)=\frac{1}{{\partial }_n-{\partial }_z}{B}_m\left(z\right)=\left(1+\frac{{\partial }_z}{{\partial }_n}+{\left(\frac{{\partial }_z}{{\partial }_n}\right)}^2+\cdots \right)n{B}_m\left(z\right)\\ =\frac{n}{1!}{B}_m\left(z\right)+\frac{n^2}{2!}{B^{\prime }}_m\left(z\right)+\cdots +\frac{n^{m+1}}{\left(m+1\right)!}{B}_m^{\left(m\right)}\left(z\right)\end{array}$ (4.11)

i.e., to the algorithm saying that

$S_m\left(z,n\right)$ is equal to $n{B}_m\left(z\right)$ plus $\frac{n^2}{2!}{B^{\prime }}_m\left(z\right)$ and so all until $\frac{n^{m+1}}{\left(m+1\right)!}{B}_m^{\left(m\right)}(\; z\; )$

For examples

$S_1\left(z,n\right)=n{B}_1\left(z\right)+\frac{n^2}{2!}{B}_0\left(z\right)=nz-\frac{n}{2}+\frac{n^2}{2}$

$S_2\left(z,n\right)=n{B}_2\left(z\right)+\frac{n^2}{2!}2B_1\left(z\right)+\frac{n^3}{3!}2B_0\left(z\right)=n\left(z^2-z+\frac{1}{6}\right)+n^2\left(z-\frac{1}{2}\right)+\frac{n^3}{3}$

$S_3\left(z,n\right)=n{B}_3\left(z\right)+\frac{n^2}{2!}3B_2\left(z\right)+n^3{B}_1\left(z\right)+\frac{n^4}{4}$

In particular, we get the recurrence relation between Bernoulli polynomials given by Roman [8]

$S_m\left(z,1\right)=z^m=B_m\left(z\right)+\frac{m}{2!}{B}_{m-1}\left(z\right)+\cdots +B_1\left(z\right)+\frac{1}{m+1}{B}_0\left(z\right)$ (4.12)

and the well-known ancient formula of Bernoulli (1713)

$S_m\left(n\right)=\frac{n}{1!}{B}_m+\frac{n^2}{2!}m{B}_{m-1}+\cdots +\frac{n^m}{m!}{B}_1+\frac{n^{m+1}}{m+1}{B}_0$ (4.13)

Lastly, because of (4.10)

$n^m=S_m\left(n+1\right)-S_m\left(n\right)={\displaystyle {\int }_n^{n+1}{B}_m\left(n\right)\text{d}n}$

we get

$z^m=\left({\text{e}}^{{\partial }_z}-1\right)S_m\left(z\right)={\displaystyle {\int }_z^{z+1}{B}_m\left(n\right)\text{d}n}$

$\frac{{\text{e}}^{{\partial }_z}-1}{{\partial }_z}{B}_m\left(z\right)={\displaystyle {\int }_z^{z+1}{B}_m\left(n\right)\text{d}n}$

and, by expanding functions into Bernoulli series, the formula found in Wikipedia

$\frac{{\text{e}}^{{\partial }_z}-1}{{\partial }_z}f\left(z\right)={\displaystyle {\int }_z^{z+1}f\left(n\right)\text{d}n}=\left(1+\frac{{\partial }_z^2}{2!}+\frac{{\partial }_z^3}{3!}+\cdots \right)f\left(z\right)$ (4.14)

We resuming the herein-before results of calculations in following Tables (Tables 3-5).

Table 3. Obtaining $B_m\left(z\right)$ and $S_m\left(n\right)$ from $B_m$ .

Table 4. Obtaining $S_m\left(z,n\right)$ from $B_m\left(z\right)$ .

Table 5. Obtaining $S_m\left(z,n\right)$ from $S_m\left(n\right)$ .

5. The Faulhaber Formulae on Power Sums of Complex Numbers

5.1. Powers Sums of Odd Order

Although the problems of powers sums and Faulhaber conjecture were treated by many authors for examples by Radermacher [12] , by Tsao in (2008) [13] , by Chen, Fu, Zhang in (2009) [14] , etc., nevertheless we would like to present hereafter one new approach about the problems.

In $S_m\left(z,n\right)$ let us replace the arguments z and n by

$Z=z\left(z-1\right)$ and $\lambda =S_1\left(z,n\right)=B_1\left(z\right)n+\frac{n^2}{2}$ (5.1)

Because

$\frac{\text{d}Z}{\text{d}z}=2B_1\left(z\right),\frac{\text{d}Z}{\text{d}n}=0,\frac{\text{d}\lambda }{\text{d}n}=B_1\left(z\right)+n,\frac{\text{d}\lambda }{\text{d}z}=n$ (5.2)

and consequently

${\partial }_n\equiv \frac{\text{d}Z}{\text{d}n}{\partial }_Z+\frac{\text{d}\lambda }{\text{d}n}{\partial }_{\lambda }=\left(B_1\left(z\right)+n\right){\partial }_{\lambda }$

${\partial }_z\equiv \frac{\text{d}Z}{\text{d}z}{\partial }_Z+\frac{\text{d}\lambda }{\text{d}z}{\partial }_{\lambda }=2B_1\left(z\right){\partial }_Z+n{\partial }_{\lambda }$

${\partial }_n-{\partial }_z\equiv B_1\left(z\right)\left({\partial }_{\lambda }-2{\partial }_Z\right)$ (5.3)

we have, regarding (4.9),

$\left({\partial }_{\lambda }-2{\partial }_Z\right)S_m\left(z,n\right)=B_1^{-1}\left(z\right)\left({\partial }_n-{\partial }_z\right)S_m\left(z,n\right)=B_1^{-1}\left(z\right)B_m\left(z\right)$ (5.4)

and the form of the formula for general power sums

$2S_m\left(z,n\right)=\left(2\lambda +\frac{{\left(2\lambda \right)}^2}{2!}{\partial }_Z+\cdots +\frac{{\left(2\lambda \right)}^m}{m!}{\partial }_Z^{m-1}\right)B_1^{-1}\left(z\right)B_m\left(z\right)$ (5.5)

that may be calculated by the following considerations.

From the property

$\begin{array}{c}\left(2k+2\right)B_{2k+1}\left(z\right)={\partial }_z{B}_{2k+2}\left(z\right)=\frac{\text{d}Z}{\text{d}z}{\partial }_Z{B}_{2k+2}\left(z\right)\\ =\left(2z-1\right){\partial }_Z{B}_{2k+2}\left(z\right)=2B_1\left(z\right){\partial }_Z{B}_{2k+2}\left(z\right)\end{array}$ (5.6)

we get, for utilization in (5.5),

$B_1^{-1}\left(z\right)B_{2k+1}\left(z\right)=\frac{1}{k+1}{\partial }_Z{B}_{2k+2}\left(z\right)$ (5.7)

and, finally,

$\begin{array}{c}2S_{2k+1}\left(z,n\right)=\left(2\lambda +\frac{{\left(2\lambda \right)}^2}{2!}{\partial }_Z+\cdots +\frac{{\left(2\lambda \right)}^{2k+1}}{\left(2k+1\right)!}{\partial }_Z^{2k}\right)\frac{1}{k+1}{\partial }_Z{B}_{2k+2}\left(z\right)\\ ={\displaystyle \underset{j=1}{\overset{2k+1}{\sum }}\frac{{\left(2\lambda \right)}^j}{j!}{\partial }_Z^{j-1}}\left({\partial }_Z\frac{B_{2k+2}\left(z\right)}{k+1}\right)\end{array}$ (5.8)

All the problem is reduced to the calculations of ${\partial }_Z{B}_{2m+2}\left(z\right)$ in function of Z which are not so difficult.

For examples:

$\begin{array}{c}2B_1^{-1}\left(z\right)B_3\left(z\right)={\partial }_Z{B}_4\left(z\right)={\partial }_Z\left(z^4-2z^3+z^2-\frac{1}{30}\right)\\ ={\partial }_Z{\left(z^2-z\right)}^2={\partial }_Z{Z}^2=2Z\end{array}$

$2S_3\left(z,n\right)=\left(2\lambda \right)Z+\frac{{\left(2\lambda \right)}^2}{2!}=2\lambda Z+2{\lambda }^2$ , $\lambda =S_1\left(z,n\right)$

$2S_3\left(n\right)=\frac{{\left(2\lambda \right)}^2}{2!}=2{\lambda }^2$

$\begin{array}{c}3B_1^{-1}\left(z\right)B_5\left(z\right)={\partial }_Z{B}_6\left(z\right)={\partial }_Z\left(z^6-3z^5+\frac{5}{2}{z}^4-\frac{1}{2}{z}^2+\frac{1}{42}\right)\\ ={\partial }_Z\left(z^6-3z^5+3z^4-z^3\right)-\frac{1}{2}{z}^4+z^3-\frac{1}{2}{z}^2\\ ={\partial }_Z\left(Z^3-\frac{1}{2}{Z}^2\right)=3Z^2-Z\end{array}$

$\begin{array}{c}{S}_5\left(z,n\right)=\frac{1}{2}Z\left(Z-\frac{1}{3}\right)\frac{2\lambda }{1!}+\frac{1}{2}\left(2Z-\frac{1}{3}\right)\frac{{\left(2\lambda \right)}^2}{2!}+\frac{1}{2}\left(2\right)\frac{{\left(2\lambda \right)}^3}{3!}\\ =\left(Z^2-\frac{1}{3}Z\right)\lambda +\left(2Z-\frac{1}{3}\right){\lambda }^2+4\frac{{\lambda }^3}{3}\end{array}$

$S_5\left(n\right)=-\frac{1}{3}{S}_1^2\left(n\right)+\frac{4}{3}{S}_1^3(\; n\; )$

$\begin{array}{c}4B_1{}^{-1}\left(z\right)B_7\left(z\right)={\partial }_Z{B}_8\left(z\right)={\partial }_Z\left(z^8-4z^7+\frac{14}{3}{z}^6-\frac{7}{3}{z}^4+\frac{2}{3}{z}^2\right)\\ ={\partial }_Z\left(Z^4-\frac{4}{3}{Z}^3+\frac{2}{3}{Z}^2\right)=4Z^3-4Z^2+\frac{4}{3}Z\end{array}$

$2S_7\left(z,n\right)=\left(Z^3-Z^2+\frac{1}{3}Z\right)\frac{2\lambda }{1!}+\left(3Z^2-2Z+\frac{1}{3}\right)\frac{{\left(2\lambda \right)}^2}{2!}+\left(6Z-2\right)\frac{{\left(2\lambda \right)}^3}{3!}$

and so all.

As corollary of the calculations of ${\partial }_Z{B}_{2k+2}\left(z\right)$ we may state that

“All $B_{2k}\left(z\right)$ and all $B_1^{-1}\left(z\right)B_{2k+1}\left(z\right)$ are polynomials of order k in Z”.

5.2. Faulhaber Formula for Even Power Sums $S_{2k}\left(z,n\right)$

By differentiating both members of (5.7) and remarking that ${\partial }_{z}Z=Z^{\prime }=\left(2z-1\right)$ , ${\lambda }^{\prime }=n$ we obtain the formula giving $S_{2m}\left(z,n\right)$

$\begin{array}{c}\left(2m+1\right)S_{2m}\left(z,n\right)={\displaystyle \underset{k=1}{\overset{2m+1}{\sum }}\frac{{\left(2\lambda \right)}^k}{k!}{B}_1\left(z\right){\partial }_Z^k}\left(\frac{1}{m+1}{\partial }_Z{B}_{2m+2}\left(z\right)\right)\\ +{\displaystyle \underset{k=1}{\overset{2m+1}{\sum }}\frac{{\left(2\lambda \right)}^{k-1}2n}{\left(k-1\right)!}{\partial }_Z^{k-1}}\left(\frac{1}{m+1}{\partial }_Z{B}_{2m+2}\left(z\right)\right)\end{array}$ (5.9)

For examples

$S_3\left(z,n\right)=Z\lambda +{\lambda }^2$

$S_5\left(z,n\right)=\left(\frac{1}{2}{Z}^2-\frac{1}{6}Z\right)\frac{2\lambda }{1!}+\left(Z-\frac{1}{6}\right)\frac{{\left(2\lambda \right)}^2}{2!}+\frac{{\left(2\lambda \right)}^3}{3!}$

$5S_4\left(z,n\right)=Z^{\prime }\left(Z-\frac{1}{6}\right)\frac{2\lambda }{1!}+\left(\frac{1}{2}{Z}^2-\frac{1}{6}Z\right)2n+Z^{\prime }\frac{{\left(2\lambda \right)}^2}{2!}+\left(Z-\frac{1}{6}\right)4\lambda n+4{\lambda }^{2}n\text{}$

The arrangement into polynomials with respect to ( $2\lambda$ ) is immediate.

Remarks and Conclusions

We subjectively think that this work is a real and effective contribution to the knowledge of Bernoulli polynomials, Bernoulli numbers and Sums of powers of entire and complex numbers, as indicated in Introduction.

The main particularity of this work is the use of the translation or shift operator ${\text{e}}^{a{\partial }_z}$ that is curiously let apart by quasi all authors although this is seen to be very useful and easy to utilize.

By the utilization of many new properties on $B_m\left(z\right)$ such as

$S_m\left(n\right)={\displaystyle \int B_m(\; n\; )}$

$B_{m+1}\left(z\right)=\left(m+1\right){\displaystyle {\int }_0^z{B}_m\left(x\right)\text{d}x}-\left(m+1\right){\displaystyle {\int }_0^1\text{d}z}{\displaystyle {\int }_0^z{B}_m\left(x\right)\text{d}x}$

$B_m\left(z+y\right)=:m\left(z+y\right)B_{m-1}\left(z+y\right)-{\left(B\left(z\right)+B\left(y+1\right)\right)}^m$

$\left(z+n\right)={\text{e}}^{z{\partial }_n}(\; n\; )$

we easily get the new key formulae

$S_m\left(z,n\right)=:{\left(S\left(n\right)+z\right)}^m$ together with $\left({\partial }_n-{\partial }_z\right)S_m\left(z,n\right)=B_m\left(z\right)$ for obtaining $S_m\left(z,n\right)$ .

We find also the miraculous symbolic formula for calculating rapidly the Bernoulli numbers

$\left(1-m\right)B_m=:{\left(B-B\right)}^m$

which together with the Lucas symbolic formula

$B_m\left(z\right)=:{\left(B+z\right)}^m$

give easily $B_m\left(z\right)$ .

Afterward by a change of arguments from z into $Z=z\left(z-1\right)$ and n into $\lambda =S_1\left(z,n\right)$ we get the relation $\left({\partial }_n-{\partial }_z\right)=B_1\left(z\right)\left({\partial }_{2\lambda }-{\partial }_Z\right)$ which together with the proof that $B_1^{-1}\left(z\right)B_{2k+1}\left(z\right)$ and $B_{2k}\left(z\right)$ are polynomials in Z gives simply rise to the Faulhaber form of $S_m\left(z,n\right)$ .

Operator calculus, which is very different from Heaviside operational calculus thus merits to be known. Moreover, it has a solid foundation and many interesting applications in the domains of Special functions, Differential equations, Fourier and other transforms, quantum mechanics [10] .

Acknowledgements

Because there are thousands of works about Bernoulli numbers and polynomials during centuries we surely have omitted to cite many references, we apologize for this and would like to receive comments from researchers in order to correct this work. Before, we thank you very much for that.

He thanks the assistant managers of Applied Mathematics, SCIRP very much noticeably Ms. Icy Yin and Ms. Nancy Ho for helping him in the procedure of publication of his former and recent works.

Finally, he reiterates his gratitude towards his lovely wife for all the cares she devoted for him during this work and also during his life.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Bernoulli, J. (1713) Ars Conjectandi.
[2]	Coen, L.E.S. (1996) Sums of Powers and the Bernoulli Numbers. Master’s Thesis, Eastern Illinois University, Charleston. http://thekeep.eiu.edu/theses/1896
[3]	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/
[4]	Beardon A.F. (1996) Sums of Powers of Integers. The American Mathematical Monthly, 103, 201-213 https://doi.org/10.1080/00029890.1996.12004725
[5]	Dilcher, K., Shula, L. and Slavutskii, I. (1991) Bernoulli Numbers Bibliography (1713-1990). Queen’s Papers in Pure and Applied Mathematics 87. Updated On-Line Version. http://www.mathstat.dal.ca/~dilcher/bernoulli.html
[6]	Si, D.T. (2019) The Powers Sums, Bernoulli Numbers, Bernoulli Polynomials Rethinked. Applied Mathematics, 10, 100-112. https://doi.org/10.4236/am.2019.103009
[7]	Roman, S. (1984) The Umbral Calculus, Pure and Applied Mathematics. Academic Press, New York, p. 111.
[8]	Castellanos D. (1988) The Ubiquitous π. Mathematics Magazine, 61, 67-98. https://doi.org/10.1080/0025570X.1988.11977350
[9]	Raabe, J.L. (1851) Zurückführung einiger Summen und bestimmten Integrale auf die Jacob-Bernoullische Function. Journal für die reine und angewandte Mathematik, 1851, 348-367. https://doi.org/10.1515/crll.1851.42.348
[10]	Do, T.S. (2016) Operator Calculus: Edification and Utilization. LAP LAMBERT Academic Publishing, Saarbrücken.
[11]	Weisstein, E.W. (n.d.) Bernoulli Polynomial. https://mathworld.wolfram.com/Euler-%20MaclaurinIntegrationFormulas.html
[12]	Rademacher, H. (1973) Topics in Analytic Number Theory (Die Grundlehren der Mathematischen Wissenschaften). Springer-Verlag, New York. https://doi.org/10.1007/978-3-642-80615-5
[13]	Tsao, H. (2008) Explicit Polynomial Expressions for Sums of Powers of an Arithmetic Progression. The Mathematical Gazette, 9, 87-92. https://doi.org/10.1017/S0025557200182610
[14]	Chen, W.Y.C., Fu, A.M. and Zhang, I.F. (2009) Faulhaber’s Theorem on Power Sums. Discrete Mathematics, 309, 2974-2981. https://doi.org/10.1016/j.disc.2008.07.027