【BZOJ】【3771】Triple

生成函数+FFT

  Orz PoPoQQQ

  这个题要算组合的方案,而且范围特别大……所以我们可以利用生成函数来算

  生成函数是一个形式幂级数,普通生成函数可以拿来算多重集组合……好吧我承认以上是在瞎扯→_→

  这个东西我也不记得是多会儿看的了……找本《组合数学》自己看看好了……或者问学数学竞赛的同学&数竞教练

 

先引用下PoPoQQQ的题解:

首先搞出这n个物品的母函数a
将a的每项的平方求和得到多项式b
将a的每项的立方求和得到多项式c
那么如果不考虑顺序和重复 那么方案数就是a+b+c
现在考虑顺序和重复后
三个物品的方案数为(a^3-3*a*b+2*c)/6
两个物品的方案数为(a^2-b)/2
一个物品的方案数为a
故最终答案为(a^3-3*a*b+2*c)/6+(a^2-b)/2+a
用FFT搞一下就好了- -

  

  我们可以这样定义一个多项式(生成函数)A=a0+a1*x+a2*(x^2)+a3*(x^3)+……(为什么要说它是形式幂级数呢?我觉得是因为我们不考虑x的值具体是多少,只考虑x^i 的系数!运算的时候就像大整数乘法那样算)

  如果我们有一把斧头的价值是w,我们就令x^w这一项的系数为1,然后我们可以这样算:

    1.A的每一项 x^i 的系数表示取1柄斧头价值为 i 的方案数

    2.A*A的每一项 x^i 的系数表示取两柄斧头总价值为 i 的方案数:

    比如我们有价值为3、4、5和6的斧头,那么在A这个多项式中,x^3、x^4、x^5和x^6的系数就为1,其他都为0,所以在A*A里面,x^9的系数为2,因为(x^3)*(x^6)=x^9,(x^4)*(x^5)=x^9;所以x^9这一项的系数为2,也就是说我们取两柄斧头有两种方案使得总价值为9.

    但是这样算出来的并不全对,因为这样会出现 (x^3)*(x^3)=(x^6)的情况,也就是一把斧头取了两遍的方案,这种是需要减去的,所以我们再定义一个多项式B,满足 x^(w[i]*2) 的系数为1,那么A*A-B的答案就是总的方案数。

    但是!!这样的方案是算了2遍的,因为(x^3)*(x^6)=(x^9)且(x^6)*(x^3)=(x^9)(注意我们是让A自己乘自己)所以答案应该是 0.5 * (A*A-B)

    3.A*A*A表示取三把斧头,同样,我们要去掉不合法的方案,并且乘以(1/6),因为三个斧头我们一共算了6次(全排列)

 

【重要】:本题使用FFT的时候,是对w[i]的范围进行多项式乘法的,而不是斧头的数量N!!!这点要注意,所以所有的FFT的范围都需是131072!!

总结:本题是利用生成函数来进行组合方案数的求解,这种方法的优势在于是用多项式的系数来表示的,所以可以用FFT加速。然而生成函数本身是用来算多重集组合数的,所以本题在使用的时候还需要进行对不合法方案的排除,这里用到了容斥原理(勉为其难地如此说吧)。

 

 1 /**************************************************************
 2     Problem: 3771
 3     User: Tunix
 4     Language: C++
 5     Result: Accepted
 6     Time:1500 ms
 7     Memory:9464 kb
 8 ****************************************************************/
 9  
10 //BZOJ 3771
11 #include<cmath>
12 #include<cstdio>
13 #include<cstring>
14 #include<cstdlib>
15 #include<iostream>
16 #include<algorithm>
17 #define rep(i,n) for(int i=0;i<n;++i)
18 #define F(i,j,n) for(int i=j;i<=n;++i)
19 #define D(i,j,n) for(int i=j;i>=n;--i)
20 using namespace std;
21 int getint(){
22     int v=0,sign=1; char ch=getchar();
23     while(!isdigit(ch)) {if(ch=='-') sign=-1; ch=getchar();}
24     while(isdigit(ch))  {v=v*10+ch-'0'; ch=getchar();}
25     return v*sign;
26 }
27 /*******************template********************/
28 const int N=131072;
29 const double pi=acos(-1),eps=1e-3;
30 struct comp{
31     double r,i;
32     comp(double _=0.0,double __=0.0):r(_),i(__){}
33     comp operator + (const comp&b)const{return comp(r+b.r,i+b.i);}
34     comp operator - (const comp&b)const{return comp(r-b.r,i-b.i);}
35     comp operator * (const comp&b)const{return comp(r*b.r-i*b.i,r*b.i+i*b.r);}
36     friend comp operator * (double x,comp b){return comp(b.r*x,b.i*x);}
37 }a[N],b[N],c[N],ans[N];
38 void FFT(comp *a,int n,int type){
39     for(int i=1,j=0;i<n-1;++i){
40         for(int s=n;j^=s>>=1,~j&s;);
41         if (i<j) swap(a[i],a[j]);
42     }
43     for(int m=1;m<n;m<<=1){
44         double u=pi/m*type; comp wm(cos(u),sin(u));
45         for(int i=0;i<n;i+=(m<<1)){
46             comp w(1,0);
47             rep(j,m){
48                 comp &A=a[i+j+m],&B=a[i+j],t=w*A;
49                 A=B-t; B=B+t; w=w*wm;
50             }
51         }
52     }
53     if (type==-1) rep(i,n) a[i].r/=n;
54 }
55  
56 int main(){
57 //  freopen("input.txt","r",stdin);
58     int n=getint(),len=N,x;
59 //  for(len=1;len<=(n<<1);len<<=1);
60     rep(i,n){
61         x=getint();
62         a[x  ]=comp(1,0);
63         b[x*2]=comp(1,0);
64         c[x*3]=comp(1,0);
65     }
66     FFT(a,len,1); FFT(b,len,1); FFT(c,len,1);
67     rep(i,len){
68         ans[i]=ans[i]+(1.0/6.0)*(a[i]*a[i]*a[i]-3.0*b[i]*a[i]+2*c[i]);
69         ans[i]=ans[i]+(0.5)*(a[i]*a[i]-b[i]);
70         ans[i]=ans[i]+a[i];
71     }
72     FFT(ans,len,-1);
73     rep(i,len){
74         x=int(ans[i].r+eps);
75         if (x) printf("%d %d\n",i,x);
76     }
77     return 0;
78 }
View Code

 

posted @ 2015-02-06 19:08  Tunix  阅读(464)  评论(0编辑  收藏  举报