【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加速。然而生成函数本身是用来算多重集组合数的，所以本题在使用的时候还需要进行对不合法方案的排除，这里用到了容斥原理（勉为其难地如此说吧）。

 

/**************************************************************
    Problem: 3771
    User: Tunix
    Language: C++
    Result: Accepted
    Time:1500 ms
    Memory:9464 kb
****************************************************************/
 
//BZOJ 3771
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define rep(i,n) for(int i=0;i<n;++i)
#define F(i,j,n) for(int i=j;i<=n;++i)
#define D(i,j,n) for(int i=j;i>=n;--i)
using namespace std;
int getint(){
    int v=0,sign=1; char ch=getchar();
    while(!isdigit(ch)) {if(ch=='-') sign=-1; ch=getchar();}
    while(isdigit(ch))  {v=v*10+ch-'0'; ch=getchar();}
    return v*sign;
}
/*******************template********************/
const int N=131072;
const double pi=acos(-1),eps=1e-3;
struct comp{
    double r,i;
    comp(double _=0.0,double __=0.0):r(_),i(__){}
    comp operator + (const comp&b)const{return comp(r+b.r,i+b.i);}
    comp operator - (const comp&b)const{return comp(r-b.r,i-b.i);}
    comp operator * (const comp&b)const{return comp(r*b.r-i*b.i,r*b.i+i*b.r);}
    friend comp operator * (double x,comp b){return comp(b.r*x,b.i*x);}
}a[N],b[N],c[N],ans[N];
void FFT(comp *a,int n,int type){
    for(int i=1,j=0;i<n-1;++i){
        for(int s=n;j^=s>>=1,~j&s;);
        if (i<j) swap(a[i],a[j]);
    }
    for(int m=1;m<n;m<<=1){
        double u=pi/m*type; comp wm(cos(u),sin(u));
        for(int i=0;i<n;i+=(m<<1)){
            comp w(1,0);
            rep(j,m){
                comp &A=a[i+j+m],&B=a[i+j],t=w*A;
                A=B-t; B=B+t; w=w*wm;
            }
        }
    }
    if (type==-1) rep(i,n) a[i].r/=n;
}
 
int main(){
//  freopen("input.txt","r",stdin);
    int n=getint(),len=N,x;
//  for(len=1;len<=(n<<1);len<<=1);
    rep(i,n){
        x=getint();
        a[x  ]=comp(1,0);
        b[x*2]=comp(1,0);
        c[x*3]=comp(1,0);
    }
    FFT(a,len,1); FFT(b,len,1); FFT(c,len,1);
    rep(i,len){
        ans[i]=ans[i]+(1.0/6.0)*(a[i]*a[i]*a[i]-3.0*b[i]*a[i]+2*c[i]);
        ans[i]=ans[i]+(0.5)*(a[i]*a[i]-b[i]);
        ans[i]=ans[i]+a[i];
    }
    FFT(ans,len,-1);
    rep(i,len){
        x=int(ans[i].r+eps);
        if (x) printf("%d %d\n",i,x);
    }
    return 0;
}