洛谷 P1829 ：[国家集训队]Crash的数字表格 / JZPTAB（莫比乌斯反演）

---

题解：
假设$n<m$，一波化简：$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)=\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast j}{gcd(i,j)}=\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast j}{d}\ast [gcd(i,j)=d]=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d\ast [gcd(i,j)=1]$$
$$设：f(d)=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d\ast [gcd(i,j)=d]$$
$$令F(x)=\sum _{x|d}f(d)=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d\ast [x|gcd(i,j)]=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{x\ast d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{x\ast d}\rfloor }i\ast j\ast d\ast x^2$$
$$反演一下：f(d)=\sum _{d|x}\mu (\frac{x}{d})\ast F(x)$$
$$f(1)=\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\mu (\frac{x}{d})\ast F(x)=\sum _{d=1}^{n}d\ast \sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }{x}^2\ast \mu (x)\ast \sum _{i=1}^{\lfloor \frac{n}{d\ast x}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{d\ast x}\rfloor }j$$
到这一步已经可以两层分块$O(n)$过掉这题 (准确来说是 $O(n^{\frac{3}{4}}$)?)
这里为什么 $x$ 的上界是 $\lfloor \frac{n}{d}\rfloor$：因为前面设了$$f(d)=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d\ast [gcd(i,j)=d]$$
$x$是d的倍数，由于 i 上限是 $\lfloor \frac{n}{d}\rfloor$（这里假设 $n<m$），因此d的倍数最大到$\lfloor \frac{n}{d}\rfloor$（准确的说是gcd(i,j)最大到$\lfloor \frac{n}{d}\rfloor$），$x$ 应该枚举到$\lfloor \frac{n}{d}\rfloor$
注意到最终式子里枚举了 $d,x$，产生了 $d\ast x$，可以令 T = $d\ast x$，枚举 T 和 d，然后可以线性筛优化一下，解决多组的情况。

注意第一层分块的时候右边界取 $min(m/(m/i),n/(n/i))$，因为后面有一个 $\lfloor \frac{\lfloor \frac{m}{d}\rfloor }{x}\rfloor$

---

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod = 20101009;
const int maxn = 1e7 + 10;
bool ispri[maxn];
int pri[maxn],mu[maxn];
typedef long long ll;
ll sum[maxn];
int n,m;
void sieve(int n) {
	ispri[0] = ispri[1] = true;
	pri[0] = 0;mu[1] = 1;
	for(int i = 2; i <= n; i++) {
		if(!ispri[i]) pri[++pri[0]] = i,mu[i] = -1;
		for(int j = 1; j <= pri[0] && i * pri[j] <= n; j++) {
			ispri[i * pri[j]] = true;
			if(i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
	sum[0] = 0;
	for(int i = 1; i <= n; i++)
		sum[i] = (sum[i - 1] + mu[i] * 1ll * i % mod * i % mod + mod) % mod;
}
ll cal(int n) {
	return 1ll * (1 + n) * n / 2 % mod;
}
int main() {
	sieve(10000000);
	scanf("%d%d",&n,&m);
	int l1,r1,l2,r2;
	ll res = 0;
	if(n > m) swap(n,m);
	for(l1 = 1; l1 <= n; l1 = r1 + 1) {
		r1 = min(n / (n / l1),m / (m / l1));
		ll tmp = 0;
		for(l2 = 1; l2 <= n / l1; l2 = r2 + 1) {
			int p1 = n / l1,p2 = m / l1;
			r2 = min(p1 / (p1 / l2),p2 / (p2 / l2));
			tmp = (tmp + cal(p1 / l2) * cal(p2 / l2) % mod * ((sum[r2] - sum[l2 - 1] + mod) % mod) % mod) % mod;
		}
		res = (res + 1ll * (l1 + r1) * (r1 - l1 + 1) / 2 % mod * tmp % mod + mod) % mod;
	}	
	printf("%lld\n",res);
	return 0;
}