题目链接

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))}
因此可以反演得到
T=1min(n,m)i=1n/kφ(ik)i=1m/kφ(ik)dTdφ(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})
可以设
g(a,b)=i=1aφ(ib) g(a,b)=\sum_{i=1}^a \varphi(ib)
这个显然可以O(nlogn)O(n\log n)预处理。

同理,设
f(T)=dTdφ(d)μ(Td) f(T)=\sum_{d|T}\frac{d}{\varphi(d)}\mu(\frac{T}{d})
这个也可以O(nlogn)O(n\log n)预处理。

那么答案就是
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)
钦定一个BB,对于i,jBi,j\leq 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)
这个可以O(nB2)O(nB^2)预处理。

考虑处理询问,对于nT,mTB\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor\leq B,整除分块,用预处理出的H(nT,mT,r)H(nT,mT,l1)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)计入答案;否则,TnBT\leq \frac{n}{B},直接枚举即可。

时间复杂度O(nlogn+nB2+T(n+nB)O(n\log n+nB^2+T(\sqrt{n}+\frac{n}{B}),观察到B=T1/3B=T^{1/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;
}