How to calculate this double sum? This occurred in a probability problem where I have to calculate the invariant $c$ which equals to $N$ divided by the following double summation:
$$\sum_{n=0}^{N} \sum_{k=0}^N |k-n|$$
 A: Exploiting the symmetry of $|n-k|$, it is equivalent to compute
$$2\sum_{n=0}^N\sum_{k=0}^n(n-k)=2\sum_{n=0}^N\sum_{k=0}^nk=\sum_{n=0}^Nk(k+1)=\frac{N(N+1)(N+2)}3.$$
(The sum of triangular numbers is a pyramidal number.)
A: We have
$$S_1=\sum_{k=0}^{n-1}(n-k)=n^2-\frac{n(n-1)}{2}$$
and 
$$S_2=\sum_{k=n+1}^N(k-n)=\frac{(N-n)(N+n+1)}{2}-(N-n)n$$
so 
$$\sum_{k=0}^N\vert n-k\vert=S_1+S_2$$
and to get the final sum use also the equality
$$\sum_{n=0}^N n^2=\frac{N(N+1)(2N+1)}{6}$$
A: This sum can be performed in a simple, straightforward manner, as follows.
$$\begin{align} \sum_{n=0}^N \sum_{k=0}^N |k-n| &= \sum_{n=0}^N \left [\sum_{k=0}^{n-1} (n-k) + \sum_{k=n+1}^N (k-n) \right ] \\ &= \sum_{n=0}^N \left (n^2 - \frac12 n (n-1) + \frac12 \left [N (N+1) - n (n+1) \right ] -n (N-n) \right )\\ &= \sum_{n=0}^N \left (\frac12 n (n+1) + \frac12 N (N+1) - \frac12 n (N+1) - n (N-n) \right)\\ &= \frac12 N (N+1)^2 - N \sum_{n=0}^N n + \sum_{n=0}^N n^2 \end{align} $$
Thus,

$$\sum_{n=0}^N \sum_{k=0}^N |k-n| = \frac13 N (N+1) (N+2) $$

ADDENDUM
You can also exploit the symmetry.  Imagine arranging the $(N+1)^2$ points in a lattice.  We can simply sum along lines parallel to the line $k=n$.  The first lines have a diff of $1$, the second lines $2$, and so on.  The sum then looks like double the following:
$$N + 2 (N-1) + 3 (N-2) + \cdots + N [N - (N-1)] = N \frac12 N (N+1) - \sum_{k=1}^N k(k-1)$$
I leave it to the OP to verify that the double the above expression is equal to that derived above.
A: Non conventional approach:
The double sum is a cubic polynomial in $N$. Indeed, the summation domain $[0,N]^2$ can be decomposed in two triangular parts with summing of $n-k$ or $k-n$. As these terms are linear, the summation on $k$ yields a quadratic polynomial in $n$, and the summation on $n$ yields a cubic polynomial in $N$.
Then the solution is the Lagrangian polynomial defined by the computed points $S_0=0,S_1=2,S_2=8,S_3=20$.
