题目描述
珈百璃刚刚学习了莫比乌斯反演,想再做一道简单题练练手。
设$sigma _{0}left ( x
ight )$为$x$ 的约数个数,给出$N$ 和$M$,求:
$sum_{i=1}^{N}sum_{j=1}^{M}sigma _{0}left ( ij
ight )$
对于100%的数据,$Tleq 50000,1leq N,M leq50000$
题解
首先有一个性质$sigma _{0}left(ij
ight)= sum_{x|i}sum_{y|j}left [ gcdleft(x,y
ight )=1
ight ]$
那么答案就是$sum_{i=1}^{N}sum_{j=1}^{M} sum_{x|i}sum_{y|j}left [ gcdleft(x,y
ight )=1
ight ]$
先枚举$x,y$,$sum_{i=1}^{N}sum_{j=1}^{M}left [ gcdleft(i,j
ight )=1
ight ]left lfloor frac{N}{i}
ight
floorleft lfloor frac{M}{j}
ight
floor$
再进行反演,$sum_{d=1}^{min(n,m)}mu (d)sum_{i=1}^{frac{N}{d}}left lfloor frac{N}{di}
ight
floorsum_{j=1}^{frac{M}{d}}left lfloor frac{M}{dj}
ight
floor$
记$f(i)=sum_{j=1}^{i}left lfloor frac{i}{j}
ight
floor$
我们要求的就是$sum_{i=1}^{min(n,m)}mu(i)*f(left lfloor frac{N}{i}
ight
floor)*f(left lfloor frac{M}{i}
ight
floor)$
预处理$f$数组就ok了,其实$f(n)$数组也就是$[1,n]$的约数和,因为可以看做是枚举约数$j$.
#include<bits/stdc++.h> using namespace std; #define ll long long const int maxn=50006; int t,n,m; int cnt,prime[maxn],mu[maxn]; ll f[maxn],ans[6006][6006]; bool not_prime[maxn]; void init(){ mu[1]=1; for(int i=2;i<=50000;i++){ if(!not_prime[i]){ prime[++cnt]=i; mu[i]=-1; } for(int j=1;prime[j]*i<=50000;j++){ not_prime[i*prime[j]]=true; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else { mu[i*prime[j]]=0; break; } } mu[i]+=mu[i-1]; } for(int i=1;i<=50000;++i) for(ll l=1,r;l<=i;l=r+1){ r=i/(i/l); f[i]+=(r-l+1)*(i/l); } } int main(){ freopen("c.in","r",stdin); freopen("c.out","w",stdout); init(); scanf("%d",&t); while(t--){ ll ret=0; scanf("%d%d",&n,&m); if(n>m) swap(n,m); if(n<=6000&&m<=6000&&ans[n][m]){ printf("%lld ",ans[n][m]); continue; } for(int l=1,r;l<=n;l=r+1){ r=min(n/(n/l),m/(m/l)); ret+=(mu[r]-mu[l-1])*f[n/l]*f[m/l]; } if(n<=6000&&m<=6000) ans[n][m]=ans[m][n]=ret; printf("%lld ",ret); } }