题解「Luogu3327 [SDOI2015]约数个数和」
首先有个东西是这题解题的关键：
\[{\rm{d}}(ij)=\sum_{x|i}\sum_{y|i}[{\rm{gcd}}(x,y)=1]
\]
有位dalao的题解证明了这个式子，可以去看看。
然后就可以开始推式子了：
\[\sum_{i=1}^n\sum_{j=1}^md(ij) =\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|i}[{\rm{gcd}}(x,y)=1]=\sum_{x=1}^n\sum_{y=1}^m[{\rm{gcd}}(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor
\]
接下来用莫比乌斯函数替换掉 \([{\rm{gcd}}(x,y)=1]\) ：
\[\sum_{d|x,d|y}\mu(d)=[{\rm{gcd}}(x,y)=1]
\]
代入：
\[\begin{align}
\sum_{x=1}^n\sum_{y=1}^m[{\rm{gcd}}(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor & =\sum_{x=1}^n\sum_{y=1}^m \big(\sum_{d|x,d|y}\mu(d)\big)\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor \\
\end{align}
\]
将这个式子由枚举 \(x,y\) 的形式转变为枚举 \(d\) ，即用 \(xd,yd\) 代替 \(x,y\) ：
\[= \sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{yd}\rfloor
= \sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\rfloor
\]
这里枚举的 \(x,y\) 是上式中的 \(\frac{x}{d},\frac{y}{d}\) 。
设 \(f(n)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor\) ，将其代入上式：
\[\sum_{i=1}^n\sum_{j=1}^md(ij)=\sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)f(\lfloor\frac{n}{d}\rfloor)f(\lfloor\frac{m}{d}\rfloor)
\]
到这里，很多题解就直接采用整除分块来计算 \(f\) 了，复杂度 \(O(n\sqrt{n})\) 。但是其实还有一种更快的做法能做到 \(O(n)\) 预处理。
不难发现 \(f(n)\) 其实是计算了 \(1 \sim n\) 每个数小于等于 \(n\) 的倍数个数之和，也就是说一个数对 \(f(n)\) 做出的贡献就是它的约数个数。那么有：
\[f(n)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor=\sum_{i=1}^n{\rm{d}}(i)
\]
由于 \({\rm{d}}\) 是积性函数，可以使用线性筛做到 \(O(n)\) 预处理。
关于线性筛 \({\rm{d}}\) ，可以看这个。
\({\rm{Code}}:\)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define maxn 50005
#define Rint register int
#define INF 0x3f3f3f3f
using namespace std;
typedef long long lxl;
template <typename T>
inline T read()
{
T x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
int prime[maxn],cnt;
bool flag[maxn];
lxl d[maxn],f[maxn],mu[maxn],sum[maxn],num[maxn];
inline void sieve()
{
d[1]=1;
mu[1]=1;
for(int i=2;i<maxn;++i)
{
if(!flag[i]) prime[++cnt]=i,d[i]=2,num[i]=1,mu[i]=-1;
for(int j=1;j<=cnt&&i*prime[j]<maxn;++j)
{
flag[i*prime[j]]=true;
if(i%prime[j])
{
num[i*prime[j]]=1;
d[i*prime[j]]=d[i]*2;
}
else
{
num[i*prime[j]]=num[i]+1;
d[i*prime[j]]=d[i]/num[i*prime[j]]*(num[i*prime[j]]+1);
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<maxn;++i)
f[i]=f[i-1]+d[i],sum[i]=sum[i-1]+mu[i];
}
inline lxl calcu(int n,int m)
{
lxl res=0;
if(n>m) swap(n,m);
for(int l=1,r=0;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
res+=f[n/l]*f[m/l]*(sum[r]-sum[l-1]);
}
return res;
}
int main()
{
// freopen("P3327.in","r",stdin);
sieve();
int T=read<int >();
while(T--)
{
int n=read<int >(),m=read<int >();
printf("%lld\n",calcu(n,m));
}
return 0;
}
