【洛谷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;
}
