给出T组询问,询问(sum_{i=1}^Nsum_{j=1}^Md(ij)),1<=N, M<=50000,1<=T<=50000。
解
显然为约数计数问题,考虑Mobius反演,但此处不存在gcd,于是考虑拆分d(ij)为d(i),d(j),设(Nleq M),所以我们有结论(d(ij)=sum_{x|i}sum_{y|j}epsilon(gcd(x,y))),代入
[ans=sum_{i=1}^Nsum_{j=1}^Msum_{x|i}sum_{y|j}epsilon(gcd(x,y))=
]
[sum_{x=1}^Nsum_{y=1}^Mepsilon(gcd(x,y))[N/x][M/y]
]
于是设
[f(k)=sum_{x=1}^Nsum_{y=1}^M(gcd(x,y)==k)[N/x][M/y]
]
[F(k)=sum_{x=1}^Nsum_{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}^Nmu(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('
');
}
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();
}