不妨设\(n<m\)
首先将原提问形式化，可以得到原式为：
$$\sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{p\in P}[gcd(x,y)==p] $$
将所有的\(x,y\)都除以\(gcd(x,y)\)，则原式可转化为：
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}[gcd(x,y)==1]$$
然后，根据莫比乌斯函数的基本性质，我们有：
$$\sum_{d|n}\mu(d)=[n==1]$$
替换入原式可以得到：
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|gcd(x,y)}\mu(d)$$
将\(d|gcd(x,y)\)一式化简可得：
$$d|gcd(x,y)=d|x\&\&d|y$$
代入原式：
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|x\&\&d|y}\mu(d)$$
将所有的\(x,y\)都除以\(d\)，则原式可转化为：
$$\sum_{x=1}^{\lfloor \frac{n}{pd} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{pd} \rfloor}\sum_{p\in P}\mu(d)$$
由于\(x,y\)对后半部分式子无影响，且\(pd\le n\)原式可变式为：
$$\sum_{p\in P}\sum_{d=1}^{\lfloor \frac{n}{p} \rfloor}\mu(d)\lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor$$
为了减少枚举的数量，不妨设\(q=pd\)，那么我们又可以将原式转化为：
$$\sum_{q=1}^{n}\lfloor \frac{n}{q} \rfloor \lfloor \frac{m}{q} \rfloor \sum_{p\in P\&\&pd==q}\mu(d)$$
又，根据狄利克雷卷积的基本式
$$g=f u \Rightarrow g(x)=\sum_{a b=x} f(a) u(b)$$
我们发现，上式靠右的如下部分可以作出变换：
$$\sum_{p\in P\&\&pd==q}\mu(d)=g(q)=\sum_{pd=q}[p\in P]\mu(d)$$
若令：
$$f(p)=[p\in P]$$
则我们可以发现：
$$ g(q)=f(p) \mu(d)$$
故而这一部分的值可以用狄利克雷卷积来\(O(nln_n)\)预处理。
一看数据范围，完了：\(10^7\)，\(O(nln_n)\)的复杂度显然是不够优秀的。我们需要将它优化到\(O(n)\)
我们考虑对这个函数\(g\)进行线性筛。
我们发现，上式可以进行如下变换：
$$g(q)=\sum_{pd=q\&\&p\in P}\mu(d)=\sum_{p\in P\&\&p|q}\mu(\frac{q}{p})$$
本质上就是枚举其质因子之和并求其除以该质因子的莫比乌斯函数值之和。
我们分三类讨论。
对于\(g(x)\)
如果\(x\in P\)，显然\(g(x)=\mu(1)=1\)
如果\(x=a^vb,a\in P,v>1\)，则显然枚举其他的质数是没有意义的。\(g(x)=\mu(b)\)
如果\(x=ab,a\in P\)，则，\(g(x)\)的值，等价于在计算\(g(b)\)的和的时候每一次计算都多乘上一个\(a\)的贡献，也就是乘上负一。然后将这个值再加上枚举\(a\)的情况，也就是\(\mu(b)\)，即\(g(x)=-g(b)+\mu(b)\)
然后大力上一个线性筛即可。
对于左半部分，上一个数论分块即可。
#include<iostream>
#include<cstdio>
const int N = 10000000+5;
int p[N/10],g[N],mu[N],n,m;
bool ip[N];
void prpr(){
p[0]=0;mu[1]=1;ip[1]=1;
for(int i=2;i<=10000000;++i){
if(!ip[i]){
p[++p[0]]=i;
mu[i]=-1;
g[i]=1;
}
for(int j=1;j<=p[0]&&p[j]*i<=10000000;++j){
ip[i*p[j]]=1;
if(!(i%p[j])){
g[i*p[j]]=mu[i];
mu[i*p[j]]=0;
break;
}else{
mu[i*p[j]]=-mu[i];
g[i*p[j]]=-g[i]+mu[i];
}
}
}
for(int i=2;i<=10000000;++i){
g[i]+=g[i-1];
}
}
long long ans;
void init(){
scanf("%d%d",&n,&m);
n>m?n^=m^=n^=m:0;
ans=0;
int k=0;
for(int i=1;i<=n;i=k+1){
k=std::min(n/(n/i),m/(m/i));
ans+=1ll*(n/i)*(m/i)*(g[k]-g[i-1]);
}
printf("%lld\n",ans);
}
int main(){
prpr();
int T;
scanf("%d",&T);
while(T--){
init();
}
return 0;
}