calculus - 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.