显然题目的答案就是$$ Ans=sum_{i=1}^Nsum_{j=1}^M[gcd(i,j)=prime]$$
我们先设设F(n)表示满足(gcd(i,j)\%t=0)的数对个数，f(t)表示满足(gcd(i,j)=t)的数对个数
$$f(t)=sum_{i=1}^Nsum_{j=1}^M[gcd(i,j)=t]$$
$$F(n)=sum_{n|t}lfloor frac{N}{n}
floor lfloor frac{M}{n}
floor$$
那么根据莫比乌斯反演的第二种形式珂以得到
$$f(n)=sum_{n|t}mu(lfloor frac{t}{n}
floor)F(t)$$
所以答案珂以变形为：
$$Ans=sum_{p in prime}sum_{i=1}^Nsum_{j-1}^M[gcd(i,j)=p)$$
$$=sum_{p in prime}f(p) qquad qquad quad$$
$$=sum_{p in prime}sum_{p|t}mu(lfloor frac{t}{p}
floor)F(t)$$
我们不枚举p，我们枚举(lfloor frac{t}{p}
floor)
$$Ans=sum_{p in prime}sum_{d=1}^{Min(frac{N}{p},frac{M}{p})}mu(t)F(tp)$$
$$qquad qquad qquad=sum_{p in prime}sum_{d=1}^{Min(frac{N}{p},frac{M}{p})}mu(t)sum_{n|t}lfloor frac{N}{tp}
floor lfloor frac{M}{tp}
floor$$
我们把tp换成T继续变形
$$Ans=sum_{T=1}^{Min(N,M)}lfloor frac{N}{T}
floor lfloor frac{M}{T}
floor(sum_{p|T,p in prime}mu(frac{T}{p}))$$
这样就珂以用整除分块求了qaq
#include <bits/stdc++.h>
#define N 10000005
#define ll long long
#define getchar nc
using namespace std;
inline char nc(){
static char buf[100000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline int read()
{
register int x=0,f=1;register char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return x*f;
}
inline void write(register ll x)
{
if(!x)putchar('0');if(x<0)x=-x,putchar('-');
static int sta[30];register int tot=0;
while(x)sta[tot++]=x%10,x/=10;
while(tot)putchar(sta[--tot]+48);
}
inline int Min(register int x,register int y)
{
return x<y?x:y;
}
int v[N],miu[N],prim[N],cnt=0,g[N];
ll sum[N];
int main()
{
miu[1]=1;
for(register int i=2;i<=N;++i)
{
if(!v[i])
{
miu[i]=-1;
prim[++cnt]=i;
}
for(register int j=1;j<=cnt&&prim[j]*i<=N;++j)
{
v[i*prim[j]]=1;
if(i%prim[j]==0)
break;
else
miu[prim[j]*i]=-miu[i];
}
}
for(register int i=1;i<=cnt;++i)
for(register int j=1;j*prim[i]<=N;++j)
g[j*prim[i]]+=miu[j];
for(register int i=1;i<=N;++i)
sum[i]=sum[i-1]+(ll)g[i];
int t=read();
while(t--)
{
int n=read(),m=read();
if(n>m)
n^=m^=n^=m;
ll ans=0;
for(register int l=1,r;l<=n;l=r+1)
{
r=Min(n/(n/l),m/(m/l));
ans+=(ll)(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
write(ans),puts("");
}
return 0;
}