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
",(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
",ans);
80     }
81     return 0;
82 }