# 题目链接

https://www.luogu.org/problemnew/show/P4240

# 题解

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

$\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)=\sum_{i=1}^a \varphi(ib)$

$f(T)=\sum_{d|T}\frac{d}{\varphi(d)}\mu(\frac{T}{d})$

$\sum_{T=1}^{\min(n,m)}g(\lfloor\frac{n}{T}\rfloor,T)g(\lfloor\frac{m}{T}\rfloor,T)f(T)$

$H(i,j,k)=\sum_{T=1}^k g(i,T)g(j,T)f(T)$

# 代码

#include <cstdio>
#include <vector>
#include <algorithm>

{
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();
while(T--)
{
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;
}