P1829 Crash 的数字表格（莫比乌斯反演）

Description

求
$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)$$

Solution

原式等价于
$$\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{gcd(i,j)}$$
老套路，令 $gcd(i,j)=d$，枚举一下 $d$
$$\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]\times \frac{ij}{d}$$
提出 $d$，$d$ 对 $\frac{ij}{d}$ 有贡献
$$\sum _{d=1}^{n}d\times \sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{m}{d}}[gcd(i,j)=1]\times ij$$
把 $[gcd(i,j)=1]$ 换掉
$$\sum _{d=1}^{n}d\times \sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{m}{d}}\sum _{k|gcd(i,j)}\mu (k)\times ij$$
套路地枚举 $k$
$$\sum _{d=1}^{n}d\times \sum _{k=1}^{\frac{n}{d}}\mu (k)\times \sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{m}{d}}[k|gcd(i,j)]\times ij$$
继续
$$\sum _{d=1}^{n}d\times \sum _{k=1}^{\frac{n}{d}}\mu (k)\times \sum _{ik=1}^{\frac{n}{d}}\sum _{jk=1}^{\frac{m}{d}}\times ij\times k^2$$

$$=\sum _{d=1}^{n}d\times \sum _{k=1}^{\frac{n}{d}}\mu (k)\times \sum _{i=1}^{\frac{n}{dk}}\sum _{j=1}^{\frac{m}{dk}}\times ij\times k^2$$

归位
$$=\sum _{d=1}^{n}d\times \sum _{k=1}^{\frac{n}{d}}\mu (k)k^2\times \sum _{i=1}^{\frac{n}{dk}}\sum _{j=1}^{\frac{m}{dk}}\times ij$$
用等差数列替换
$$=\sum _{d=1}^{n}d\times \sum _{k=1}^{\frac{n}{d}}\mu (k)k^2\times \frac{\frac{n}{dk}(\frac{n}{dk}+1)}{2}\times \frac{\frac{m}{dk}(\frac{m}{dk}+1)}{2}$$
筛一下 $\mu (k)k^2$，然后数论分块套数论分块。

复杂度约为 $O(n)$，要到处膜。

Code

#include <bits/stdc++.h>
using namespace std;
#define re register
#define F first
#define S second
typedef long long ll;
typedef pair<int, int> P;
const int N = 1e7 + 5, M = 1e7 + 5, mod = 20101009;
const int INF = 0x3f3f3f3f;
inline int read() {
    int X = 0,w = 0; char ch = 0;
    while(!isdigit(ch)) {w |= ch == '-';ch = getchar();}
    while(isdigit(ch)) X = (X << 3) + (X << 1) + (ch ^ 48),ch = getchar();
    return w ? -X : X;
}
int p[N], mu[N];
ll sum[N];
bool vis[N];
int n, m;
void init(){
	int cnt = 0, k = min(n, m); mu[1] = 1;
	for (int i = 2; i <= k; i++){
		if (!vis[i]) p[++cnt] = i, mu[i] = -1;
		for (int j = 1; j <= cnt && i * p[j] <= k; j++){
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) {
				mu[i * p[j]] = 0; break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
	for (int i = 1; i <= N; i++) sum[i] = (sum[i - 1] + 1ll * i * i % mod * (mu[i] + mod) % mod) % mod;
}
ll get(int n, int m){
	return (1ll * n * (n + 1) / 2 % mod)  * (1ll * m * (m + 1) / 2 % mod) % mod;
}
ll solsum(int n, int m){
	ll ans = 0;
	for (int i = 1, j; i <= min(n, m); i = j + 1){
		j = min(n / (n / i), m / (m / i));
		ans = (ans + 1ll * (sum[j] - sum[i - 1] + mod) % mod * get(n / i, m / i) % mod) % mod;
	}
	return ans;
}
ll solve(int n, int m){
	ll ans = 0;
	for (int i = 1, j; i <= min(n, m); i = j + 1){
		j = min(n / (n / i), m / (m / i));
		ans = (ans + 1ll * (j - i + 1) * (i + j) / 2 % mod * solsum(n / i, m / i) % mod) % mod;
	}
	return ans;
}
int main(){
	n = read(), m = read();
	init();
	printf("%lld\n", solve(n, m));
	return 0;
}