# Bernoulli numbers in the general series-expansion formula for sums of powers?

We know the series-expansion for these "sum of powers":

$$\sum_{i=0}^n i = \frac{n(n+1)}{2} = \frac{1}{2}n^2+\frac{1}{2}n$$ $$\sum_{i=0}^n i^2 = \frac{n(n+1)(n+2)}{6} = \frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n$$ $$\sum_{i=0}^n i^3 = \frac{n^2(n+1)^2}{4} = \frac{1}{4}n^4+\frac{1}{2}n^3+\frac{1}{4}n^2$$ $$...$$ But can we come up with a general formula? $$\sum_{i=0}^n i^P =$$ "Compact form" = "Series form"
Using Mathematica, I generated sums of higher powers to look for overall patterns in the coefficients and the polynomials. I could not find anything significant in the "compact form" but I found several in the "series form." Take a look: Series expansion formula for sum of powers from P=1 to P=15 ;; Coefficients only
So what I found was that the "sum of powers" formula has this general equation $$S(n)=An^{P+1}+Bn^P+Cn^{P-1}+Dn^{P-2}+...+En$$. If we include the fact that some of the coefficients are $$0$$, then there are $$P+1$$ terms in this polynomial equation. I was tackling this problem with the assumption that $$P$$ was a natural number but it turns out that this general equation is valid for any real number $$P$$. Also, the general series expansion formula is equivalent to the series expansion for the generalized harmonic number (https://en.wikipedia.org/wiki/Harmonic_number).
Anyway, here are the patterns I found in the coefficients: $$A = \frac{1}{P+1}, B=\frac{1}{2}, C=\frac{P}{12}$$ $$D=0$$ and every second letter after D is $$0$$.
$$E=\frac{-P(P-1)(P-2)}{720}, G=\frac{P(P-1)(P-2)(P-3)(P-4)}{720},$$ $$I=\frac{-P(P-1)(P-2)(P-3)(P-4)(P-5)(P-6)}{1209600},...$$
And so on. It appears that for every coefficient letter, this is the numerator: $$\frac{P!(-1)^{0.5k+1.5}}{(P-k)!}$$ where k is the $$k^{th}$$ letter after B with k starting at 1 (and k is odd in this formula since every second letter after D is $$0$$). But where does the denominator come from? If you examine more coefficient letters, the denominator seems to be growing exponentially (factorial?), but that is all I could find. 
Update! I later found the generalized formula here: https://en.wikipedia.org/wiki/Faulhaber%27s_formula. This is what the denominator equals: $$\frac{(k+1)!}{|B(k+1)|}$$ Amazing how the Bernoulli numbers belong in this formula. That explains the alternating sums and the coefficients of $$0$$. So here is the complete equation: $$\sum_{i=0}^n i^P =\frac{1}{P+1}n^{P+1}+\frac{1}{2}n^P+\sum_{k=1}^{P-1} \frac{B(k+1)P!}{(P-k)!(k+1)!}n^{P-k}$$ This can be further condensed, assuming $$B(1)=+\frac{1}{2}$$: $$\sum_{i=0}^n i^P =\sum_{k=-1}^{P-1} \frac{B(k+1)P!}{(P-k)!(k+1)!}n^{P-k}$$

Interesting. The equation kind of resembles the binomial expansion formula (https://en.wikipedia.org/wiki/Binomial_theorem). But why do the Bernoulli numbers appear in this formula? I looked at several sites (https://en.wikipedia.org/wiki/Bernoulli_number, http://math.ucr.edu/~res/math153/s12/bernoulli-numbers.pdf, https://ncatlab.org/nlab/show/Bernoulli+number), but the explanations are not intuitive. From what I see it, "they just do" but I do not understand why. Let me know if there are any mistakes with my equations or if there is anything I am missing.

Consider the identity $$1+e^{x} +e^{2x}+\ldots+e^{nx}=\frac{e^{(n+1)x}-1}{e^x-1}= \frac{e^{(n+1)x}-1}{x}\frac{x}{e^x-1}$$ for power series and note that $$\frac{x}{e^x-1}=\sum_{k=0}^\infty\frac{B_k}{k!}x^k.$$