pku3904 容斥原理的运用,好题!

解法发现网上一个大牛写的很好,就转载过来吧。可怜的我,开始没想到算法就算了,想到算法后又刻意依赖STL结果o(N)写成o(logN)成功葬送。
解法

题意:给你n个数,求GCD(gcd(a,b),gcd(c,d))=1的方案数,注意a,b,c,d并不一定两两互质,有多组数据

本来这个题的题解有一大把,但没有讲题解的只有贴代码的。

本来一直在做题,没有什么时间发题解,但这个题加深了我对容斥原理的理解,所以讲一下

这个题的算法是个伪多项式

首先,先将算法流程说一下,原理后面会有

Step 1:用筛法筛出10000内的素数表

Step 2:组合素数,用bool数组记录由m(m<=5)个素数相乘所得数的m的奇偶

例如:182=2*7*13 ———> m=3 so,bool[182]=false

            2=2 ——>m=1 so,bool[2]=false

            91=7*13 ——>m=2 so,bool[91]=true

            注意12=2^2*3等这种是B=p1^k1*p2^k2+...这种有质因数指数不为一的合数不要处理因为不会用到

以上是预处理

Step 3:读入当前数为a,将a分解为质因数形式,注意每个质因数只被记录一次

例如:12=2*2*3 则 12会被分为2*3,多余的2被消去了

Step 4:将分出的素数进行组合,将组合出的a的约数的time+1

例如:12=2^2*3——>12=2*3——>time[2]++,time[3]++,time[6]++

Step 5:处理之后,ans赋值为C(n,4)即n*(n-1)*(n-2)*(n-3)div 24

Step 6:for i 2-->10000

                if bool[i] then ans:=ans-C(time[i],4)

                               else ans:=and+C(time[i],4);

Step 7:输出ans


原理:首先考虑一个问题,1000以内6,7,8,9的倍数有多少个?答案是

1000div6+1000div7+1000div8+1000div9

-1000div(6*7)-1000div(6*8)-1000div(6*9)-1000div(7*8)-1000div(7*9)-1000div(8*9)

+1000div(6*7*8)+1000div(6*8*9)+1000div(7*8*9)

-1000div(6*7*8*9)

这是容斥原理的一个最简单的应用,类比这道题,Step3到4其实将每个数a的不重复约数记录了下来,有公共约数的四个数的方案要从ans中减去,多减的要加上

显然m为奇时要减,m为偶时要加,这和”1000以内6,7,8,9的倍数有多少个?“这个问题奇偶是反的,由于a最大为10000,所以m最大可以有5 (2*3*5*7*11<10000,2*3*5*7*11*13>10000)

至于12=2*2*3这种约数不处理因为可以分为2*6,而2和6会算一次,所以不须再算。


代码:
 1 # include <cstdio>
 2 # include <map>
 3 # include <cstring>
 4 # include <algorithm>
 5 using namespace std;
 6 int pa[2000],pp=0,sa[10],sp=0;
 7 int refer[5][10001];
 8 void make_prime()
 9 {
10     bool used[10001];
11     memset(used,true,sizeof(used));
12     for(int i=2;i<=10000;i++)
13          if(used[i])
14          {
15              pa[pp++]=i;
16              for(int j=2*i;j<=10000;j+=i)
17                 used[j]=false;
18          }
19 }
20 void spilt(int n)
21 {
22     sp=0;
23     for(int i=0;i<pp&&n!=1;i++)
24         if(n%pa[i]==0)
25         {
26             sa[sp++]=pa[i];
27             while(n%pa[i]==0)
28                 n/=pa[i];
29         }
30 }
31 void putmap(int n)
32 {
33     spilt(n);
34     for(int i=1;i<(1<<sp);i++)
35     {
36         int n1=0,n2=1;
37         for(int j=0;j<5;j++)
38             if(i&(1<<j))
39                 n1++,n2*=sa[j];
40         refer[n1-1][n2]++;
41     }
42 }
43 long long c(int n)
44 {
45     return (long long)n*(n-1)*(n-2)*(n-3)/4/3/2;
46 }
47 long long getans(int n)
48 {
49     long long ans=c(n);
50     for(int i=0;i<5;i++)
51     {
52         bool flag=false;
53         for(int j=1;j<=10000;j++)
54             if(refer[i][j]>=4)
55                 flag=true,
56                 ans+=c(refer[i][j])*(i%2?1ll:-1ll);
57         if(!flag)break;
58     }
59     return ans;
60 }
61 int main()
62 {
63     int n;
64     make_prime();
65     while(scanf("%d",&n)!=EOF)
66     {
67         memset(refer,0,sizeof(refer));
68         for(int i=0;i<n;i++)
69         {
70             int t;
71             scanf("%d",&t);
72             putmap(t);
73         }
74         printf("%lld\n",getans(n));
75     }
76     return 0;
77 }

posted on 2012-02-17 02:03 yzhw 阅读(263) 评论(0)  编辑 收藏 引用 所属分类: combination math


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理


<2012年2月>
2930311234
567891011
12131415161718
19202122232425
26272829123
45678910

导航

统计

公告

统计系统

留言簿(1)

随笔分类(227)

文章分类(2)

OJ

最新随笔

搜索

积分与排名

最新评论

阅读排行榜