How to get formula for sums of powers? Assuming I have Bernoulli numbers:
$B = [\frac{1}{1},\frac{1}{2},\frac{1}{6},\frac{0}{1},-\frac{1}{30}, \frac{0}{1}, \frac{1}{42}, ...]$
How can I get the coefficients of the sums of powers formulas?
For example the sum of squares is $(1/3)n^3 + (1/2)n^2 + 1/6n$
 A: The exponential generating function for the Bernoulli Numbers is
$$
\frac{x}{1-e^{-x}}\tag1
$$
Formally, Taylor's Theorem says
$$
f(x+a)=e^{aD}f(x)\tag2
$$
where
$$
\begin{align}
e^{aD}f(x)
&=\sum_{k=0}^\infty\frac{a^k}{k!}D^kf(x)\\
&=\sum_{k=0}^\infty\frac{a^k}{k!}f^{(k)}(x)\tag3
\end{align}
$$
and $D^{-1}$ is an anti-derivative.
For a polynomial, the right side of $(2)$ is finite.
Define
$$
F(m)=\sum_{k=1}^mf(k)\tag4
$$
Applying $(2)$ to $(4)$ yields
$$
\begin{align}
f(m)
&=F(m)-F(m-1)\\
&=\left(1-e^{-D}\right)\!F(m)\tag5
\end{align}
$$
We can formally invert $(5)$ to get
$$
\begin{align}
\frac{D}{1-e^{-D}}D^{-1}x^n
&=\sum_{k=0}^\infty\frac{B_k}{k!}D^{k-1}x^n+C\\
&=\sum_{k=0}^\infty\frac{B_k}{k!}\frac{n!}{(n-k+1)!}\,x^{n-k+1}+C\\
&=\frac1{n+1}\sum_{k=0}^{n+1}B_k\binom{n+1}{k}x^{n-k+1}+C\tag6
\end{align}
$$
Note that these formal operations are exact on polynomials. Furthermore, $C$ is the constant of integration from $D^{-1}x^n=\frac1{n+1}x^{n+1}+C$, which comes from the $k=0$ term.
Since the $k=n+1$ term in $(6)$ is a constant and $F(0)=0$, $C$ cancels the $k=n+1$ term, and we are left with Faulhaber's Formula:
$$
\sum_{k=1}^mk^n=\frac1{n+1}\sum_{k=0}^nB_k\binom{n+1}{k}m^{n-k+1}\tag7
$$
A: What you're looking for is Faulhaber's formula (Wikipedia link):
$$\sum_{k=1}^n k^p = {1 \over p+1} \sum_{j=0}^p (-1)^j{p+1 \choose j} B_j n^{p+1-j},\qquad \mbox{where}~B_1 = -\frac{1}{2}$$
expressing the sum of the first $n$ $p$th powers as a polynomial in $n$ whose coefficients involve the Bernoulli numbers and some binomial coefficients.
To get the example you listed,
$$\begin{align*}
\sum_{k=1}^nk^2&=\frac{1}{3}\sum_{j=0}^2(-1)^j\binom{3}{j}B_jn^{3-j}\\\\
&=\left(\frac{1}{3}(-1)^0\binom{3}{0}B_0\right)n^3+\left(\frac{1}{3}(-1)^1\binom{3}{1}B_1\right)n^2+\left(\frac{1}{3}(-1)^2\binom{3}{2}B_2\right)n\\\\
&=\left(\frac{1}{3}\cdot 1\cdot 1\cdot 1\right)n^3+\left(\frac{1}{3}\cdot (-1)\cdot 3\cdot \left(-\frac{1}{2}\right)\right)n^2+\left(\frac{1}{3}\cdot 1\cdot 3\cdot \frac{1}{6}\right)n\\\\
&=\left(\frac{1}{3}\right)n^3+\left(\frac{1}{2}\right)n^2+\left(\frac{1}{6}\right)n
\end{align*}$$
A: Not as a good an answer but might be easier to program in some cases:
$$f_p(n)=\sum_{k=0}^nk^p$$
$$f_p(n)=\frac{(n+1)^{p+1}-(n+1)-\sum_{j=1}^{p-1} {{p+1} \choose j}f_j(n)}{p+1}$$
$p>0$.
$f_1(n)=\frac{n(n+1)}{2}.$
A: Stiring Numbers of the Second Kind is also good for building sum of powers formula. 
$x^n = \sum _{k=0} ^n S(n,k) (x)_k$
Example, for sum of x^2 formula, "integrate" the falling factorial form.
$$x^2 = (1)(x)_1 + (1)(x)_2$$
$$\sum_{k=0}^{n-1} x^2 
= \frac{(n)_2}{2}+\frac{(n)_3}{3}
= \frac{n^3}{3} - \frac{n^2}{2} + \frac{n}{6}$$
A: A vey-well rememberable scheme is perhaps the following.      
Use the vector of Bernoulli-numbers in the first row of a little scheme and multiply through each column (the result per column is in the last row of the table):
$$ \begin{bmatrix}\frac{1}{1}&\frac{1}{2}&\frac{1}{6}&\frac{0}{1}&-\frac{1}{30}& \frac{0}{1}& \frac{1}{42}& ... \\
*x^4 & *x^3&*x^2&*x& & & &  \\ 
/4 & /3& /2&/1& & & &  \\ 
*1 & *3& *3&*1& & & &  \\ \hline
1/4x^4 & 1/2x^3& 1/4x^2& 0x& & & & 
\end{bmatrix}$$
Sum of cubes is $1/4x^4 + 1/2x^3+ 1/4x^2$
$$ \begin{bmatrix}\frac{1}{1}&\frac{1}{2}&\frac{1}{6}&\frac{0}{1}&-\frac{1}{30}& \frac{0}{1}& \frac{1}{42}& ... \\
*x^5 & *x^4&*x^3&*x^2&*x & & &  \\ 
/5 & /4& /3&/2&/1 & & &  \\ 
*1 & *4& *6&*4&*1 & & &  \\ \hline
1/5x^5 & 1/2x^4& 1/3x^3& 0x^2&-1/30x 
\end{bmatrix}$$
sum-of-$4$'th powers is $ 1/5x^5 + 1/2x^4 + 1/3x^3  -1/30x $
