Description
共(T(Tleq5 imes10^4))组数据。给出(n,m(n,mleq5 imes10^4))，求$$sum_{i=1}nsum_{j=1}msigma_0(ij)$$
Solution
首先有结论：(sigma_0(xy)=sum_{d_1|x}sum_{d_2|y}[gcd(d_1,d_2)=1])。下面先证明一下这个结论。
将(x,y)分解质因数，得到(x=prod_{i=1}^kp_i^{a_i})，(y=prod_{i=1}^kp_i^{b_i})。那么若(gcd(d_1,d _2)=1)，则说明对于一个质因数(p_i)，其只存在于(d_1,d_2)之中的零或一个中。若存在于(d_1)中则有(a_i)种可能性，若存在于(d_2)中则有(b_i)种可能性，再加上都不存在的一种共计(a_i+b_i+1)种。那么共有(prod_{i=1}^k(a_i+b_i+1))对互质的((d_1,d_2))，而这个数恰等于(sigma_0(xy))。
于是就可以正常推导啦。
[egin{align*}
ans &= sum_{i=1}^n sum_{j=1}^m sum_{d_1|i} sum_{d_2|j}[gcd(d_1,d_2)=1] \
&= sum_{d_1=1}^{+∞} sum_{d_2=1}^{+∞} sum_{d_1|i}^n sum_{d_2|j}^m [gcd(d_1,d_2)=1] \
&= sum_{d_1=1}^{+∞} sum_{d_2=1}^{+∞} ⌊frac{n}{d_1}⌋ ⌊frac{m}{d_2}⌋ sum_{d_3|gcd(d_1,d_2)} mu(d_3) \
&= sum_{d_3=1}^{+∞} mu(d_3) sum_{d_3|d_1} sum_{d_3|d_2} ⌊frac{n}{d_1}⌋ ⌊frac{m}{d_2}⌋ \
&= sum_{d_3=1}^{+∞} mu(d_3) sum_{t_1=1}^{+∞} ⌊frac{⌊frac{n}{d_3}⌋}{t_1}⌋ sum_{t_2=1}^{+∞} ⌊frac{⌊frac{m}{d_3}⌋}{t_2}⌋ \
end{align*}$$ $O(nsqrt n)$预处理出$f(x)=sum_{i=1}^x ⌊dfrac{x}{i}⌋$，就有$$ans=sum_{d=1}^{+∞} mu(d)f(⌊dfrac{n}{d}⌋)f(⌊dfrac{m}{d}⌋)$$ 就可以在$O(sqrt n)$的时间内求解啦。
##Code
```cpp
//[SDOI2015]约数个数和
#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long lint;
inline char gc()
{
static char now[1<<16],*s,*t;
if(s==t) {t=(s=now)+fread(now,1,1<<16,stdin); if(s==t) return EOF;}
return *s++;
}
inline int read()
{
int x=0; char ch=gc();
while(ch<'0'||'9'<ch) ch=gc();
while('0'<=ch&&ch<='9') x=x*10+ch-'0',ch=gc();
return x;
}
const int N=5e4+10;
int pCnt,pr[N]; bool pNot[N];
int mu[N],pre[N]; int f[N];
void init(int n)
{
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!pNot[i]) pr[++pCnt]=i,mu[i]=-1;
for(int j=1;j<=pCnt;j++)
{
int x=i*pr[j]; if(x>n) break;
pNot[x]=true;
if(i%pr[j]) mu[x]=-mu[i];
else {mu[x]=0; break;}
}
}
for(int i=1;i<=n;i++) pre[i]=pre[i-1]+mu[i];
for(int x=1;x<=n;x++)
for(int L=1,R;L<=x;L=R+1)
{
int v=x/L; R=x/v;
f[x]+=(R-L+1)*v;
}
}
int main()
{
int task=read(); init(5e4);
while(task--)
{
int n=read(),m=read();
if(n>m) swap(n,m);
lint ans=0;
for(int L=1,R;L<=n;L=R+1)
{
int v1=n/L,v2=m/L; R=min(n/v1,m/v2);
ans+=1LL*(pre[R]-pre[L-1])*f[v1]*f[v2];
}
printf("%lld
",ans);
}
return 0;
}
```]