首页 > 代码库 > hdu4746 Mophues
hdu4746 Mophues
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=4746
题意:给出n, m, p,求有多少对a, b满足gcd(a, b)的素因子个数<=p,(其中1<=a<=n, 1<=b<=m)
分析:
设A(d):gcd(a, b)=d的有多少种
设B(j): gcd(a, b)是j的倍数的有多少种,易知B(j) = (n/j)*(m/j)
则由容斥原理得:(注:不同行的μ是不相同的,μ为莫比乌斯函数)
A(1) = μ(1)*B(1) + μ(2)*B(2) + μ(3)*B(3) + ... + μ(p1*p2...)*B(p1*p2...)
A(2) = μ(1)*B(1*2) + μ(2)*B(2*2) + μ(3)*B(3*2) + ... + μ(p1*p2..)*B(p1*p2..*2)
...
A(d) = μ(1)*B(1*d) + μ(2)*B(2*d) + μ(3)*B(3*d) + ... + μ(p1*p2..)*B(p1*p2..*d)
ans = A(1)+A(2)+...+A(d) = F(1)*B(1) + F(2)*B(2) + ... + F(p1*p2..)*B(p1*p2..)
于是可以枚举公约数i{表示A(i)},利用筛法找出i的倍数j,i对B(j)的贡献系数为:F(j)+=μ(j/i)
总之,求出B(j)的总贡献系数F(j)即可得答案:F(1)*B(1)+F(2)*B(2)+...+F(n)*B(n)
上面没有限制gcd的素因子个数,要限制其实不难,给系数加多一维即可:
F(d)(p)表示:素因子个数<=p时,对B(d)的贡献系数
代码:
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #include<cstring> 5 6 using namespace std; 7 const int maxn=500010; 8 9 int n,m,p; 10 int vis[maxn]; 11 int prime[maxn]; 12 int cnt; 13 int num[maxn]; 14 int mu[maxn]; 15 int f[maxn][20]; 16 17 void init() 18 { 19 memset(vis,0,sizeof(vis)); 20 cnt=0; 21 mu[1]=1; 22 for(int i=2;i<maxn;i++) 23 { 24 if(!vis[i]) 25 { 26 prime[cnt++]=i; 27 mu[i]=-1; 28 num[i]=1; 29 } 30 for(int j=0;j<cnt&&i*prime[j]<maxn;j++) 31 { 32 vis[i*prime[j]]=1; 33 num[i*prime[j]]=num[i]+1; 34 if(i%prime[j]) 35 mu[i*prime[j]]=-mu[i]; 36 else 37 { 38 mu[i*prime[j]]=0; 39 break; 40 } 41 } 42 } 43 memset(f,0,sizeof(f)); 44 for(int i=1;i<maxn;i++) 45 for(int j=i;j<maxn;j+=i) 46 f[j][num[i]]+=mu[j/i]; 47 48 for(int i=0;i<maxn;i++) 49 for(int j=1;j<20;j++) 50 f[i][j]+=f[i][j-1]; 51 52 for(int i=1;i<maxn;i++) 53 for(int j=0;j<20;j++) 54 f[i][j]+=f[i-1][j]; 55 } 56 57 int main() 58 { 59 init(); 60 int T; 61 scanf("%d",&T); 62 while(T--) 63 { 64 scanf("%d%d%d",&n,&m,&p); 65 if(p>=20) 66 { 67 printf("%lld\n",(long long)n*m); 68 continue; 69 } 70 if(n>m) 71 swap(n,m); 72 int last; 73 long long ans=0; 74 for(int i=1;i<=n;i=last+1) 75 { 76 last=min(n/(n/i),m/(m/i)); 77 ans+=(long long)(f[last][p]-f[i-1][p])*(n/i)*(m/i); 78 } 79 printf("%lld\n",ans); 80 } 81 return 0; 82 }
hdu4746 Mophues