The sum of squares of the first $n$ natural numbers. My basic question is this: how to find the sum of squares of the first $n$ natural numbers?
My thoughts have led me to an interesting theorem: Faulhaber's formula. It is known that $$1^k+2^k+\ldots+n^k=P_{k+1}(n)$$ is a polynomial of degree $n$ $(k+1)$ (!). For my problem: $$1^2+2^2+\ldots+n^2=a+bn+cn^2+dn^3.$$ Further resolving an uncomplicated system of linear equations:
$$\left\{
\begin{aligned}
0=&a\\
1^2=&a+b\cdot1+c\cdot1^2+d\cdot1^3\\
1^2+2^2=&a+b\cdot2+c\cdot2^2+d\cdot2^3\\
1^2+2^2+3^2=&a+b\cdot3+c\cdot3^2+d\cdot3^3\\
\end{aligned}
\right.$$ 
Thus we get: $a=0,\,b=\frac16,\,c=\frac12,\,d=\frac13$, i.e$$P_3(n)=\frac{n(n+1)(2n+1)}{6}.$$
My further questions:
1) What are some ways to find the sum of the squares of the first n natural numbers?
2) (!) How to prove that the sum of $1^k+2^k+\ldots+n^k$ is a polynomial of degree $n$ $k+1$?
 A: An integer-valued polynomial is a polynomial with complex coefficients taking values in $\mathbb{Z}$ when all the variables take integer values.  For example, $\frac{x^2+3x}{2}$ and $\frac{13x^3+5xy^2}{6}$ are integer-valued polynomials.  Clearly, the set of integer-valued polynomials with variables $x_1,x_2,\ldots,x_n$ form a subring of $\mathbb{Q}\left[x_1,x_2,\ldots,x_n\right]$.  A result by Polya states that the ring of integer-valued polynomials in one variable $x$ is a free abelian group with basis elements $\binom{x}{k}=\frac{x(x-1)(x-2)\cdots(x-k+1)}{k!}$ for $k=0,1,2,\ldots$.
To answer your question, $x^k$ is an integer-valued polynomial.  Therefore, $x^k=\sum_{r=0}^k \,a_r\binom{x}{r}$ for some $a_0,a_1,\ldots,a_n\in\mathbb{Z}$ (obviously, $a_n\neq 0$).  Now, $\sum_{m=0}^n\,m^k=\sum_{m=0}^n\,\sum_{r=0}^k\,a_r\binom{m}{r}=\sum_{r=0}^k\,a_k\,\sum_{m=0}^n\,\binom{m}{r}$.  By the Hockey-Stick Identity (see http://www.artofproblemsolving.com/wiki/index.php/Combinatorial_identity#Hockey-Stick_Identity), $\sum_{m=0}^n\,m^k=\sum_{r=0}^k\,a_k\,\binom{n+1}{r+1}$.  Hence, $\sum_{m=0}^n\,m^k$ is a polynomial in $n$ of degree $k+1$, as the coefficient of $n^{k+1}$ is $\frac{a_k}{(k+1)!}\neq 0$.  (In fact, $a_k=k!$, so we know that $\sum_{m=0}^n\,m^k=\frac{n^{k+1}}{k+1}+\mathcal{O}\left(n^k\right)$.)
A: There is a simple combinatorial identity that can be helpful here:
$$
\sum_{m=0}^n \binom{m+k}{k} = \binom{n+k+1}{k+1}
$$
The right-hand side counts the number of subsets of $\{1,\ldots,n+k+1\}$ of size $k+1$. The left-hand side counts them by their maximal element $m+k+1$.
Simple linear algebra shows that every polynomial of degree $d$ is a linear combination of $\binom{n+e}{e}$ for $e \leq d$ (in a unique way). In particular, we can represent $n^d$ in this way, and so the combinatorial identity implies that $\sum_{m=0}^n m^d$ is a polynomial of degree $d+1$.
