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.
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} $$