注意到：
\[\varphi\left(\dfrac{\mathrm{lcm}(i,j)}{\gcd(i,j)}\right)=\varphi\left(\dfrac{ij}{\gcd^2(i,j)}\right)=\varphi\left(\dfrac{i}{\gcd(i,j)}\right)\varphi\left(\dfrac{j}{\gcd(i,j)}\right)
\]
莫反可以得到：
\[\sum_{T=1}^n\sum_{d\mid d}\mu (d)\sum_{i=1}^{\left\lfloor n/T\right\rfloor}\varphi(id)\sum_{j=1}^{\left\lfloor m/T\right\rfloor}\varphi(jd)
\]
设 \(f(d,i)=\sum_{j=1}^i \varphi(jd)\)，可以 \(O(n\log n)\) 处理。
注意到不能直接整除分块，考虑分治，设阈值 \(B\)。
设 \(g(i,j,k)=\sum_{T=1}^k \sum_{d\mid T} \mu(d) f(d,i) f(d,j)\)，其中 \(j\le i\le B\)，状态数 \(O(nB)\)，这样对于 \(\left\lfloor\dfrac{n}{T}\right\rfloor\le B\) 的情况可以整除分块，区间 \([l,r]\) 的答案是 \(g\left(\left\lfloor\dfrac{n}{l}\right\rfloor,\left\lfloor\dfrac{m}{l}\right\rfloor,r\right)-g\left(\left\lfloor\dfrac{n}{l}\right\rfloor,\left\lfloor\dfrac{m}{l}\right\rfloor,l-1\right)\)，其余 \(\left\lfloor\dfrac{n}{T}\right\rfloor>B\) 的情况暴力计算。
总复杂度 \(O\left(n\log n+nB\log n+T\sqrt{n}+T\dfrac{n}{B}\log n\right)\)，取 \(B=170\) 即可。
点击查看代码
int pr[maxn];
bool vis[maxn];
int mu[maxn],phi[maxn];
vector<int> f[maxn],g[maxB][maxB];
inline void linear_sieve(){
mu[1]=1,phi[1]=1;
for(int i=2;i<=lim;++i){
if(!vis[i]) pr[++pr[0]]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=pr[0]&&i*pr[j]<=lim;++j){
vis[i*pr[j]]=1;
mu[i*pr[j]]=-mu[i],phi[i*pr[j]]=phi[i]*phi[pr[j]];
if(i%pr[j]==0){
mu[i*pr[j]]=0,phi[i*pr[j]]=phi[i]*pr[j];
break;
}
}
}
for(int d=1;d<=lim;++d){
f[d].resize(lim/d+1);
f[d][0]=0;
for(int i=1;i<=lim/d;++i){
f[d][i]=(f[d][i-1]+phi[i*d])%mod;
}
}
for(int j=1;j<=B;++j){
for(int k=1;k<=j;++k){
g[j][k].resize(lim/j+1);
for(int d=1;d*j<=lim;++d){
if(!mu[d]) continue;
for(int i=d;i*j<=lim;i+=d){
g[j][k][i]=(g[j][k][i]+1ll*(mu[d]+mod)*f[d][j]%mod*f[d][k]%mod)%mod;
}
}
for(int i=1;i*j<=lim;++i){
g[j][k][i]=(g[j][k][i-1]+g[j][k][i])%mod;
}
}
}
}
int t;
int n,m;
int ans;
int main(){
linear_sieve();
t=read();
while(t--){
n=read(),m=read();
ans=0;
for(int d=1;n/d>B;++d){
for(int i=d;n/i>B;i+=d){
if(!mu[d]) continue;
ans=(ans+1ll*(mu[d]+mod)%mod*f[d][n/i]%mod*f[d][m/i]%mod)%mod;
}
}
for(int l=1,r;l<=m;l=r+1){
r=min(n/(n/l),m/(m/l));
if(n/l<=B) ans=(ans+(g[n/l][m/l][r]-g[n/l][m/l][l-1]+mod)%mod)%mod;
}
printf("%d\n",ans);
}
return 0;
}