洛谷P4240 毒瘤之神的考验 【莫比乌斯反演 + 分块打表】

题目链接

洛谷P4240

题解

式子不难推，分块打表真的没想到
首先考虑如何拆开\(\varphi(ij)\)
考虑公式

\[\varphi(ij) = ij\prod\limits_{p | ij}\frac{p - 1}{p} \]

而

\[\begin{aligned} \varphi(i)\varphi(j) &= i\prod\limits_{p | i}\frac{p - 1}{p} j \prod\limits_{p | j}\frac{p - 1}{p} \\ \varphi(i)\varphi(j)&= ij \prod\limits_{p | ij}\frac{p - 1}{p} \prod\limits_{p | (i,j)}\frac{p - 1}{p} \\ \varphi(i)\varphi(j)(i,j)&= ij \prod\limits_{p | ij}\frac{p - 1}{p} (i,j)\prod\limits_{p | (i,j)}\frac{p - 1}{p} \\ \varphi(i)\varphi(j)(i,j)&= \varphi(i,j)\varphi((i,j)) \\ \varphi(i,j) &= \frac{\varphi(i)\varphi(j)(i,j)}{\varphi((i,j))} \end{aligned} \]

所以我们有

\[\begin{aligned} ans &= \sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m} \frac{\varphi(i)\varphi(j)(i,j)}{\varphi((i,j))} \\ &= \sum\limits_{d = 1}^{n}\frac{d}{\varphi(d)} \sum\limits_{d|i}^{n} \varphi(i) \sum\limits_{d|j}^{m} \varphi(j) [(i,j) = d] \\ &= \sum\limits_{z = 1}^{n}\frac{z}{\varphi(z)} \sum\limits_{z | d} \mu(\frac{d}{z}) \sum\limits_{d|i}^{n} \varphi(i) \sum\limits_{d|j}^{m} \varphi(j) \\ &= \sum\limits_{z = 1}^{n}\frac{z}{\varphi(z)} \sum\limits_{z | d} \mu(\frac{d}{z}) \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \varphi(id) \sum\limits_{j = 1}^{\lfloor \frac{m}{d} \rfloor} \varphi(jd) \\ &= \sum\limits_{d = 1}^{n}\sum\limits_{z | d} \frac{z}{\varphi(z)}\mu(\frac{d}{z}) \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \varphi(id) \sum\limits_{j = 1}^{\lfloor \frac{m}{d} \rfloor} \varphi(jd) \end{aligned} \]

隐约感觉已经差不多了
令

\[F(d) = \sum\limits_{z | d} \frac{z}{\varphi(z)}\mu(\frac{d}{z}) \]

\[G(n,d) = \sum\limits_{k = 1}^{n}\varphi(kd) \]

那么

\[ans = \sum\limits_{d = 1}^{n}F(d)G(\lfloor \frac{n}{d} \rfloor,d)G(\lfloor \frac{m}{d} \rfloor,d) \]

注意到\(\lfloor \frac{n}{d} \rfloor d \le n\)，\(G(n,d)\)的状态数是\(O(n)\)的，可以预处理
\(F(d)\)可以\(O(nlogn)\)预处理
那么现在我们就可以\(O(n)\)计算辣~
但这还不够，却已经无法从式子上进行优化了
分块！

我们设

\[T(n,x,y) = \sum\limits_{d = 1}^{n}F(d)G(x,d)G(y,d) \]

我们就可以利用整除分块\(T(nxt,\lfloor \frac{n}{i} \rfloor,\lfloor \frac{m}{i} \rfloor) - T(i - 1,\lfloor \frac{n}{i} \rfloor,\lfloor \frac{m}{i} \rfloor)\)进行\(O(\sqrt{n})\)的计算
所以我们只需预处理出\(T(n,x,y)\)
但状态数很多，考虑只打出\(B\)项\(x,y\)
考虑到\(nx,ny \le 10^5\)，时间空间复杂度\(O(nB)\)
\(y \ge B\)的\(n\)只有\(\lfloor \frac{n}{B} \rfloor\)项，暴力计算
所以询问时复杂度\(O(\sqrt{n} + \lfloor \frac{n}{B} \rfloor)\)
令\(T\frac{n}{B} = nB\)，则\(B = \sqrt{T}\)
\(B\)取\(100\)左右就差不多了

复杂度\(O(nlogn + T(\sqrt{n} + \lfloor \frac{n}{B} \rfloor))\)

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<vector>
#include<queue>
#include<cmath>
#include<map>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define cls(s,v) memset(s,v,sizeof(s))
#define mp(a,b) make_pair<int,int>(a,b)
#define cp pair<int,int>
using namespace std;
const int maxn = 100005,maxm = 100005,INF = 0x3f3f3f3f;
inline int read(){
	int out = 0,flag = 1; char c = getchar();
	while (c < 48 || c > 57){if (c == '-') flag = 0; c = getchar();}
	while (c >= 48 && c <= 57){out = (out << 1) + (out << 3) + c - 48; c = getchar();}
	return flag ? out : -out;
}
const int B = 100,N = 100000,P = 998244353;
int p[maxn],pi,isn[maxn],mu[maxn],phi[maxn],inv[maxn];
int F[maxn],*G[maxn],*T[B + 1][B + 1];
void init(){
	inv[0] = inv[1] = 1;
	mu[1] = phi[1] = 1;
	for (int i = 2; i <= N; i++) inv[i] = 1ll * (P - P / i) * inv[P % i] % P;
	for (int i = 2; i <= N; i++){
		if (!isn[i]) p[++pi] = i,mu[i] = P - 1,phi[i] = i - 1;
		for (int j = 1; j <= pi && i * p[j] <= N; j++){
			isn[i * p[j]] = true;
			if (i % p[j] == 0){
				phi[i * p[j]] = p[j] * phi[i];
				break;
			}
			phi[i * p[j]] = (p[j] - 1) * phi[i];
			mu[i * p[j]] = (P - mu[i]) % P;
		}
	}
	for (int i = 1; i <= N; i++)
		for (int j = i; j <= N; j += i)
			F[j] = (F[j] + 1ll * i * inv[phi[i]] % P * mu[j / i] % P) % P;
	for (int d = 1; d <= N; d++){
		int len = N / d + 1;
		G[d] = new int[len]; G[d][0] = 0;
		for (int i = 1; i < len; i++)
			G[d][i] = (G[d][i - 1] + phi[i * d]) % P;
	}
	for (int x = 1; x < B; x++)
		for (int y = x; y < B; y++){
			int len = N / y + 1;
			T[x][y] = new int[len]; T[x][y][0] = 0;
			for (int i = 1; i < len; i++)
				T[x][y][i] = (T[x][y][i - 1] + 1ll * F[i] * G[i][x] % P * G[i][y] % P) % P;
		}
}
void work(){
	int n = read(),m = read(),ans = 0;
	if (n > m) swap(n,m);
	int E = m / B;
	for (int i = 1; i <= E; i++)
		ans = (ans + 1ll * F[i] * G[i][n / i] % P * G[i][m / i] % P) % P;
	for (int i = E + 1,nxt; i <= n; i = nxt + 1){
		nxt = min(n / (n / i),m / (m / i));
		ans = ((ans + T[n / i][m / i][nxt]) % P + P - T[n / i][m / i][i - 1]) % P;
	}
	printf("%d\n",ans);
}
int main(){
	init();
	int T = read();
	while (T--) work();
	return 0;
}