Description
共\(T(T\leq5\times10^4)\)组数据。给出\(n,m(n,m\leq5\times10^4)\)，求$$\sum_{i=1}n\sum_{j=1}m\sigma_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)\)。
于是就可以正常推导啦。
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(n\sqrt 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\n",ans);
}
return 0;
}
```\]