可以得到：
证明考虑 \(\varphi\) 的展开式。
选取中间的式子带进去化简。
于是无脑反演可以得到：
里面有关于 \(T\) 的部分，不能直接预处理某个函数接上整除分块。
设 \(f(n)=\sum_{k\mid n} \mu(n/k)\dfrac{k}{\varphi(k)}\) 暴力卷积 \(O(n\log n)\) 可以预处理。
再设 \(g(n,k)=\sum_{i=1}^n \varphi(ik)\)，由于 \(nk\le 10^5\)，所以合法状态数是 \(O(n\log n)\) 的，可以 vector 开空间预处理，递推式 \(g(n,k)=g(n-1,k)+\varphi(nk)\)。
但是这样必须逐个枚举 \(T\)，但我们希望整除分块。
之后比较优美之处在于，上界 \(\left\lfloor n/T\right\rfloor\) 和枚举的 \(T\) 看起来可以平衡，也就是当 \(T\) 较大时，\(\left\lfloor n/T\right\rfloor\) 较小，可以暴力预处理，而 \(T\) 较小时暴力逐个算就完了。
这时候我们不再关心 \(n,m\) 了，主要是把所有 \((\left\lfloor n/T\right\rfloor,\left\lfloor m/T\right\rfloor)\) 的值都列出来，此时还需要知道 \(T\) 的值，于是设 \(h(n,a,b)=\sum_{i=1}^n f(i)\times g(a,i)\times g(b,i)\)，那么有递推式 \(h(n,a,b)=h(n-1,a,b)+f(n)\times g(a,n)\times g(b,n)\)。
这样就只要设定一个阈值 \(B\)，当 \(\left\lfloor n/T\right\rfloor\le B\) 时使用 \(h\) 计算，具体是形如 \(h(r,\left\lfloor n/l\right\rfloor,\left\lfloor m/l\right\rfloor)-h(l-1,\left\lfloor n/l\right\rfloor,\left\lfloor m/l\right\rfloor)\)，其他枚举 \(T\) 用 \(f,g\) 计算。
点击查看代码
int t;
int n,m,B=40;
int pr[maxn],mu[maxn],phi[maxn],mn[maxn],mnpw[maxn];
bool vis[maxn];
int inv[maxn],f[maxn];
vector<int> g[maxn];
vector<int> h[41][41];
inline void linear_sieve(){
mu[1]=1,phi[1]=1,mn[1]=1,mnpw[1]=1;
for(int i=2;i<=lim;++i){
if(!vis[i]) pr[++pr[0]]=i,mu[i]=-1,phi[i]=i-1,mn[i]=i,mnpw[i]=i;
for(int j=1;j<=pr[0]&&i*pr[j]<=lim;++j){
vis[i*pr[j]]=1;
mu[i*pr[j]]=-mu[i],mn[i*pr[j]]=pr[j];
if(mn[i]==pr[j]) mnpw[i*pr[j]]=mnpw[i]*pr[j];
else mnpw[i*pr[j]]=pr[j];
if(i*pr[j]==mnpw[i*pr[j]]) phi[i*pr[j]]=phi[i]*pr[j];
else phi[i*pr[j]]=phi[i*pr[j]/mnpw[i*pr[j]]]*phi[mnpw[i*pr[j]]];
if(i%pr[j]==0){
mu[i*pr[j]]=0;
break;
}
}
}
inv[1]=1;
for(int i=2;i<=lim;++i) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
for(int i=1;i<=lim;++i){
for(int j=1;i*j<=lim;++j){
int now=(1ll*mu[i]*j*inv[phi[j]]%mod+mod)%mod;
f[i*j]=(f[i*j]+now)%mod;
}
}
for(int i=0;i<=lim;++i) g[0].push_back(0);
for(int i=1;i<=lim;++i){
g[i].push_back(0);
for(int j=1;i*j<=lim;++j){
g[i].push_back(0);
g[i][j]=(g[i-1][j]+phi[i*j])%mod;
}
}
for(int j=1;j<=B;++j){
for(int k=1;k<=B;++k){
h[j][k].push_back(0);
for(int i=1;i<=lim&&i*j<=lim&&i*k<=lim;++i){
h[j][k].push_back(0);
h[j][k][i]=(h[j][k][i-1]+1ll*f[i]*g[j][i]%mod*g[k][i]%mod)%mod;
}
}
}
}
int ans;
int main(){
linear_sieve();
t=read();
while(t--){
n=read(),m=read();
if(n>m) swap(n,m);
ans=0;
for(int i=1;i<=min(m/B,n);++i){
ans=(ans+1ll*f[i]*g[n/i][i]%mod*g[m/i][i]%mod)%mod;
}
for(int l=m/B+1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+(h[n/l][m/l][r]-h[n/l][m/l][l-1]+mod)%mod)%mod;
}
printf("%d\n",ans);
}
return 0;
}