【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加速。然而生成函数本身是用来算多重集组合数的,所以本题在使用的时候还需要进行对不合法方案的排除,这里用到了容斥原理(勉为其难地如此说吧)。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
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 }