bzoj 4659 Lcm 莫比乌斯反演

题面

题目传送门

解法

还是写得详细一点比较好

• 我们可以比较显然地得出式子：${\sum }_{i=1}^n{\sum }_{j=1}^n\mu (gcd(i,j){)}^{2}lcm(i,j)$
• 然后下面就是推导过程了：
  $$=\sum _{d=1}^n\mu (d{)}^{2}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)==1]ij$$
• 通过反演，我们可以得到下面的式子：
  $$=\sum _{d=1}^n\mu (d{)}^{2}d\sum _{d_1}^{\lfloor \frac{n}{d}\rfloor }\mu (d_1)d_1^2\sum _{i=1}^{\lfloor \frac{n}{d{d}_1}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d{d}_1}\rfloor }ij$$
• 那么我们令$s(n)={\sum }_{i=1}^{n}i$，就可以得到下面的式子：
  $$=\sum _{d=1}^n\mu (d{)}^{2}d\sum _{d_1}^{\lfloor \frac{n}{d{d}_1}\rfloor }\mu (d_1)d_1^2\ast s(\lfloor \frac{n}{d{d}_1}\rfloor {)}^2$$
• 然后枚举$d{d}_1=k$，可以得到：
  $$=\sum _{k=1}^{n}s(\lfloor \frac{n}{k}\rfloor {)}^2\sum _{d|k}\mu (d{)}^{2}d\ast \mu (\frac{k}{d})(\frac{k}{d}{)}^2$$
• 我们令$f(k)={\sum }_{d|k}\mu (d{)}^{2}d\ast \mu (\frac{k}{d})(\frac{k}{d}{)}^2$，因为$n\le 4\ast 10^6$，所以使用$O(n\log n)$的筛法是不能通过的。
• 那么我们考虑如何线性筛这个函数$f(k)$。设$k=p_1^{k_1}\dots p_m^{k_m}$，显然，如果$\exists k_i>2$，那么$f(k)=0$。
• 可以发现，$f$是一个积性函数。那么我们只要处理出$f(p^k)$即可
• 因为$k\le 2$，所以我们对于$k=1/2$分别讨论
• 显然可以发现，$f(p)=-p^2+p,f(p^2)=-p^3$，所以我们是可以线性筛的。
• 时间复杂度：$O(n+T\sqrt{n})$
• 因为是对$2^{30}$取模，所以我们可以直接自然溢出，然后$\&(2^{30}-1)$即可。

代码

#include <bits/stdc++.h>
#define Mod ((1 << 30) - 1)
#define N 4000010
using namespace std;
template <typename T> void chkmax(T &x, T y) {x = x > y ? x : y;}
template <typename T> void chkmin(T &x, T y) {x = x > y ? y : x;}
template <typename T> void read(T &x) {
	x = 0; int f = 1; char c = getchar();
	while (!isdigit(c)) {if (c == '-') f = -1; c = getchar();}
	while (isdigit(c)) x = x * 10 + c - '0', c = getchar(); x *= f;
}
bool f[N]; int g[N], p[N], mx[N], mns[N], tmp[N];
int s(int n) {return n * (n + 1) / 2;}
void add(int &x, int y) {x += y;}
void sieve(int n) {
	memset(f, true, sizeof(f)); int len = 0; g[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (f[i]) {
			p[++len] = i, g[i] = i - 1ll * i * i;
			mns[i] = mx[i] = tmp[i] = 1;
		}
		for (int j = 1; j <= len && i * p[j] <= n; j++) {
			int k = i * p[j]; f[k] = false;
			if (i % p[j] == 0) {
				tmp[k] = tmp[i], mns[k] = mns[i] + 1;
				mx[k] = max(mx[tmp[k]], mns[k]);
				if (mx[k] > 2) {g[k] = 0; continue;}
				g[k] = g[tmp[k]] * (-p[j] * p[j] * p[j]);
				break;
			} else {
				tmp[k] = i, mns[k] = 1, mx[k] = max(mx[i], 1);
				if (mx[k] > 2) {g[k] = 0; continue;}
				g[k] = g[i] * (p[j] - p[j] * p[j]);
			}
		}
	}
	for (int i = 1; i <= n; i++) g[i] += g[i - 1];
}
int solve(int n, int m) {
	int ret = 0, x = 0;
	for (int i = 1; i <= n; i = x + 1) {
		x = min(n / (n / i), m / (m / i));
		add(ret, (g[x] - g[i - 1]) * s(n / i) * s(m / i));
	}
	return ret & Mod;
}
int main() {
	sieve(4e6); int T; read(T);
	while (T--) {
		int n, m; read(n), read(m);
		if (n > m) swap(n, m);
		cout << solve(n, m) << "\n";
	}
	return 0;
}