【洛谷P4240】毒瘤之神的考验

题目

题目链接：https://www.luogu.com.cn/problem/P4240
\(Q\) 次询问，每次询问给出 \(n,m\)，求

\[\left(\sum^{n}_{i=1}\sum^{m}_{j=1}\varphi(ij)\right)\bmod 998244353 \]

\(n,m\leq 10^5\)，\(Q\leq 10^4\)。

思路

\[\varphi(i)\varphi(j)=i\prod_{p|i}\frac{p-1}{p}\times j\prod_{p|j}\frac{p-1}{p} \]

\[\varphi(i)\varphi(j)=ij\times \prod_{p|ij}\frac{p-1}{p}\times\prod_{p|\gcd(i,j)}\frac{p-1}{p} \]

\[\varphi(i)\varphi(j)\gcd(i,j)=ij\times \prod_{p|ij}\frac{p-1}{p}\times \gcd(i,j)\prod_{p|\gcd(i,j)}\frac{p-1}{p} \]

\[\varphi(ij)=\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))} \]

所以

\[\sum^{n}_{i=1}\sum^{m}_{j=1}\varphi(ij)=\sum^{n}_{i=1}\sum^{n}_{j=1}\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))} \]

\[=\sum^{n}_{d=1}\frac{d}{\varphi(d)}\sum^{n}_{i=1}\sum^{n}_{j=1}\varphi(i)\varphi(j)\times[\gcd(i,j)=d] \]

\[=\sum^{n}_{d=1}\frac{d}{\varphi(d)}\sum^{\lfloor\frac{n}{d}\rfloor}_{k=1}\mu(k)\sum^{\lfloor\frac{n}{dk}\rfloor}_{i=1}\sum^{\lfloor\frac{m}{dk}\rfloor}_{j=1}\varphi(idk)\varphi(jdk) \]

\[=\sum^{n}_{k=1}\left(\sum_{d|k}\frac{d\mu(\frac{k}{d})}{\varphi(d)}\times \sum^{\lfloor\frac{n}{k}\rfloor}_{i=1}\varphi(ik)\times \sum^{\lfloor\frac{m}{k}\rfloor}_{j=1}\varphi(jk)\right) \]

令 \(f(i)=\sum_{d|i}\frac{d\mu(\frac{i}{d})}{\varphi(d)}\)，\(g(k,i)=\sum^{i}_{j=1}\varphi(jk)\)。前者枚举倍数可以 \(O(n\log n)\) 预处理出来，后者由于 \(k\times i\leq n\)，也可以在 \(O(n\log n)\) 内预处理。
然后原式

\[=\sum^{n}_{k=1}f(k)\times g(k,\lfloor\frac{n}{k}\rfloor)\times g(k,\lfloor\frac{m}{k}\rfloor) \]

整除分块，令 \(h(a,b,k)=\sum^{n}_{k=1}f(k)\times g(k,a)\times g(k,b)\)。当然这个不可能全部预处理出来，考虑平衡规划。
设定一个阈值 \(B\)。预处理出 \(\max(a,b)\leq B\) 时所有的 \(h\)，当 \(\max(a,b)>B\) 时，说明 \(k\leq \lfloor\frac{n}{B}\rfloor\)。直接暴力求即可。
时间复杂度 \(O(n\log n+nB^2+Q(\sqrt{n}+\frac{n}{B}))\)。空间复杂度 \(O(nB^2)\)。取 \(B=20\) 即可。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=100010,B=20,MOD=998244353;
int Q,n,m,mu[N],phi[N],prm[N];
bool v[N];
ll ans,f[N],h[B+2][B+2][N];
vector<ll> g[N];

void findprm(int n)
{
	mu[1]=phi[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,phi[i]=i-1,mu[i]=-1;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1;
			mu[i*prm[j]]=-mu[i]; phi[i*prm[j]]=phi[i]*(prm[j]-1);
			if (!(i%prm[j]))
			{
				mu[i*prm[j]]=0; phi[i*prm[j]]=phi[i]*prm[j];
				break;
			}
		}
	}
}

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

int main()
{
	findprm(N-1);
	for (int i=1;i<N;i++)
	{
		ll inv=fpow(phi[i],MOD-2);		
		for (int j=i;j<N;j+=i)
			f[j]=(f[j]+1LL*i*mu[j/i]*inv)%MOD;
		g[i].push_back(0LL);
		for (int j=1;i*j<N;j++)
			g[i].push_back((g[i][j-1]+phi[i*j])%MOD);
	}
	for (int i=1;i<=B;i++)
		for (int j=1;j<=B;j++)
			for (int k=1;k<N;k++)
				if (i*k<N && j*k<N)
					h[i][j][k]=(h[i][j][k-1]+f[k]*g[k][i]%MOD*g[k][j])%MOD;
	scanf("%d",&Q);
	while (Q--)
	{
		scanf("%d%d",&n,&m);
		if (n>m) swap(n,m);
		ans=0;
		for (int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			if (m/l<=B) ans=(ans+h[n/l][m/l][r]-h[n/l][m/l][l-1])%MOD;
			else
			{
				for (int i=l;i<=r;i++)
					ans=(ans+f[i]*g[i][n/i]%MOD*g[i][m/i])%MOD;
			}
		}
		cout<<(ans+MOD)%MOD<<"\n";
	}
	return 0;
}