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.
$\def\peq{\mathrel{\phantom{=}}{}}$Suppose $M, N &gt; 0$ and $d = (M, N)$, $M = md$, $N = nd$.\begin{align*}
\sum_{\substack{0\leqslant k\leqslant M\\0\leqslant j\leqslant N}} |kN - jM| &amp;= \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1}} |kN - jM| + \sum_{k = 0}^{M - 1} (NM - kN) + \sum_{j = 0}^{N - 1} (MN - jM)\\
&amp;= \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1}} |kN - jM| + \frac{1}{2} NM(M + 1) + \frac{1}{2} MN(N + 1). \tag{1}
\end{align*}
Note that$$
\sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1}} (kN - jM) = \frac{1}{2} NM(M - 1) - \frac{1}{2} MN(N - 1),
$$
then\begin{align*}
\sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1}} |kN - jM| &amp;= \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN\geqslant jM}} (kN - jM) - \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN&lt;jM}} (kN - jM)\\
&amp;= 2\sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN\geqslant jM}} (kN - jM) - \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1}} (kN - jM)
\end{align*}
which implies$$
(1) = 2\sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN\geqslant jM}} (kN - jM) + MN(N + 1). \tag{2}
$$

Now for each $0 \leqslant k \leqslant M - 1$, define $$r_k = kn - m\left[ \dfrac{kn}{m} \right],$$ i.e. $r_k$ is the remainder of $kn$ modulo $m$.
\begin{align*}
&amp;\peq \sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN\geqslant jM}} (kN - jM) = \sum_{\substack{0\leqslant k\leqslant dm-1\\0\leqslant j\leqslant dn-1\\kdn\geqslant jdm}} (kdn - jdm) = d \sum_{\substack{0\leqslant k\leqslant dm-1\\0\leqslant j\leqslant dn-1\\kn\geqslant jm}} (kn - jm)\\
&amp;= d\sum_{k = 0}^{dm - 1} \sum_{j = 0}^{\left[ \frac{kn}{m} \right]} (kn - jm) = d\sum_{k = 0}^{dm - 1} \left( kn \left( \left[ \frac{kn}{m} \right] + 1 \right) - m·\frac{1}{2} \left[ \frac{kn}{m} \right] \left( \left[ \frac{kn}{m} \right] + 1 \right) \right)\\
&amp;= d\sum_{k = 0}^{dm - 1} \left( kn - \frac{m}{2} \left[ \frac{kn}{m} \right] \right) \left( \left[ \frac{kn}{m} \right] + 1 \right) = d\sum_{k = 0}^{dm - 1} \left( kn - \frac{m}{2}·\frac{kn - r_k}{m} \right) \left( \frac{kn - r_k}{m} + 1 \right)\\
&amp;= d\sum_{k = 0}^{dm - 1} \left( \frac{1}{2} kn + \frac{1}{2} r_k \right) \left( \frac{kn}{m} + 1 - \frac{r_k}{m} \right) = \frac{d}{2m} \sum_{k = 0}^{dm - 1} (kn(kn + m) + m r_k - r_k^2). \tag{3}
\end{align*}

Note that $(m, n) = 1$ implies for each $0 \leqslant s \leqslant d - 1$,$$\{kn \mid sm \leqslant k \leqslant (s + 1)m - 1\}$$ is a complete residue system modulo $m$, therefore$$
\{r_{sm}, \cdots, r_{(s + 1)m - 1}\} = \{ 0, \cdots, m - 1\}.
$$
Thus$$
\sum_{k = 0}^{dm - 1} r_k = \sum_{s = 0}^{d - 1} \sum_{k = sm}^{(s + 1)m - 1} r_k = \sum_{s = 0}^{d - 1} \sum_{k = 0}^{m - 1} k = \frac{1}{2} m(m - 1)·d,\\
\sum_{k = 0}^{dm - 1} r_k^2 = \sum_{s = 0}^{d - 1} \sum_{k = 0}^{m - 1} r_k^2 = \sum_{s = 0}^{d - 1} \sum_{k = 0}^{m - 1} k^2 = \frac{1}{6} (2m - 1)(m - 1)m·d,
$$
then\begin{align*}
(3) &amp;= \frac{d}{2m} \sum_{k = 0}^{dm - 1} kn(kn + m) + \frac{d}{2} \sum_{k = 0}^{dm - 1} r_k - \frac{d}{2m} \sum_{k = 0}^{dm - 1} r_k^2\\
&amp;= \frac{dn^2}{2m} \sum_{k = 0}^{dm - 1} k^2 + \frac{1}{2} dn \sum_{k = 0}^{dm - 1} k + \frac{d}{2} × \frac{1}{2} m(m - 1)·d - \frac{d}{2m} × \frac{1}{6} (2m - 1)(m - 1)m·d\\
&amp;= \frac{dn^2}{2m} × \frac{1}{6} (2dm - 1)(dm - 1)dm + \frac{1}{2} dn × \frac{1}{2} (dm - 1)dm\\
&amp;\peq + \frac{d}{2} × \frac{1}{2} m(m - 1)·d - \frac{d}{2m} × \frac{1}{6} (2m - 1)(m - 1)m·d\\
&amp;= \frac{1}{6} d^4 m^2 n^2 + \frac{1}{4} d^3 mn(m - n) + \frac{1}{12} d^2 (m^2 + n^2 - 3mn) - \frac{1}{12} d^2.
\end{align*}

Combining with (2),\begin{align*}
\sum_{\substack{0\leqslant k\leqslant M\\0\leqslant j\leqslant N}} |kN - jM| &amp;= 2\sum_{\substack{0\leqslant k\leqslant M-1\\0\leqslant j\leqslant N-1\\kN\geqslant jM}} (kN - jM) + MN(N + 1)\\
&amp;= 2\left( \frac{1}{6} d^4 m^2 n^2 + \frac{1}{4} d^3 mn(m - n) + \frac{1}{12} d^2 (m^2 + n^2 - 3mn) - \frac{1}{12} d^2 \right)\\
&amp;\peq + d^2 mn(dn + 1)\\
&amp;= \frac{1}{3} d^4 m^2 n^2 + \frac{1}{2} d^3 mn(m + n) + \frac{1}{6} d^2 (m^2 + n^2 + 3mn) - \frac{1}{6} d^2\\
&amp;=\frac{1}{6} \left( MN(2MN + 3(M + N + 1)) + M^2 + N^2-(M,N)^2\right).
\end{align*}
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?