[SDOI2015]约数个数和
给出T组询问，询问\(\sum_{i=1}^N\sum_{j=1}^Md(ij)\),1<=N, M<=50000，1<=T<=50000。
解
显然为约数计数问题，考虑Mobius反演，但此处不存在gcd，于是考虑拆分d(ij)为d(i),d(j),设\(N\leq M\),所以我们有结论\(d(ij)=\sum_{x|i}\sum_{y|j}\epsilon(gcd(x,y))\),代入
\[ans=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}\epsilon(gcd(x,y))=
\]
\[\sum_{x=1}^N\sum_{y=1}^M\epsilon(gcd(x,y))[N/x][M/y]
\]
于是设
\[f(k)=\sum_{x=1}^N\sum_{y=1}^M(gcd(x,y)==k)[N/x][M/y]
\]
\[F(k)=\sum_{x=1}^N\sum_{y=1}^M(k|gcd(x,y))[N/x][M/y]=
\]
\[\sum_{x=1}^{[N/k]}\sum_{y=1}^{[M/k]}[N/kx][M/ky]
\]
由Mobius反演定理，我们有
\[f(k)=\sum_{k|d}F(d)\mu(d/k)=\sum_{k|d}\mu(d/k)\sum_{x=1}^{[N/d]}\sum_{y=1}^{[M/d]}[N/dx][M/dy]
\]
\[ans=f(1)=\sum_{d=1}^N\mu(d)\sum_{x=1}^{[N/d]}[N/dx]\sum_{y=1}^{[M/d]}[M/dy]
\]
由我们的约数的性质，我们可以知道后式为约数前缀和的形式，我们可以递推出来，至于前半部分只要预处理出\(\mu\)的值，维护前缀和，在对后式利用整除分块即可，易知时间复杂度\(O(\sqrt{n})\)。
参考代码：
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int xdk[50001],prime[50001],pt,
mb[50001],cjx[50001];
il void read(int&);
il int min(int,int);
void pen(ll),prepare(int);
int main(){
int lsy,a,b,i,j;ll ans;
read(lsy),prepare(50000);
while(lsy--){
read(a),read(b),ans&=0;
if(a>b)swap(a,b);
for(i=1;i<=a;i=j+1)
j=min(a/(a/i),b/(b/i)),
ans+=(ll)xdk[a/i]*xdk[b/i]*(mb[j]-mb[i-1]);
pen(ans),putchar('\n');
}
return 0;
}
il int min(int a,int b){
return a<b?a:b;
}
void prepare(int n){
int i,j;
mb[1]|=xdk[1]|=check[1]|=cjx[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mb[i]=-1,xdk[i]=2,cjx[i]=1;
for(j=1;prime[j]*i<=n&&j<=pt;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j])){
cjx[i*prime[j]]=cjx[i]+1;
xdk[i*prime[j]]=xdk[i]*(cjx[i*prime[j]]+1)/(cjx[i]+1);
break;
}
mb[i*prime[j]]=-mb[i],cjx[i*prime[j]]|=true;
xdk[i*prime[j]]=xdk[i]<<1;
}
}for(i=1;i<=n;++i)mb[i]+=mb[i-1],xdk[i]+=xdk[i-1];
}
void pen(ll x){
if(x>9)pen(x/10);putchar(x%10+48);
}
il void read(int &x){
x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}
