标签:
题意: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 __int64 const 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; }
版权声明:本文为博主原创文章,未经博主允许不得转载。
标签:
原文地址:http://blog.csdn.net/trq1995/article/details/47208195