[bzoj3771]Triple——生成函数+容斥原理+FFT

题目大意：

我们讲一个悲伤的故事。
从前有一个贫穷的樵夫在河边砍柴。
这时候河里出现了一个水神，夺过了他的斧头，说：
“这把斧头，是不是你的？”
樵夫一看：“是啊是啊！”
水神把斧头扔在一边，又拿起一个东西问：
“这把斧头，是不是你的？”
樵夫看不清楚，但又怕真的是自己的斧头，只好又答：“是啊是啊！”
水神又把手上的东西扔在一边，拿起第三个东西问：
“这把斧头，是不是你的？”
樵夫还是看不清楚，但是他觉得再这样下去他就没法砍柴了。
于是他又一次答：“是啊是啊！真的是！”
水神看着他，哈哈大笑道：
“你看看你现在的样子，真是丑陋！”
之后就消失了。
樵夫觉得很坑爹，他今天不仅没有砍到柴，还丢了一把斧头给那个水神。
于是他准备回家换一把斧头。
回家之后他才发现真正坑爹的事情才刚开始。
水神拿着的的确是他的斧头。
但是不一定是他拿出去的那把，还有可能是水神不知道怎么偷偷从他家里拿走的。
换句话说，水神可能拿走了他的一把，两把或者三把斧头。
樵夫觉得今天真是倒霉透了，但不管怎么样日子还得过。
他想统计他的损失。
樵夫的每一把斧头都有一个价值，不同斧头的价值不同。总损失就是丢掉的斧头价值和。
他想对于每个可能的总损失，计算有几种可能的方案。
注意：如果水神拿走了两把斧头a和b，(a,b)和(b,a)视为一种方案。拿走三把斧头时，(a,b,c),(b,c,a),(c,a,b),(c,b,a),(b,a,c),(a,c,b)视为一种方案。

思路：

题面真好玩
考虑斧头价值的生成函数，自己和自己做卷积，然后发现二次多项式和三次多项式里面都有一些不合法的，直接容斥减掉即可。
代码实现上，可以先把多项式全部都变成点值表示，最后再逆变换回来即可。

/*=======================================
 * Author : ylsoi
 * Time : 2019.1.30
 * Problem : bzoj3771
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("bzoj3771.in","r",stdin);
    freopen("bzoj3771.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T fl=1; char ch=getchar();
	for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
	for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
	_*=fl;
}

const int maxn=240000+10;
const double pi=acos(-1);

struct cp{
	double x,y;
	cp(double xx=0,double yy=0){
		x=xx,y=yy;
	}
};
cp operator + (cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator - (cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator * (cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp operator / (cp a,double b){return cp(a.x/b,a.y/b);}
cp operator * (cp a,double b){return cp(a.x*b,a.y*b);}

cp g[maxn],ig[maxn];
int lim,cnt,dn[maxn];

void fft(cp *A,int ty){
	REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
	for(int len=1;len<lim;len<<=1){
		cp w= ty==1 ? g[len<<1] : ig[len<<1];
		for(int L=0;L<lim;L+=len<<1){
			cp wk=cp(1,0);
			REP(i,L,L+len-1){
				cp u=A[i],v=A[i+len]*wk;
				A[i]=u+v;
				A[i+len]=u-v;
				wk=wk*w;
			}
		}
	}
}

int n,ans[maxn];
cp a[maxn],b[maxn],c[maxn],y[maxn],z[maxn];

void init(){
	lim=1,cnt=0;
	while(lim<=12e4)lim<<=1,++cnt;
	if(!cnt)cnt=1;
	REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
	g[lim]=cp(cos(pi*2.0/lim),sin(pi*2.0/lim));
	ig[lim]=cp(cos(pi*2.0/lim),-sin(pi*2.0/lim));
	for(int i=lim>>1;i;i>>=1){
		g[i]=g[i<<1]*g[i<<1];
		ig[i]=ig[i<<1]*ig[i<<1];
	}
}

void work(){
	read(n);
	int t;
	REP(i,1,n)read(t),a[t].x+=1;
	REP(i,0,lim-1)if(a[i].x!=0){
		b[i*2].x+=a[i].x*a[i].x;
		c[i*3].x+=a[i].x*a[i].x*a[i].x;
	}

	fft(a,1),fft(b,1),fft(c,1);
	REP(i,0,lim-1){
		y[i]=(a[i]*a[i]-b[i])/2;
		z[i]=(a[i]*a[i]*a[i]-a[i]*b[i]*3+c[i]*2)/6;
	}
	fft(a,-1),fft(y,-1),fft(z,-1);
	REP(i,0,lim-1){
		a[i].x/=lim;
		y[i].x/=lim;
		z[i].x/=lim;
		ans[i]+=(int)(a[i].x+0.5)+(int)(y[i].x+0.5)+(int)(z[i].x+0.5);
	}
	REP(i,0,lim-1)if(ans[i])
		printf("%d %d
",i,ans[i]);
}

int main(){
	File();
	init();
	work();
	return 0;
}