生成函数+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=,sign=; char ch=getchar();
while(!isdigit(ch)) {if(ch=='-') sign=-; ch=getchar();}
while(isdigit(ch)) {v=v*+ch-''; ch=getchar();}
return v*sign;
}
/*******************template********************/
const int N=;
const double pi=acos(-),eps=1e-;
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=,j=;i<n-;++i){
for(int s=n;j^=s>>=,~j&s;);
if (i<j) swap(a[i],a[j]);
}
for(int m=;m<n;m<<=){
double u=pi/m*type; comp wm(cos(u),sin(u));
for(int i=;i<n;i+=(m<<)){
comp w(,);
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==-) 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(,);
b[x*]=comp(,);
c[x*]=comp(,);
}
FFT(a,len,); FFT(b,len,); FFT(c,len,);
rep(i,len){
ans[i]=ans[i]+(1.0/6.0)*(a[i]*a[i]*a[i]-3.0*b[i]*a[i]+*c[i]);
ans[i]=ans[i]+(0.5)*(a[i]*a[i]-b[i]);
ans[i]=ans[i]+a[i];
}
FFT(ans,len,-);
rep(i,len){
x=int(ans[i].r+eps);
if (x) printf("%d %d\n",i,x);
}
return ;
}