HDU 4746 Mophues 莫比乌斯第三弹

分类：数论

题意：1<=x<=n,1<=y<=m,使得gcd(x,y)=k，k的素因数个数小于等于p

例：24=2*2*2*3，k=4

解：设f[n]为gcd(a,b)=n的对数

F[d]为d|gcd(a,b)的对数

f[n]=sigema(mu[i],F[i*n])：

f[1]=mu[1]*F[1]+mu[2]*F[1*2]+…+mu[n]*F[1*n]

f[2]=mu[2]*F[2]+mu[2]*F[2*2]+…+mu[n]*F[2*n]

……

sum=f[1]+f[2]+…+f[n]=G[1]*F[1]+G[2]*F[2]+…+G[n]*F[n]

枚举每一个i，则i的倍数j为G[j]提供了mu[j/i]的贡献，即G[j]+=mu[j/i]

因为所求为k的素因数个数，所以将G[j]开成二维数组G[j][p]表示j对素因数个数为p的贡献

需要使用分块加速的方法，否则还是要爆炸

#include <stdio.h>#include <string.h>#define MIN(a,b) ((a)<(b)?(a):(b))#define MAX(a,b) ((a)<(b)?(b):(a))#define ll __int64const int maxn=500005;int num[maxn];int prime[maxn];int mu[maxn];int factor[maxn];int mbs[maxn][20];void mobius(){memset(num,0,sizeof(num));int all=0;mu[1]=1;factor[1]=0;for(int i=2;i<maxn;i++){if(!num[i]){prime[all++]=i;mu[i]=-1;factor[i]=1;//记录素因数个数}for(int j=0;j<all&&i*prime[j]<maxn;j++){num[i*prime[j]]=1;factor[i*prime[j]]=factor[i]+1;if(i%prime[j]){mu[i*prime[j]]=-mu[i];}else{mu[i*prime[j]]=0;break;}}}return ;}void inti(){memset(mbs,0,sizeof(mbs));for(int i=1;i<maxn;i++)for(int j=i;j<maxn;j+=i)mbs[j][factor[i]]+=mu[j/i];//每个j在factor[i]个素因数中的贡献/*下面是为了分块加速求和*/for(int i=1;i<maxn;i++)for(int j=0;j<19;j++)mbs[i][j]+=mbs[i-1][j]; for(int i=0;i<maxn;i++)for(int j=1;j<19;j++)mbs[i][j]+=mbs[i][j-1];return ;}int main(){int t,n,m,p;mobius();inti();ll sum;while(scanf("%d",&t)!=-1){while(t–){scanf("%d%d%d",&n,&m,&p);if(p>=19)//因为2^19>500000，所以超过了19就是全体均满足要求{printf("%I64d\n",(ll)n*m);continue;}if(n>m){int te=n;n=m;m=te;}sum=0;for(int i=1,last;i<n;i=last+1){last=MIN(n/(n/i),m/(m/i)); //分块加速，因为[n/i][m/i]在i递加过程中具有重复部分，跳掉这些i可是简化计算，，是复杂度降低至sqrt(n)sum+=((ll)(n/i)*(m/i)*(mbs[last][p]-mbs[i-1][p]));}printf("%I64d\n",sum);}}return 0;}

上一篇HDU 1695 GCD 莫比乌斯第二发下一篇HDU 4675 GCD of Sequence

顶0踩0

变幻原是永恒，我们唯有用永恒的诺言制约世事的变幻。

相关文章：

你感兴趣的文章：

标签云：