Description：

求： $$ \sum_{i=1}^n\sum_{j=1}^m(n\% i)(m\% j)(i\ne j) $$ $1\leqslant n,m\leqslant 10^9$

Solution：

$$ \begin{align} 原式&=\sum_{i=1}^n\sum_{j=1}^m(n\%i)(m\%j)-\sum_{i=1}^{\min(n,m)}(n\%i)(m\%i)\\ \end{align} $$

分两部分考虑，前半部分： $$ \begin{align} \sum_{i=1}^n\sum_{j=1}^m(n\%i)(m\%j)&=\sum_{i=1}^n\sum_{j=1}^mnm-\sum_{i=1}^n\sum_{j=1}^mm\lfloor\frac n i\rfloor i-\sum_{i=1}^n\sum_{j=1}^mn\lfloor\frac m j\rfloor j+\sum_{i=1}^n\sum_{j=1}^m\lfloor\frac n i\rfloor i\lfloor\frac m j\rfloor j\\ &=n^2m^2-m^2\sum_{i=1}^n\lfloor\frac n i\rfloor i-n^2\sum_{j=1}^m\lfloor\frac m j\rfloor j+(\sum_{i=1}^n\lfloor\frac n i\rfloor i)\times(\sum_{j=1}^m\lfloor\frac m j\rfloor j) \end{align} $$ 这个就是一堆除法分块。

后半部分： $$ \begin{align} \sum_{i=1}^{\min(n,m)}(n\%i)(m\%i)&=\sum_{i=1}^{\min(n,m)}(nm-m\lfloor\frac n i\rfloor i-n\lfloor\frac m i\rfloor i+\lfloor\frac n i\rfloor\lfloor\frac m i\rfloor i^2)\\ &=\min(n,m)nm+\sum_{i=1}^{\min(n,m)}(\lfloor\frac n i\rfloor\lfloor\frac m i\rfloor i^2-m\lfloor\frac n i\rfloor i-n\lfloor\frac m i\rfloor i) \end{align} $$ 后面那个式子之所以不拆开是因为拆开了 $\min(n,m)$ 就不好处理了。

都可以除法分块做，总复杂度 $O(\sqrt n)$ 。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,m;

#define MOD 19940417

int sum1(int k){return 1ll * k * (k + 1) / 2 % MOD;}

int calc1(int n)

{

	int ans = 0;

	for(int l = 1,r;l <= n;l = r + 1)

	{

		r = n / (n / l);

		ans = (ans + 1ll * (n / l) * (sum1(r) - sum1(l - 1) + MOD) % MOD) % MOD;

	}

	return ans;

}

#define INV6 3323403

int sum2(int k){return 1ll * k * (k + 1) % MOD * (2 * k + 1) % MOD * INV6 % MOD;}

int calc2(int n,int m)

{

	int ans1 = 0,ans2 = 0;

	for(int l = 1,r;l <= min(n,m);l = r + 1)

	{

		r = min(n / (n / l),m / (m / l));

		ans1 = (ans1 + 1ll * (n / l) * (m / l) % MOD * (sum2(r) - sum2(l - 1) + MOD) % MOD) % MOD;

		int s1 = 1ll * m * (n / l) % MOD * (sum1(r) - sum1(l - 1) + MOD) % MOD;

		int s2 = 1ll * n * (m / l) % MOD * (sum1(r) - sum1(l - 1) + MOD) % MOD;

		ans2 = (ans2 - (s1 + s2) % MOD + MOD) % MOD;

	}

	return (ans1 + ans2) % MOD;

}

int main()

{

	cin >> n >> m;

	int ans1 = 0,ans2 = 0;

	ans1 = (ans1 + 1ll * n * n % MOD * m % MOD * m % MOD) % MOD;

	ans1 = (ans1 - 1ll * m * m % MOD * calc1(n) % MOD + MOD) % MOD;

	ans1 = (ans1 - 1ll * n * n % MOD * calc1(m) % MOD + MOD) % MOD;

	ans1 = (ans1 + 1ll * calc1(n) * calc1(m) % MOD) % MOD;

	ans2 = (ans2 + 1ll * min(n,m) * n % MOD * m % MOD);

	ans2 = (ans2 + calc2(n,m)) % MOD;

	cout << (ans1 - ans2 + MOD) % MOD << endl;

	return 0;

}