[国家集训队]Crash的数字表格 / JZPTAB

link
又是喜闻乐见的推式子环节。
以下令 \(n\le m\)。
\(ans=\sum_{i=1}^n \sum_{j=1}^m \mathrm {lcm}(i,j)\)
\(=\sum_{i=1}^n \sum_{j=1}^m \frac {i\times j}{\gcd(i,j)}\)
\(=\sum_{d=1}^n d\times \sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[\gcd(i,j)=1]\times i\times j\)
\(=\sum_{d=1}^n d\times\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{k|\gcd(i,j)}\mu_k\times i\times j\)

Way 1：把 \(\mu\) 推到最后
令 \(1\sim x = 1+2+3+..+x\)。
\(=\sum_{d=1}^n d\times \sum_{g=1}^{n} g^2 \times (1\sim \lfloor n/d/g \rfloor)(1\sim \lfloor m/d/g \rfloor) * \sum_{k|g}\mu_k\)
令 \(T=dg\)。
\(=\sum_{T=1}^{n} \lfloor n/T\rfloor \lfloor m/T \rfloor \times T\times \sum_{d|T} d\times \mu_{\frac T d}\times (\frac T d)^2\)
后面那两个 \(\sum\) 可以用类似埃筛的方法预处理，时间复杂度：\(\mathcal {O(nlog_2(n))}\)。虽然这个A不了本题，但却是很常见的一种方法，可以做到单次 \(\mathcal {O(\sqrt n)}\)。

update 2022.2.11：由狄利克雷卷积知 \(\sum_{d|T} d\times \mu_{\frac T d}\times (\frac T d)^2\) 是一个积性函数，于是直接线性筛出来即可。\(\mathcal {O}(n)+\mathcal {O}(\sqrt n+\sqrt m)\)。

Way2：把 \(\mu\) 推到前面。
\(=\sum_{d=1}^n d\times\sum_{k=1}^{n} \mu_k \times k^2 \times(1\sim n/d/k) \times (1\sim m/d/k)\)
细心的读者可以发现，这里可以数论分块套数论分块做到单次 \(\mathcal {O}(n+m)\)。但预处理是 \(\mathcal {O}(n)\)的。
代码咕咕咕咕咕~做文化课去了~~~~~
Code Way2（注意取模相减后是负数的情况）

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#define uint unsigned int
#define LL long long
using namespace std;
const int Mod = 20101009, MAXN = 1e7 + 5;
int prime[MAXN / 10], tot, n, m;
LL mu[MAXN], pre[MAXN];
bool vis[MAXN];
void Prime() {
	mu[1] = 1;
	for(int i = 2; i <= 10000000; i ++) {
		if(!vis[i]) prime[++ tot] = i, mu[i] = -1;
		for(int j = 1; j <= tot && i * prime[j] <= 10000000; j ++) {
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) { mu[i * prime[j]] = 0; break; }
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= 10000000; i ++) pre[i] = pre[i - 1] + mu[i] * i * i, pre[i] %= Mod;
} 
LL Getsum(int l, int r) {
	return (LL)(l + r) * (r - l + 1) / 2 % Mod;
}
LL Calc1(int u, int v) {
	LL ans = 0;
	for(int l = 1, r; l <= u && l <= v; l = r + 1) {
		r = min(u / (u / l), v / (v / l));
		ans = (ans + (pre[r] - pre[l - 1]) * Getsum(1, u / l) % Mod * Getsum(1, v / l) % Mod) % Mod;
	}
	return ans;
}
LL Calc() {
	LL ans = 0;
	for(int l = 1, r; l <= n && l <= m; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += Getsum(l, r) * Calc1(n / l, m / l); ans %= Mod;
	}
	ans = (ans + Mod) % Mod;
	return ans;
}
int main() {
	Prime(); scanf("%d%d", &n, &m);
	printf("%lld", Calc());
	return 0;
}