题目
题目链接：https://www.luogu.com.cn/problem/P4240
(Q) 次询问，每次询问给出 (n,m)，求
[left(sum^{n}_{i=1}sum^{m}_{j=1}varphi(ij)
ight)mod 998244353
]
(n,mleq 10^5)，(Qleq 10^4)。
思路
[varphi(i)varphi(j)=iprod_{p|i}frac{p-1}{p} imes jprod_{p|j}frac{p-1}{p}
]
[varphi(i)varphi(j)=ij imes prod_{p|ij}frac{p-1}{p} imesprod_{p|gcd(i,j)}frac{p-1}{p}
]
[varphi(i)varphi(j)gcd(i,j)=ij imes prod_{p|ij}frac{p-1}{p} imes 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) imes[gcd(i,j)=d]
]
[=sum^{n}_{d=1}frac{d}{varphi(d)}sum^{lfloorfrac{n}{d}
floor}_{k=1}mu(k)sum^{lfloorfrac{n}{dk}
floor}_{i=1}sum^{lfloorfrac{m}{dk}
floor}_{j=1}varphi(idk)varphi(jdk)
]
[=sum^{n}_{k=1}left(sum_{d|k}frac{dmu(frac{k}{d})}{varphi(d)} imes sum^{lfloorfrac{n}{k}
floor}_{i=1}varphi(ik) imes sum^{lfloorfrac{m}{k}
floor}_{j=1}varphi(jk)
ight)
]
令 (f(i)=sum_{d|i}frac{dmu(frac{i}{d})}{varphi(d)})，(g(k,i)=sum^{i}_{j=1}varphi(jk))。前者枚举倍数可以 (O(nlog n)) 预处理出来，后者由于 (k imes ileq n)，也可以在 (O(nlog n)) 内预处理。
然后原式
[=sum^{n}_{k=1}f(k) imes g(k,lfloorfrac{n}{k}
floor) imes g(k,lfloorfrac{m}{k}
floor)
]
整除分块，令 (h(a,b,k)=sum^{n}_{k=1}f(k) imes g(k,a) imes g(k,b))。当然这个不可能全部预处理出来，考虑平衡规划。
设定一个阈值 (B)。预处理出 (max(a,b)leq B) 时所有的 (h)，当 (max(a,b)>B) 时，说明 (kleq lfloorfrac{n}{B}
floor)。直接暴力求即可。
时间复杂度 (O(nlog 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<<"
";
}
return 0;
}