题目链接
https://www.luogu.org/problemnew/show/P4240
题解
容易发现
φ(ij)=φ(i)φ(j)gcd(i,j)φ(gcd(i,j)) \varphi(ij)=\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))} φ(ij)=φ(gcd(i,j))φ(i)φ(j)gcd(i,j)
因此可以反演得到
∑T=1min(n,m)∑i=1⌊n/k⌋φ(ik)∑i=1⌊m/k⌋φ(ik)∑d∣Tdφ(d)μ(Td) \sum_{T=1}^{\min(n,m)} \sum_{i=1}^{\lfloor n/k\rfloor}\varphi(ik)\sum_{i=1}^{\lfloor m/k\rfloor }\varphi(ik)\sum_{d|T}\frac{d}{\varphi(d)}\mu(\frac{T}{d}) T=1∑min(n,m)i=1∑⌊n/k⌋φ(ik)i=1∑⌊m/k⌋φ(ik)d∣T∑φ(d)dμ(dT)
可以设
g(a,b)=∑i=1aφ(ib) g(a,b)=\sum_{i=1}^a \varphi(ib) g(a,b)=i=1∑aφ(ib)
这个显然可以O(nlogn)O(n\log n)O(nlogn)预处理。
同理,设
f(T)=∑d∣Tdφ(d)μ(Td) f(T)=\sum_{d|T}\frac{d}{\varphi(d)}\mu(\frac{T}{d}) f(T)=d∣T∑φ(d)dμ(dT)
这个也可以O(nlogn)O(n\log n)O(nlogn)预处理。
那么答案就是
∑T=1min(n,m)g(⌊nT⌋,T)g(⌊mT⌋,T)f(T) \sum_{T=1}^{\min(n,m)}g(\lfloor\frac{n}{T}\rfloor,T)g(\lfloor\frac{m}{T}\rfloor,T)f(T) T=1∑min(n,m)g(⌊Tn⌋,T)g(⌊Tm⌋,T)f(T)
钦定一个BBB,对于i,j≤Bi,j\leq Bi,j≤B,设
H(i,j,k)=∑T=1kg(i,T)g(j,T)f(T) H(i,j,k)=\sum_{T=1}^k g(i,T)g(j,T)f(T) H(i,j,k)=T=1∑kg(i,T)g(j,T)f(T)
这个可以O(nB2)O(nB^2)O(nB2)预处理。
考虑处理询问,对于⌊nT⌋,⌊mT⌋≤B\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor\leq B⌊Tn⌋,⌊Tm⌋≤B,整除分块,用预处理出的H(⌊nT⌋,⌊mT⌋,r)−H(⌊nT⌋,⌊mT⌋,l−1)H(\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor,r)-H(\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor,l-1)H(⌊Tn⌋,⌊Tm⌋,r)−H(⌊Tn⌋,⌊Tm⌋,l−1)计入答案;否则,T≤nBT\leq \frac{n}{B}T≤Bn,直接枚举即可。
时间复杂度O(nlogn+nB2+T(n+nB)O(n\log n+nB^2+T(\sqrt{n}+\frac{n}{B})O(nlogn+nB2+T(n+Bn),观察到B=T1/3B=T^{1/3}B=T1/3最优。
代码
#include <cstdio>
#include <vector>
#include <algorithm>
int read()
{
int x=0,f=1;
char ch=getchar();
while((ch<'0')||(ch>'9'))
{
if(ch=='-')
{
f=-f;
}
ch=getchar();
}
while((ch>='0')&&(ch<='9'))
{
x=x*10+ch-'0';
ch=getchar();
}
return x*f;
}
const int maxn=100000;
const int maxm=65;
const int mod=998244353;
int p[maxn+10],prime[maxn+10],cnt,mu[maxn+10],phi[maxn+10],inv[maxn+10],f[maxn+10];
std::vector<int> g[maxn+10],h[maxm+2][maxm+2];
int getprime()
{
inv[0]=inv[1]=1;
for(int i=2; i<=maxn; ++i)
{
inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
p[1]=mu[1]=phi[1]=1;
for(int i=2; i<=maxn; ++i)
{
if(!p[i])
{
prime[++cnt]=i;
mu[i]=mod-1;
phi[i]=i-1;
}
for(int j=1; (j<=cnt)&&(i*prime[j]<=maxn); ++j)
{
int x=i*prime[j];
p[x]=1;
if(i%prime[j]==0)
{
mu[x]=0;
phi[x]=phi[i]*prime[j];
break;
}
mu[x]=mod-mu[i];
phi[x]=phi[i]*(prime[j]-1);
}
}
for(int i=1; i<=maxn; ++i)
{
int v=1ll*i*inv[phi[i]]%mod;
for(int j=1; j<=maxn/i; ++j)
{
f[i*j]=(f[i*j]+1ll*mu[j]*v)%mod;
}
}
for(int i=1; i<=maxn; ++i)
{
g[i].push_back(0);
for(int j=1; j<=maxn/i; ++j)
{
int k=((i==1)?0:g[i-1][j])+phi[i*j];
if(k>=mod)
{
k-=mod;
}
g[i].push_back(k);
}
}
for(int i=1; i<=maxm; ++i)
{
for(int j=1; j<=maxm; ++j)
{
h[i][j].push_back(0);
for(int k=1; k<=std::min(maxn/i,maxn/j); ++k)
{
h[i][j].push_back((h[i][j][k-1]+1ll*g[i][k]*g[j][k]%mod*f[k])%mod);
}
}
}
return 0;
}
int T,n,m;
int main()
{
getprime();
T=read();
while(T--)
{
n=read();
m=read();
int ans=0,k=std::max(n,m)/maxm;
for(int i=1; i<=k; ++i)
{
ans=(ans+1ll*g[n/i][i]*g[m/i][i]%mod*f[i])%mod;
}
for(int l=k+1,r; l<=std::min(n,m); l=r+1)
{
r=std::min(n/(n/l),m/(m/l));
ans+=h[n/l][m/l][r]-h[n/l][m/l][l-1];
if(ans>=mod)
{
ans-=mod;
}
if(ans<0)
{
ans+=mod;
}
}
printf("%d\n",ans);
}
return 0;
}