( ext{Description})
( ext{Solution})
首先有：
[sigma(x imes y)=sum_{i|x}sum_{j|y} [gcd(i,j)=1]
]
具体证明戳这。
柿子变成了：
[sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j} [gcd(x,y)=1]
]
而常见形式是：
[sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=1]
]
尝试将 (x,y) 提出去，向常见形式靠拢。
其实可以发现 (sum_{i=1}^nsum_{j=1}^m) 就是个限制条件，它决定了次数。
[sum_{x=1}^nsum_{y=1}^m[gcd(x,y)=1] imes left lfloorfrac{n}{x}
ight
floor imes left lfloorfrac{m}{y}
ight
floor
]
设，
[g(k)=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=k] imes left lfloorfrac{n}{i}
ight
floor imes left lfloorfrac{m}{j}
ight
floor
]
[f(k)=sum_{k|d}g(d)
]
所以，
[ ext{Ans}=g(1)=sum_{i=1}^{min{n,m}}mu(i) imes f(i)
]
我们需要快速求出 (f(i))（根据经验，应该是可以整除分块的）。
不过这次似乎无法像 (sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=1],sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=1] imes i imes j) 这样的柿子猛拆二重循环，因为有 (left lfloorfrac{n}{i} ight floor imes left lfloorfrac{m}{j} ight floor) 碍事的。
尝试直接化柿子：
[f(k)=sum_{i=1}^nsum_{j=1}^m[k|gcd(i,j)] imes left lfloorfrac{n}{i}
ight
floor imes left lfloorfrac{m}{j}
ight
floor
]
显然我们无法简化 (i,j)，显然不能枚举 (gcd)。
其实可以将 (k) 提出来，这样 (gcd) 肯定是 (k) 的倍数：
[=sum_{i=1}^{left lfloorfrac{n}{k}
ight
floor}sum_{j=1}^{left lfloorfrac{m}{k}
ight
floor}left lfloorfrac{n}{i imes k}
ight
floor imes left lfloorfrac{m}{j imes k}
ight
floor
]
预处理 (sum_{i=1}^n left lfloorfrac{n}{i} ight floor) 就行了。
时间复杂度 (mathcal O((T+n) imes sqrt n))。
( ext{Code})
#include <cstdio>
#define rep(i,_l,_r) for(register signed i=(_l),_end=(_r);i<=_end;++i)
#define fep(i,_l,_r) for(register signed i=(_l),_end=(_r);i>=_end;--i)
#define erep(i,u) for(signed i=head[u],v=to[i];i;i=nxt[i],v=to[i])
#define efep(i,u) for(signed i=Head[u],v=to[i];i;i=nxt[i],v=to[i])
#define print(x,y) write(x),putchar(y)
template <class T> inline T read(const T sample) {
T x=0; int f=1; char s;
while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
return x*f;
}
template <class T> inline void write(const T x) {
if(x<0) return (void) (putchar('-'),write(-x));
if(x>9) write(x/10);
putchar(x%10^48);
}
template <class T> inline T Max(const T x,const T y) {if(x>y) return x; return y;}
template <class T> inline T Min(const T x,const T y) {if(x<y) return x; return y;}
template <class T> inline T fab(const T x) {return x>0?x:-x;}
template <class T> inline T gcd(const T x,const T y) {return y?gcd(y,x%y):x;}
template <class T> inline T lcm(const T x,const T y) {return x/gcd(x,y)*y;}
typedef long long ll;
const int maxn=5e4+5;
int n,m,u[maxn],p[maxn],pc;
bool is[maxn];
ll f[maxn];
void init() {
u[1]=1;
rep(i,2,maxn-5) {
if(!is[i]) p[++pc]=i,u[i]=-1;
rep(j,1,pc) {
if(p[j]*i>maxn-5) break;
is[p[j]*i]=1;
if(i%p[j]==0) break;
u[i*p[j]]=-u[i];
}
}
rep(i,2,maxn-5) u[i]+=u[i-1];
int r; ll ret;
rep(i,1,maxn-5) {
ret=0;
for(int l=1;l<=i;l=r+1) r=i/(i/l),ret+=1ll*(r-l+1)*(i/l);
f[i]=ret;
}
}
ll Query() {
int lim=Min(n,m),r; ll ret=0;
for(int l=1;l<=lim;l=r+1) {
r=Min(n/(n/l),m/(m/l));
ret+=1ll*(u[r]-u[l-1])*f[n/l]*f[m/l];
}
return ret;
}
int main() {
init();
for(int T=read(9);T;--T) {
n=read(9),m=read(9);
print(Query(),'
');
}
return 0;
}