Finding a closed form for $\sum\limits_{\substack{0\le n\le N\\0\le m\le M}}\left|nM-Nm\right|$ I am trying to figure out if there is a closed form for the following sum:
$$
\sum_{\substack{0\le n\le N \\ 0\le m\le M}}\left|(N-n)(M-m)-nm\right|=\sum_{\substack{0\le n\le N \\ 0\le m\le M}}\left|nM-Nm\right|.
$$
Clearly, from symmetry, if we remove the absolute values, the sum evaulates to $0$.
Using the symmetry, I tried to evaluate the sum as
$$
2\cdot\sum_{\substack{0\le n\le N \\ 0\le m\le M \\ (N-n)(M-m)\gt nm}}\left((N-n)(M-m)-nm\right).
$$
However, that expression ended up containing sums of floor functions, for which I did not know a closed form.
Is there another way to obtain a closed form for the above sum?

Edit: If no closed form can be provided an efficient algorithm to compute the above sum given $N,M$ would also be good.

Edit 2: Since the time I have placed the bounty, I was able to find a closed form for the sum:

$$\frac{1}{6}\left[MN(2 M N + 3 (M + N + 1)) + M^2 + N^2-\gcd(M,N)^2\right].$$

So now I change the question to a challange: derive the above form from the sum. The prettiest derivation (if there are any) will get the bounty.
 A: For completeness, here's my own derivation:
First, given $d:=\gcd(M,N)$, $\,N':= N/d$, we have the following identity:
$$
\sum_{n=0}^{N}f\left(\left\{n\frac{M}{N}\right\}\right)
=f\left(0\right)+d\sum_{n=0}^{N'-1}f\left(\frac{n}{N'}\right)
$$
Removing the absolute value will give a sum that equals $0$. Hence, the positive part of the sum is equal to the negative part. Thus:
$$
\begin{equation}
\begin{split}
\sum_{n=0}^N\sum_{m=0}^M\left|nM-Nm\right| &= 2\sum_{n=0}^N\sum_{m=0}^{\left\lfloor n\frac{M}{N}\right\rfloor}\left(nM-Nm\right) \\
&=\sum_{n=0}^N\left(2\left(\left\lfloor n\frac{M}{N}\right\rfloor+1\right)nM-N\left(\left\lfloor n\frac{M}{N}\right\rfloor^2+\left\lfloor n\frac{M}{N}\right\rfloor\right)\right) \\
&=\sum_{n=0}^N\left(\left\lfloor n\frac{M}{N}\right\rfloor+1\right)\left(2nM-N\left\lfloor n\frac{M}{N}\right\rfloor\right) \\
&=\sum_{n=0}^N\left(n\frac{M}{N}-\left\{ n\frac{M}{N}\right\}+1\right)\left(nM+N\left\{ n\frac{M}{N}\right\}\right) \\
&=\sum_{n=0}^N\left(n^2\frac{M^2}{N}+nM-N\left\{ n\frac{M}{N}\right\}^2+N\left\{ n\frac{M}{N}\right\}\right) \\
&=\small{\frac{1}{6}\left((N+1)(2N+1)M^2+3N(N+1)M-Nd\frac{(N'-1)(2N'-1)}{N'}+3Nd(N'-1)\right)} \\
&=\small{\frac{1}{6}\left((N+1)(2N+1)M^2+3N(N+1)M-(N-d)(2N-d)+3N(N-d)\right)} \\
&=\frac{1}{6}\left(MN(2MN + 3(M+N+1)) +M^2 +N^2-d^2\right)
\end{split}
\end{equation}
$$
A: Hint:
The problem is governed by the solutions of $Mn=Nm$, the main diagonal of the rectangle. WLOG $M\ge N$ and let us assume for now that $M,N$ are relative primes, so that equality only occurs at the corners.
By symmetry, we just look below the diagonal and for a given $m$,
$$0\le n\le\left\lfloor{\frac{Nm}M}\right\rfloor=\left\lfloor{Qm}\right\rfloor.$$
Then, 
$$S:=\sum_{m=0}^{M-1}\sum_{n=0}^{\left\lfloor{Qm}\right\rfloor}(Nm-Mn)=M\sum_{m=0}^{M-1}\sum_{n=0}^{\left\lfloor{Qm}\right\rfloor}(Qm-n)\\
=M\sum_{m=0}^{M-1}\left(Qm-\frac12\left\lfloor{Qm}\right\rfloor\right)\left(\left\lfloor{Qm}\right\rfloor+1\right)\\
=\frac M2\sum_{m=0}^{M-1}\left(Qm+\{Qm\}\right)\left(Qm-\{Qm\}+1\right)\\
=\frac M2\sum_{m=0}^{M-1}\left((Qm)^2+Qm-\{Qm\}^2+\{Qm\}\right).$$ 
As in the sums of fractional parts, all fractions from $0/M$ to $(M-1)/M$ appear (in disorder),
$$12S=6M\left(Q^2\frac{(M-1)M(2M-1)}6+Q\frac{(M-1)M}2-\frac{(M-1)M(2M-1)}{6M^2}+\frac{(M-1)M}{2M}\right)\\
=2M^2N^2+M^2-3MN^2-3MN+N^2+3NM^2-1.$$
To this, we need to add the omitted term ($m=M$), which is found to be
$$T:=\sum_{n=0}^{N}(NM-Mn)=\frac{MN(N+1)}2,$$
and finally
$$2(S+T)=\frac{2M^2N^2+M^2+9MN^2+9MN+N^2+3NM^2-1}6.$$
Now if $M,N$ are not relative primes, one may split the summation in $G=\gcd(M,N)$ subsums of $\dfrac MG$ terms and a final $M^{th}$ term. This amounts to replacing the final $-1$ by $-G^2$.
A: Partial solution:
We can break the expression into four parts $E_1 + E_2 + E_3 + E_4$
where $$E_1 = \sum_{\substack{0\le n\le N/2 \\ 0\le m\le M/2}}(MN -mN-nM) = $$
$$E_2 = \sum_{\substack{0\le n\le N/2 \\ M/2\le m\le M}}|MN -mN-nM|$$
$$E_3 = \sum_{\substack{N/2\le n\le N \\ 0\le m\le M/2}}|MN -mN-nM|$$
$$E_4 = \sum_{\substack{N/2\le n\le N \\ M/2\le m\le M}}(mN + nM - MN) = \sum_{\substack{0\le n\le N/2 \\ 0\le m\le M/2}}(mN + nM)$$
$E_1 + E_4 = M^2N^2/4$
Both $E_2$ and $E_3$ can be rewritten as 
$$ E_2 = E_3 = \sum_{\substack{0\le n\le N/2 \\ 0\le m\le M/2}}|MN/2 -mN-nM|$$
so $$E_2 + E_3 = \sum_{\substack{0\le n\le N/2 \\ 0\le m\le M/2}}|MN -2mN-2nM|$$
So $$LHS = M^2N^2/4 + \sum_{\substack{0\le n\le N/2 \\ 0\le m\le M/2}}|MN -2mN-2nM|$$ 
and may be we can recursively move forward?
