洛谷P4240 毒瘤之神的考验

传送门

这种题第一步只能是用一些经典推法推柿子。但是这个 \(\varphi (ij)\) 显得无从下手。

所以这里必须知道一个经典结论：\(\varphi (ij)=\dfrac{\varphi (i)\varphi (j)\gcd(i,j)}{\varphi (gcd(i,j))}\)。

然后是一些经典套路。（下文钦定 \(n\leqslant m\)）

\(\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi (ij) &=\sum_{i=1}^{n}\sum_{j=1}^{m} \dfrac{\varphi (i)\varphi (j)\gcd(i,j)}{\varphi (gcd(i,j))} \\&=\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor }\sum_{j=1}^{\left\lfloor \frac{m}{d} \right\rfloor}\varphi (id)\varphi (jd)\,[ \gcd(i,j)= 1 ] \\&=\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{t=1}^{\left\lfloor \frac{n}{d} \right\rfloor }\mu (t)\sum_{i=1}^{\left \lfloor \frac{n}{dt} \right \rfloor } \sum_{j=1}^{\left \lfloor \frac{m}{dt} \right \rfloor } \varphi (idt)\,\varphi (jdt) \\&=\sum_{T=1}^{n}\sum_{d|T}^{}\frac{d\cdot \mu (\frac{T}{d} )}{\varphi (d)} \sum_{i=1}^{\left \lfloor \frac{n}{T} \right \rfloor }\varphi(iT)\sum_{j=1}^{\left \lfloor \frac{m}{T} \right \rfloor }\varphi(jT)\quad(令T=dt)\end{aligned}\)

讲个笑话：笔者只能算到这一步了。

观察最后的柿子，考虑函数封装起来。

设 $\displaystyle f(x)=\sum_{d|x}^{}\frac{d\cdot \mu (\frac{x}{d} )}{\varphi (d)}\ \ (x\leqslant n)\ $ 和 \(\ \displaystyle g(x,y)=\sum_{i=1}^{y}\varphi(xi)\ \ (xy\leqslant n)\)

这两个函数的预处理容易做到调和级数 \(O(n \ln n)\) 。

然后整个柿子表示为

\(\begin{aligned} S(y,z,x)&=\sum_{T=1}^{x}\sum_{d|T}^{}\frac{d\cdot \mu (\frac{T}{d} )}{\varphi (d)} \sum_{i=1}^{y}\varphi(iT)\sum_{j=1}^{z}\varphi(jT) \\&=S(y,z,x-1)+f(x)\cdot g(x,y) \cdot g(x,z)\end{aligned}\)

那么对于连续的一段区间满足 $\left \lfloor \frac{n}{l} \right \rfloor = \left \lfloor \frac{n}{r} \right \rfloor $ 且 $\left \lfloor \frac{m}{l} \right \rfloor = \left \lfloor \frac{m}{r} \right \rfloor \ $, 我们就要求出 $S(\left \lfloor \frac{n}{l} \right \rfloor ,\left \lfloor \frac{m}{l} \right \rfloor ,r) - S(\left \lfloor \frac{n}{l} \right \rfloor ,\left \lfloor \frac{m}{l} \right \rfloor ,l-1) $。

观察到 \(r\) 越大，即要算的 \(S(y,z,x)\) 中的 \(x\) 越大，那么显然 \(\left \lfloor \frac{n}{r} \right \rfloor ,\left \lfloor \frac{m}{r} \right \rfloor\) 也就是 \(y,z\) 就越小。
所以当 \(\left \lfloor \frac{m}{r} \right \rfloor\) 小到一定程度，也就是 \(l\) 和 \(r\) 大到一定程度后其实就可以通过预处理代替暴力减小复杂度。

所以设一个阙值 \(B\) ，表示我们考虑预处理出全部的 \(S(1,1,1)\) 至 \(S(B,B,lim)\ (lim=\left \lfloor \frac{n}{B} \right \rfloor)\) 。在求答案的时候，如果我们的 \(\left \lfloor \frac{m}{r} \right \rfloor \leqslant B\) ，那么考虑直接调用预处理数组 \(O(1)\) 求贡献。否则就有 \(r\leqslant \left \lfloor \frac{m}{B} \right \rfloor\) ，这时候直接枚举 \(T\) 暴力计算。

容易求出总复杂度为 \(O(n \log n+nB+T(\sqrt{n}+\dfrac{n}{B}))\) 。 \(B\) 取 \(\sqrt{T}\) 即 \(100\) 是最优的。(傻瓜注释：这里的 \(T\) 是数据组数，非上文柿子中的 \(T\) 。)

看到这里一定会有人问预处理 \(S\) 为什么是 \(O(nB)\) 而不是 \(O(nB^2)\) 。

注意到 \(B\cdot lim \leqslant n\)，所以复杂度为：

\(\displaystyle O(\sum_{i=1}^{B}\sum_{j=i}^{B}\frac{n}{j})=O(\sum_{j=1}^{B}j\cdot \frac{n}{j})=O(nB)\) 。

code

//writer:Oier_szc

#include <bits/stdc++.h>
//#include <windows.h>
#define ED cerr<<endl;
#define TS cerr<<"I AK IOI"<<endl;
#define cr(x) cerr<<x<<endl;
#define cr2(x,y) cerr<<x<<" "<<y<<endl;
#define cr3(x,y,z) cerr<<x<<" "<<y<<" "<<z<<endl;
#define cr4(x,y,z,w) cerr<<x<<" "<<y<<" "<<z<<" "<<w<<endl;
#define pii pair<int,int>
#define fi first
#define se second
#define int long long
//#define ull unsigned long long
using namespace std;
const int N=1e5+5,B=100,INF=2e9,mod=998244353;
int T,n,m;
int primes[N],st[N],mu[N],phi[N],len=0;
int f[N],inv[N];
vector<int> s[N],g[N];
vector<int> t[B+5][B+5];

int qpow(int a,int b) {
	int res=1;
	while(b) {
		if(b&1) res=res*a%mod;
		a=a*a%mod,b>>=1;
	}
	return res;
}

void init(int up) {
	mu[1]=phi[1]=1;
	for(int i=2;i<=up;++i) {
		if(!st[i]) primes[++len]=i,mu[i]=-1,phi[i]=i-1;
		for(int j=1;j<=len&&primes[j]<=up/i;++j) {
			st[i*primes[j]]=1;
			if(i%primes[j]==0) {
				mu[i*primes[j]]=0;
				phi[i*primes[j]]=phi[i]*primes[j];
				break;
			}
			else {
				mu[i*primes[j]]=-mu[i];
				phi[i*primes[j]]=phi[i]*(primes[j]-1);
			}
		}
	}
	inv[1]=1;
	for(int i=2;i<=up;++i) {
		inv[i]=(mod-mod/i)*inv[mod%i]%mod;
	}
	for(int i=1;i<=up;++i) {
		g[i].resize(up/i+2);
		for(int j=i,k=1;j<=up;j+=i,++k) {
			s[j].emplace_back(i);
			f[j]=(f[j]+i*mu[k]%mod*inv[phi[i]]%mod)%mod;
			g[i][k]=(g[i][k-1]+phi[j])%mod;
		}
	}
	for(int i=1;i<=B;++i) {
		for(int j=i;j<=B;++j) {
			t[i][j].resize(up/j+2);
			for(int k=1;k<=up/j;++k) {
				t[i][j][k]=(t[i][j][k-1]+f[k]*g[k][i]%mod*g[k][j]%mod)%mod;
			}
		}
	}
}

int sol(int n,int m) {
	if(n>m) swap(n,m);
	int res=0;
	for(int i=1;i<=m/B;++i) res=(res+f[i]*g[i][n/i]%mod*g[i][m/i]%mod)%mod;
	for(int l=m/B+1,r;l<=n;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		res=(res+t[n/l][m/l][r]-t[n/l][m/l][l-1]+mod)%mod;
	}
	return res;
}

signed main()
{
	init(1e5);
	scanf("%lld",&T);
	while(T--) {
		scanf("%lld%lld",&n,&m);
		printf("%lld\n",sol(n,m));
	}
	return 0;
}