Calculating $\sum\limits_{i=1}^{n-1}\sum\limits_{j=1}^{m-1} |mi-nj|$ I'd like to calculate $f(n,m)=\sum\limits_{i=1}^{n-1}\sum\limits_{j=1}^{m-1} |mi-nj|$ for all $1 \leq n \leq N,\ 1 \leq m \leq M$. Straightforward brute force method runs in $O(N^2M^2)$ which is too slow.  How to calculate all values in $O(NM)$?
 A: 
We show the following is valid for positive integers $n,m$:
  \begin{align*}
\sum_{i=1}^{n-1}\sum_{j=1}^{m-1}|mi-nj|=\frac{1}{6}\left(2m^2n^2-3m^2n+m^2-3mn^2+3mn+n^2-\left(\gcd(m,n)\right)^2\right)
\end{align*}

In the following we denote with   $d=\gcd(m,n)$.

We obtain
  \begin{align*}
\color{blue}{\sum_{i=1}^{n-1}}&\color{blue}{\sum_{j=1}^{m-1}|mi-nj|}\tag{2}\\
&=2\sum_{i=1}^{n-1}\sum_{j=1}^{\lfloor mi/n\rfloor}(mi-nj)\tag{3}\\
&=2m\sum_{i=1}^{n-1}i\sum_{j=1}^{\lfloor mi/n\rfloor}1-2n\sum_{i=1}^{n-1}\sum_{j=1}^{\lfloor mi/n\rfloor}j\\
&=2m\sum_{i=1}^{n-1}i\left\lfloor\frac{m}{n}i\right\rfloor
-n\sum_{i=1}^{n-1}\left\lfloor\frac{m}{n}i\right\rfloor\left(\left\lfloor\frac{m}{n}i\right\rfloor+1\right)\tag{4}\\
&=\sum_{i=1}^{n-1}\left(2mi-n\left\lfloor\frac{m}{n}i\right\rfloor-n\right)\left\lfloor\frac{m}{n}i\right\rfloor\tag{5}\\
&=\sum_{i=1}^{n-1}\left(2mi-n\left(\frac{m}{n}i-\left\{\frac{m}{n}i\right\}\right)-n\right)\left(\frac{m}{n}i-\left\{\frac{m}{n}i\right\}\right)\tag{6}\\
&=n\sum_{i=1}^{n-1}\left(\frac{m^2}{n^2}i^2-\left\{\frac{m}{n}i\right\}^2-\frac{m}{n}i+\left\{\frac{m}{n}i\right\}\right)\\
&=\frac{m^2}{n}\sum_{i=1}^{n-1}i^2-n\sum_{i=1}^{n-1}\left\{\frac{m}{n}i\right\}^2-m\sum_{i=1}^{n-1}i+n\sum_{i=1}^{n-1}\left\{\frac{m}{n}i\right\}\\
&=\frac{m^2}{n}\frac{1}{6}(n-1)n(2n-1)-nd\sum_{i=0}^{n/d-1}\left(\frac{d}{n}i\right)^2\\
&\qquad -m\frac{1}{2}(n-1)n+nd\sum_{i=0}^{n/d-1}\left(\frac{d}{n}i\right)\tag{7}\\
&=\frac{1}{6}m^2(n-1)(2n-1)-\frac{d^3}{n}\frac{1}{6}\left(\frac{n}{d}-1\right)\frac{n}{d}\left(\frac{2n}{d}-1\right)\\
&\qquad-\frac{1}{2}mn(n-1)+d^2\frac{1}{2}\left(\frac{n}{d}-1\right)\frac{n}{d}\tag{8}\\
&\,\,\color{blue}{=\frac{1}{6}\left(2m^2n^2-3m^2n+m^2-3mn^2+3mn+n^2-d^2\right)}
\end{align*}
  and the claim (1) follows.

Comment:


*

*In (3) we use that positive and negative parts in (2) correspond to each other.

*In (4) we expand the inner sums.

*In (5) we rearrange the terms and factor out $\left\lfloor\frac{m}{n}i\right\rfloor$.

*In (6) we rewrite the expression using the fractional part $\{x\}=x-\lfloor x\rfloor$ of $x$.

*In (7) we expand the sums with linear and quadratic terms and we apply the identity 
\begin{align*}
\sum_{i=1}^{n}f\left(\left\{\frac{m}{n}i\right\}\right)=d\sum_{i=0}^{n/d-1}f\left(\frac{d}{n}i\right)
\end{align*}
where $d=\gcd(m,n)$.

*In (8) we expand the sums and simplify the expression in the final step.
A: Given $n$ and $m\geq n$ you can determine the numbers
$$r_{nm}(j):=\left\lfloor{m j\over n}\right\rfloor\qquad(1\leq j\leq n-1)$$in $O(n)$ steps, and another $O(n)$ steps then compute
$$f(n,m)=\sum_{j=1}^{n-1}r_{nm}(j)\bigl(2 m j-n(r_{nm}(j)+1)\bigr)\ .\tag{1}$$
One arrives at this formula after observing that by symmetry it is sufficient to consider the lattice points $(j,k)$ below the line $y={m\over n}x$. Look at these lattice points on the ordinate $x=j$. Their $y$-coordinates $k$  run through the set $[r_{nm}(j)]$. The"integrand" $p(j,k):=|m j-n k|=mj -n k$ produces an arithmetic progression on these lattice points. Therefore the sum of the $p$-values along the ordinate $x=j$ is their number $r_{nm}(j)$, times the arithmetic mean of the first and the last $p$-values. This leads to formula $(1)$.
