# Compilation of proofs for the summation of natural squares and cubes

I want to know different proofs for the following formulas,

$$\sum_{i=1}^n{i^2} = \frac{(n)(n+1)(2n+1)}{6}$$

$$\sum_{i=1}^n{i^3} = \frac{n^2(n+1)^2}{2^2}$$

Please do not mark this as duplicate, since what I specifically want is to be exposed to a variety of proofs using different techniques (I did not find such a compilation anywhere on the net)

I am only familiar with two proofs, one which uses expansion of $(x+1)^2 - x^2$ and $(x+1)^3 - x^3$ and the other which uses induction. I have provided link for the induction proof in a self-answer.

I am particularly interested in proof without words and proofs which use a unrelated mathematical concept (higher level math upto class 12 level is acceptable).

Also,

(1) Don't think I am being rude or anything, it is out of genuine interest that I am asking this question.

(2) Someone marked this as a duplicate of Methods to compute $\sum_{k=1}^nk^p$ without Faulhaber's formula

My question is different in three ways:

(i) I want to focus only on these two summation and not the general case,

(ii) Hence, it follows that the proofs which I am looking for a simpler than the ones provided in that link and are simpler (using images, pictures or high school algebra). What I want is to study new proofs. I believe it is a good practice when learning math to so this.

(iii) Since the proofs in the link are given for the general case, they are complicated and I am finding it hard to understand them. If someone is able to use the same method to the two cases in my question, then it would probably become much simpler and easier to digest.

Appendix

Calculus

Basic Binomial Expansion

Coordinate Geometry

Algebra (upto what 18 year olds learn)

Taylor Series Expansions

Geometry (18 year old level)

basically.............math which 18 year old's learn on Earth.

If you want to err, then err on the higher math side:)

1. By Newton series
2. By Sterling Numbers
3. By Induction
4. From the book Generatingfunctionology by Herbert Wilf
5. By generalizing the following pattern \begin{align} &\ \,4\cdot5\cdot6\cdot7\\=&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6+4\cdot5\cdot6\cdot7)\\-&(0\cdot1\cdot2\cdot3+1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6)\\ =&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot4+3\cdot4\cdot5\cdot4+4\cdot5\cdot6\cdot4)\\ \end{align}

6. By Lagrangian Interpolation

7. By Formal Differentiation
8. By the Euler-Maclaurin Summation Formula
9. By Assuming that the expression is a polynomial of degree $2$.
10. A Proof Without Words for the cube case
11. By integrating and assuming a error term.
12. SimplyBeatifulArt's Personal Approach
• There is a standard proof using the fact that $(a+1)^3-a^3=(a+1)^2+(a+1)+1$. The left side telescopes (all but the first and last terms cancel out), and the right side is your desired sum + other things that you can calculate (like the the sum of numbers from $1$ to $n$). – fierydemon Jul 24 '17 at 16:55
• Possible duplicate of Methods to compute $\sum_{k=1}^nk^p$ without Faulhaber's formula because it has a very large variety of proofs of these formulas. – Simply Beautiful Art Jul 24 '17 at 16:57
• Please, we're told it isn't a duplicate: if the OP asks it, it's something special. – Professor Vector Jul 24 '17 at 16:59
• @ProfessorVector Just saying, it does contain a very large variety of proofs of these things... perhaps the OP just didn't see it :-) – Simply Beautiful Art Jul 24 '17 at 17:00
• I think so myself, I've seen at least a half dozen methods to prove such things (probably more). And if somebody needs a compilation: they can just make one. – Professor Vector Jul 24 '17 at 17:03

For the sum of cubes there is a nice illustration in the wikipedia article Squared Triangular Number Illustrating Nicomachus's theorem.

$$\sum_{i=1}^ni^3=\left(\sum_{i=1}^ni\right)^2$$ By integrating and assuming an error term

Supposing known: $$\sum_{i=1}^{n}{1}=n \tag{*}$$ and $$\sum_{i=1}^{n}{i}=\frac{n(n+1)}{2}\tag{**}$$

Let $\sigma$ be: $$\sigma(n)=\sum_{i=1}^{n}{i^2}$$ This is similar to the integral: $$\Sigma(n)=\int_{0}^{n}{x^2\,dx}$$ up to an error $\epsilon(n)$. Therefore: $$\sigma(n) = \Sigma(n) +\epsilon(n) \tag{1}$$ This error may be obtained by the fact that the curve $x^2$ is always below of the actual value $i^2$, hence applying $(*)$ and $(**)$: $$\epsilon(n) = \sum_{i=1}^{n}{\left[i^2-\int_{i-1}^{i}{x^2\,dx}\right]}=\sum_{i=1}^{n}{\left[i-\frac{1}{3}\right]}=\frac{n(n+1)}{2}-\frac{n}{3}$$ Integrating $\Sigma(n)$ and adding them up in $(1)$ the desired result is reached: $$\sigma(n) = \frac{n^3}{3}+\frac{n(n+1)}{2}-\frac{n}{3}=\frac{n(2n+1)(n+1)}{6}$$

The one for your second required proof involves the previous result, and it is done in a similar fashion. Try it!

• Recommend looking up the Euler-Maclaurin formula, which is one way to derive Faulhaber's formula, as robjohn describes. – Simply Beautiful Art Jul 24 '17 at 21:15
• @SimplyBeautifulArt I don't know what you want me to do. Please would you mind being more specific?. Btw, the OP asked a simple way to calculate these sums (please do not misinterpret my comment). I am not a matematician, and therefore I think that sometimes, the simplicity and the easy-to-follow reasonings are better than those generalized formulas that are definitely not understandeable by anyone who has not a Phd in mathematics. In my humble opinion. – HBR Jul 24 '17 at 21:44
• Of course. I was simply noting that your $\epsilon(n)$ has another form, which is given by the Euler-Maclaurin formula, an extension of trapezoidal rule, and it estimates the difference between an integral and corresponding Riemann sum with $\Delta x=1$ – Simply Beautiful Art Jul 24 '17 at 23:48

Let $f(x,p)$ be a polynomial of $x$ with degree $p+1$ such that for $n\in\mathbb N$, we have

$$f(n,p)=\sum_{k=1}^nk^p$$

This function satisfies the recursive relation:

$$f(x,p)=a_px+\int_0^xf(t,p-1)~\mathrm dt$$

where

$$a_p=1-p\int_0^1f(t,p-1)~\mathrm dt$$

For example, if know that

$$f(x,1)=\frac{x(x+1)}2$$

Then it follows that $a_1=\frac12$ and

$$f(x,2)=\frac12x+\int_0^x\frac{t(t+1)}2~\mathrm dt=\frac{x(x+1)(2x+1)}6$$

which is the formula for $\sum k^2$.

As described in this answer, we have

$$\sum_{k=1}^nk^p=\sum_{k=1}^p{p\brace k}\frac{(n+1)^{\underline{k+1}}}{k+1}$$

Where ${p\brace k}$ is a stirling number of the second kind and $a^{\underline b}=a(a-1)(a-2)\dots(a-b+1)$ is a falling factorial. For your cases,

$$\sum_{k=1}^nk^2={2\brace1}\frac{(n+1)(n)}2+{2\brace2}\frac{(n+1)(n)(n-1)}3\\=\frac{n(n+1)}2+\frac{n(n+1)(n-1)}3\\=\frac{n(n+1)(2n+1)}6$$

$$\sum_{k=1}^nk^3={3\brace1}\frac{(n+1)(n)}2+{3\brace2}\frac{(n+1)(n)(n-1)}3+{3\brace3}\frac{(n+1)(n)(n-1)(n-2)}4\\=\frac{n(n+1)}2+n(n+1)(n-1)+\frac{n(n+1)(n-1)(n-2)}4\\=\left[\frac{n(n+1)}2\right]^2$$

Of course, for normal purposes, just use the table of values to evaluate basic stirling numbers.

I'm surprised that no one has pointed out the elementary Newton series, which I'll describe below.

The forward difference operator $Δ$ defined as: $\def\nn{\mathbb{N}} \def\zz{\mathbb{Z}} \def\lfrac#1#2{{\large\frac{#1}{#2}}} \def\lbinom#1#2{{\large\binom{#1}{#2}}}$

$Δ = ( \text{function$f$on$\zz$} \mapsto ( \text{int$n$} \mapsto f(n+1) - f(n) ) )$

Namely, for any function $f$ on $\zz$ and $n \in \zz$, we have $Δ(f)(n) = f(n+1) - f(n)$.

If you think of the functions as sequences (infinite in both directions), then taking the forward difference means replacing each term with the value of the next term minus itself. For example if you repeatedly take the forward difference of the sequence of cubes:

...,-27,-8,-1, 0, 1, 8,27,...
..., 19, 7, 1, 1, 7,19,37,...
...,-12,-6, 0, 6,12,18,24,...
...,  6, 6, 6, 6, 6, 6, 6,...
...,  0, 0, 0, 0, 0, 0, 0,...
...,  0, 0, 0, 0, 0, 0, 0,...


This powerful abstraction makes it easy to get a lot of things. For instance, the numbers obtained here can be easily used to obtain the general formula for sum of cubes, as you desire.

General method for indefinite summation

The key is that:

$Δ\left( \text{int$n$} \mapsto \lbinom{n}{k+1} \right) = \left( \text{int$n$} \mapsto \lbinom{n}{k} \right)$ for any $k \in \zz$.

This is to be expected because it follows directly from extending Pascal's triangle naturally, namely if we define $\lbinom{n}{k}$ by the recurrence:

$\lbinom{n}{0} = 1$ for any $n \in \zz$.

$\lbinom{0,k+1}{0} = 0$ for any $k \in \nn$.

$\lbinom{n+1}{k+1} = \lbinom{n}{k+1} + \lbinom{n}{k}$ for any $k \in \nn$ and $n \in \zz$.

Now consider any function $f$ on $\zz$ such that $f(n) = \sum_{k=0}^d a_k \lbinom{n}{k}$ for any $n \in \zz$. Then we have for any $m \in \nn_{\le d}$:

$Δ^m(f)(n) = \sum_{k=0}^{d-m} a_{k+m} \lbinom{n}{k}$ for any $n \in \zz$.

And hence:

$Δ^m(f)(0) = a_m$.

Which immediately gives Newton's series:

$f(n) = \sum_{k=0}^d Δ^k(f)(0) \lbinom{n}{k}$ for any $n \in \zz$.

From a high-level perspective, this is the discrete version of the Taylor series.

This works for any polynomial function $f$ on $\zz$, since $D^k(f)$ is the zero function once $k$ is larger than the degree of $f$, so we can use it to immediately find the series for $(\text{int n} \mapsto n^3)$, and then just take the anti-difference by shifting the coefficients of the series the other way. The undetermined constant that appears will drop out once we perform a definite sum like if we want the sum of the first $m$ cubes.

Sum of $p$ powers

For example if we want $\sum_{k=1}^{n-1} k^3$ we first find the iterated forward differences of the sequence of cubes $( n^3 )_{n \in \zz}$:

..., 0, 1, 8,27,64,...
..., 1, 7,19,37,...
..., 6,12,18,...
..., 6, 6,...
..., 0,...


So we immediately get $n^3 = 0 \binom{n}{0} + 1 \binom{n}{1} + 6 \binom{n}{2} + 6 \binom{n}{3}$ and hence $\sum_{k=0}^{n-1} k^3 = 0 \binom{n}{1} + 1 \binom{n}{2} + 6 \binom{n}{3} + 6 \binom{n}{4} = \lfrac{n(n-1)}{2} \Big( 1 + \lfrac{6(n-2)}{3} \big( 1 + \lfrac{n-3}{4} \big) \Big) = \Big( \lfrac{n(n-1)}{2} \Big)^2$.

Computation efficiency

This is far more efficient than the method you mentioned in your question (namely by taking summation on both sides of $(n+1)^3-n^3 = 3n^2+3n+1$ and telescoping) because the Newton series is easy to manipulate and easy to compute. For sum of $p$-powers we only need $O(p^2)$ arithmetic operations to find the forward-differences and then $O(p^2)$ more to simplify the series form into a standard polynomial form. In contrast, the other method requires $O(p^3)$ arithmetic operations.

• Curios to know whether the computation time is for humans or computers? – Agile_Eagle Jul 26 '17 at 13:14
• @PrashantGokhale: For both. It was in terms of arithmetic operations, which are additions, subtractions or multiplications of integers. This is reasonable if you do not need arbitrary precision (such as if 64-bit arithmetic is enough for you). If you need exact answers, then you need more accurate accounting, since arithmetic on larger integers obviously takes more time. The computational efficiency of using the Newton series over other methods is still there and should be even more pronounced, but I didn't really analyze it. – user21820 Jul 26 '17 at 14:45
• @PrashantGokhale: If you really want to compare hardness for humans, just try out all the methods for some suitably large $p$, say $p = 7$, and time yourself. I can easily bet that you will find this method a breeze but may even give up on the other methods. – user21820 Jul 26 '17 at 14:47

By formal differentiation:

We express the sum as a polynomial in $x$ and evaluate it at $x=1$ (taking the limit, by L'Hospital).

$$S_0(n)=\sum_{k=1}^nx^k=\frac{x^{n+1}-x}{x-1}\to n.$$

$$S_1(n)=\sum_{k=1}^nkx^k=xS_0'(n)=\frac{nx^{n+2}-(n+1)x^{n+1}+x}{(x-1)^2}\to\frac{n^2+n}2.$$

$$S_2(n)=\sum_{k=1}^nk^2x^k=xS_1'(n)=\frac{n^2x^{n+2}-(2n^2+2n-1)x^{n+1}+(n^2+2n)x^n-x-1}{(x-1)^3}\to\frac{2n^3+3n^2+n}6.$$

$$S_3(n)=\sum_{k=1}^nk^3x^k=xS_2'(n)=\frac{(n^3-n^2)x^{n+3}-(3n^3-5n+2)x^{n+1}+(3n^3+3n^2-5n-1)x^n-(n^3+2n^2)x^{n-1}+2n+4}{(x-1)^4}\to\frac{n^4+2n^3+n^2}4.$$

The computation can be somewhat eased by the change of variable $x=t+1$, such that you keep working with polynomials (and the Binomial theorem).

$$S_0(n)=\sum_{k=0}^{n-1}x^k=\frac{x^n-1}{x-1}=\frac{(t+1)^n-1}t=\sum_{i=1}^n\binom nit^{i-1}.$$

$$S_1(n)=xS_0'(n)=\sum_{k=0}^{n-1}kx^k=(t+1)\sum_{i=2}^n\binom ni(i-1)t^{i-2}=\sum_{i=2}^n\binom ni(i-1)(t^{i-1}+t^{i-2})\to\binom n1.$$

$$S_2(n)=xS_1'(n)=\sum_{k=0}^{n-1}kx^k=(t+1)\sum_{i=2}^n\binom ni(i-1)((i-1)t^{i-2}+(i-2)t^{i-3})\\ =\sum_{i=2}^n\binom ni(i-1)((i-1)t^{i-1}+(2i-3)t^{i-2}+(i-2)t^{i-3})\to\binom n2+2\binom n3.$$

• Pretty tedious ! – Yves Daoust Jul 25 '17 at 6:13

There is a neat proof for the sum of squares using generating functions in the book Generatingfunctionology by Herbert Wilf, pages 35 - 36. First take the generating function $\sum_{n=0}^N x^n = \frac{x^{N+1}-1}{x-1}$. Then, we convert it into the generating function $\sum_{n=0}^N n^2x^n$ by applying $(xD)^2$ to both sides (this is a standard technique that you can learn about in the book): \begin{align*} (xD)^2\sum_{n=0}^N x^n &= (xD)\sum_{n=0}^N x\cdot nx^{n-1} = (xD)\sum_{n=0}^N nx^n = \sum_{n=0}^N x\cdot n^2x^{n-1} = \sum_{n=0}^N n^2x^n = (xD)^2\frac{x^{N+1}-1}{x-1} \end{align*} Now evaluate both sides at $x=1$ (take the limit as $x\to1$) and massage the resulting expression until you get what you want! Bashing out the derivative and limit on the right hand side takes forever though, I didn't have the guts to do it out.

Of course, you could do the same for $n^3$ but that would be even more bashing than for $n^2$.

• +1 if you can edit the answer and make it a little more digestible – Agile_Eagle Jul 24 '17 at 17:39
• What parts would you like me to expand on? – Tai Jul 24 '17 at 17:39
• Here's an even better approach: Let $x=e^t$ so that we get $$\sum_{n=0}^Ne^{nt}$$Now there's no need to multiply by $x$, as we may simply use$$\sum_{n=0}^Nn^ke^{nt}=\frac{d^k}{dt^k}\sum_{n=0}^Ne^{nt}$$See [Albert Steppi's answer](math.stackexchange.com/a/2040960/272831). Of course, this approach has the advantage of Faa Di Bruno's formula, but for small values of $k$, this solution is simple enough. – Simply Beautiful Art Jul 24 '17 at 17:40
• @TaisukeYasuda 1. The formatting 2. The language 3. The part of bashing and all 4. and if possible the standard technique – Agile_Eagle Jul 24 '17 at 17:42
• @PrashantGokhale can you be more specific? How would you like the formatting and language? And what parts/how much of the bashing would you like to see? – Tai Jul 24 '17 at 17:50

\begin{align} &\ \,4\cdot5\cdot6\cdot7\\=&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6+4\cdot5\cdot6\cdot7)\\-&(0\cdot1\cdot2\cdot3+1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6)\\ =&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot4+3\cdot4\cdot5\cdot4+4\cdot5\cdot6\cdot4)\\ \end{align}

By generalizing the pattern, this is enough to establish

$$n^{(m)}=\sum_{k=1}^n k^{(m)}-\sum_{k=0}^{n-1} k^{(m)}=(m+1)\sum_{k=1}^nk^{(m-1)},$$ where $a^{(b)}$ denotes the raising factorial.

From this,

$$n=\sum_{k=1}^n 1=S_0(n)$$

$$n(n+1)=2\sum_{k=1}^n k=S_1(n)$$ $$n(n+1)(n+2)=3\sum_{k=1}^n k(k+1)=S_2(n)$$ $$n(n+1)(n+2)(n+3)=4\sum_{k=1}^n k(k+1)(k+2)=S_3(n).$$

They are the "linear", triangular, tetrahedral, "hypertetrahedral" numbers.

Then

$$\sum_{k=1}^n 1=S_0(n)=n$$

$$\sum_{k=1}^n k=\frac12S_1(n)=\frac{n^2+n}2$$

$$\sum_{k=1}^n k^2=\frac13S_2(n)-\sum_{k=1}^n k=\frac{2n^3+3n^2+n}6$$

$$\sum_{k=1}^n k^3=\frac14S_3(n)-3\sum_{k=1}^n k^2-2\sum_{k=1}^n k=\frac{n^4+2n^3+n^2}4$$

One can continue in a systematic way.

• Aka hockey stick formula? – Simply Beautiful Art Jul 24 '17 at 21:11
• Nice, nice this looks good, hehe:) – Agile_Eagle Jul 25 '17 at 4:37

While it might be considered overkill, use of the Euler-Maclaurin Summation Formula facilitates evaluation of the sum $\sum_{i=1}^n i^m$ for any integer $m\ge 0$.

Invoking the EMSF for $f(i)=i^2$ yields

\begin{align} \sum_{i=1}^n i^2&=\int_0^n x^2\,dx +\frac12 (n^2)+\frac1{12}\left(2n\right)\\\\ &=\frac{n^3}{3}+\frac12 n^2 +\frac16 n\\\\ &=\frac n6 (2n+1)(n+1) \end{align}

and invoking the EMSF for $f(i)=i^3$ yields

\begin{align} \sum_{i=1}^n i^3&=\int_0^n x^3\,dx +\frac12 (n^3)+\frac1{12}\left(3n^2\right)\\\\ &=\frac{n^4}{4}+\frac{n^3}{2}+\frac{n^2}{4}\\\\ &=\frac{n^2}4 (n+1)^2 \end{align}

as expected!

• I've deduced sth. similar in a post below yours. Btw, I did't know that this was indeed the Euler-Mclaurin Summation Formula – HBR Jul 24 '17 at 18:46
• @HBR Your methodology is not the same as the EMSF, which involves the Bernoulli numbers. And your solution did not provide an answer for the sum of cubes of integers. ;-)) – Mark Viola Jul 24 '17 at 18:54
• That's the reason I wrote "similar". It seems interesting how by yours they can be directly obtained without the need of previous results like in my approach. Btw. I left the cubes for the reader... their sum is obtained in a similar manner. Greetings! – HBR Jul 24 '17 at 19:00

By Lagrangian interpolation.

The expression of the sum of the cubes must be a quartic polynomial, because its first order difference is cubic. Furthermore, it has no constant term because the sum of no numbers is zero.

Hence, the average of the cubes is a cubic polynomial, by the four points

$$(1,1),\left(2,\frac{1+2^3}2\right),\left(3,\frac{1+2^3+3^3}3\right),\left(4,\frac{1+2^3+3^3+4^4}4\right).$$

The resquested polynomial is the Lagrangian interpolant.

Average of ones (constant, by $(1,1)$):

$$\overline S_0(n)=1$$

Average of naturals (linear, by $(1,1)$ and $(2,\frac32)$):

$$\overline S_1(n)=\frac{\frac32-1}{2-1}(n-1)+1=\frac{n+1}2$$

Average of squares (quadratic, by $(1,1)$, $(2,\frac52)$ and $(3,\frac{14}3)$):

$$\overline S_2(n)=\cdots=\frac{2n^2+3n+1}2$$

Average of cubes (cubic, by $(1,1)$, $(2,\frac92)$, $(3,\frac{26}3)$ and $(4,\frac{45}2)$):

$$\overline S_3(n)=\cdots=\frac{n^3+2n^2+n}4$$

An explicit formula can be obtained for any degree, by means of the Lagrangian interpolation formula for equidistant points.

Newtonian interpolation could be easier.

One technique is the Newton series, which permits termwise summation just as the Taylor series permits termwise integration. There is also the more powerful but less efficient method of indefinite summation, which is unsurprisingly the discrete version of integration by parts. $\def\zz{\mathbb{Z}} \def\lbinom#1#2{{\large\binom{#1}{#2}}}$

Let $R$ be the right-shift operator defined as:

$R = ( \text{function$f$on$\zz$} \mapsto ( \text{int$n$} \mapsto f(n+1) ) )$

Namely, for any function $f$ on $\zz$ and $n \in Z$, we have $R(f)(n) = f(n+1)$.

For convenience we also define the summation operator:

$Σ = ( \text{function$f$on$\zz$} \mapsto ( \text{int$n$} \mapsto \sum_{k=0}^{n-1} f(k) ) )$

Then we have the important property that $ΔΣ(f) = f$ for any function $f$ on $\zz$, analogous to the fundamental theorem of calculus.

Now we can derive summation by parts from the product rule in the same manner as for the continuous version. Consider any functions $f,g$ on $\zz$:

$Δ(f·g)(n) = f(n+1) g(n+1) - f(n) g(n) = f(n+1) Δ(g)(n) - Δ(f)(n) g(n)$ for any $n \in \zz$.

And hence the discrete product rule (with the usual pointwise sum and product of functions):

$Δ(f·g) = R(f)·Δ(g) + Δ(f)·g$.

Now by substituting $f$ with $Σ(f)$ and taking summation on both sides we get summation by parts:

$Δ(Σ(f)·g) = R(Σ(f))·Δ(g) + f·g$.

$Σ(f·g) = Σ(f)·g - Σ(R(Σ(f))·Δ(g)) + c$ for some constant function $c$ on $\zz$.

For convenience we let "$Δ_n E(n)$" denote "$E(n+1)-E(n)$" and "$Σ_n E(n)$" denote "$\sum_{k=0}^{n-1} E(k)$".

We can therefore compute $\sum_{k=0}^{n-1} k^3$ as follows:

$\sum_{k=0}^{n-1} k^3 = Σ_n\Big(\lbinom{n}{0}·n^3\Big)$   [Recall this definition of binomial coefficients.]

$\ = \lbinom{n}{1}·n^3 - Σ_n\Big(\lbinom{n+1}{1}·(3n^2+3n+1)\Big) + c$, where $c$ is some constant (independent of $n$),

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + Σ_n\Big(\lbinom{n+2}{2}·(6n+6)\Big) + c'$,

where $c'$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - Σ_n\Big(\lbinom{n+3}{3}·6\Big) + c''$,

where $c''$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - \lbinom{n+3}{4}·6 + c''$.

By substituting $n = 0$ we readily get $c'' = 0$ and hence simplifying gives the desired formula.

Using this we can also easily compute things like $\sum_{k=0}^{n-1} k^3 3^k$ by applying it three times, each time reducing the degree of the polynomial part. Incidentally, there is another purely discrete technique to compute this sum, namely subtracting it from $3 \sum_{k=1}^n k^3 3^k = n^3 3^{n+1} + \sum_{k=1}^n (k-1)^3 3^k$, to reduce the polynomial degree, and repeating three times. But interestingly this technique does not work for $\sum_{k=1}^n k^3$.

By indeterminate coefficients.

The sum of the cubes must be a quartic polynomial with no independent term, so

$$S_3(n)=an^4+bn^3+cn^2+dn.$$

With $n=1,2,3,4$,

$$\begin{cases}a+b+c+d&=1,\\16a+8b+4c+2d&=9,\\ 81a+27b+9c+3d&=36,\\ 256a+64b+16c+4d&=100.\\ \end{cases}$$

$$\begin{cases}14a+6b+2c&=7,\\ 78a+24b+6c&=33,\\ 252a+60b+12c&=96.\\ \end{cases}$$

$$\begin{cases} 36a+6b&=12,\\ 168a+24b&=54.\\ \end{cases}$$

$$\begin{cases} 24a&=6.\\ \end{cases}$$

Then $b=\dfrac12,c=\dfrac14,d=0$.

There is a faster way, using

$$S_3(n)-S_3(n-1)=n^3\\ =a(n^4-(n-1)^4)+b(n^3-(n-1)^3)+c(n^2-(n-1)2)+d(n-(n-1))\\ =4an^3+(-6a+3b)n^2+(4a-3b-2c)n+(-a+b-c+d).$$

Then by identifying the coefficients, you get an easy triangular system:

$$\begin{cases}-a+b-c+d&=0,\\4a-3b+2c&=0,\\-6a+3b&=0,\\4a&=1.\end{cases}$$

Notice that the system is easy to write, as the columns are the rows of Pascal's triangle with alternating signs and the last coefficient missing.

This is not the most elegant way to do it, but there is a direct visualization of a related sum which directly implies that first formula.

# Starting one step lower...

Consider the following visual proof for the sum of counting numbers: they arrange into a triangle which can be overlaid upon itself to form a square, duplicating the diagonal: The implication is that the sum of counting numbers, being visually represented as one of these triangles, must be $(n^2 + n)/2$ which can afterwards be simplified to $n(n+1)/2,$ and the key focus that I want to highlight is this double-counting of stuff on the diagonal.

This is the basic mechanism that we're going to use to directly visualize the sum of the first $m$ odd squares.

# Going from odds to general sums.

There is also a slightly more direct visualization also that the sum of the first $n$ odd numbers is the $n^\text{th}$ square number: to see this just color the top left corner of an $n\times n$ square green, then color the three cells beside it blue, then color the 5 cells beside those cells red, and so on... one immediately sees that each color highlights a subsequent odd number. They are odd, because they are made up of two equal sides plus the one square on the diagonal.

One can then use this sum of odd numbers to prove the same formula. Let me only consider the case where we end on an even number: sum up the counting numbers from $1$ to $2m$; this consists of the first $m$ odd numbers plus the first $m$ even numbers (if the "first" even number is regarded as 2, not 0). Subtracting 1 from each even number gives the corresponding odd number, so the sum is $2m^2 + m = n^2/2 + n/2.$ One can make a similar argument for odds or just add $n+1$ to this form, which is very similar to a proof by induction but lacks the inductive step.

# Extending this directly

The obvious 3D generalization of that first overlay diagram concerns the sum of odd squares, as we build a staggered pyramid with a $1\times1$ square for its topmost layer supported by a $3\times3$ square for its second layer supported by a $5\times5$ square for its third layer, and so on, up to a $(2n-1)\times(2n-1)$ square at the bottom. I don't have a great 3D visualization package for this sort of thing so the quickest I could mock up looks like this: We can use 6 of these staggered pyramids to complete a supercube with sides of length $2n-1$, if we can count the overlaid squares. We start by counting every unit subcube of the overall supercube once, $(2n-1)^3$, and this gets most of what we need. Looking at just one of the pyramids above we see that there is the "topmost" square (in blue) which needs to be counted 6 times total, beside it each of the edges of length $n-1$ (in gold) needs to be counted 3 times total, and flanked by those edges there is an "inner triangle" which needs to be counted 2 times total. Since we've already counted each of these once we need 5 more for the teal cube, 2 more of each gold cube, and one more of these grassy cubes forming the inner triangle on the outside. Of these, the inner triangles are the only ones that look daunting -- until you realize that they are a sum of the first $n-1$ odd numbers and therefore are just $(n-1)^2,$ so subtracting the one time we've already counted each of these squares we need to add $5 + 2\cdot(n-1)\cdot(\text{# of edges}) + 1\cdot(n-1)^2(\text{# of inner triangles})$ yielding the expression: \begin{align}6\sum_{k=1}^{n} (2k-1)^2 =& (2n-1)^3 + 12 (n - 1)^2 + 16 (n-1) + 5\\ =&~(2n-1)~2n~(2n+1).\end{align}Now if we have a sum of the squares from $1$ to $4m^2$ we find that this is a sum of the first $m$ odd squares plus the first $m$ even squares, but each even square is a subsequent odd number above an odd square: $3 + 7 + 11 + \dots + (4m-1).$ Subtracting $1$ from all $m$ of these and dividing by $2$ gives this as $2(1 + 3 + 5 + \dots + (2m-1)) + m = 2m^2 + m,$ so our final result is $$\sum_{k=1}^{2m} k^2 = \frac{(2m-1)~2m~(2m+1)}{3} + 2m^2 + m = \frac13 m (4m+1)(2m + 1),$$ from which the result follows.