bzoj3771: Triple(容斥+生成函数+FFT)

传送门

咳咳忘了容斥了……

\(A(x)\)为斧头的生成函数,其中第\(x^i\)项的系数为价值为\(i\)的斧头个数,那么\(A(x)+A^2(x)+A^3(x)\)就是答案(于是信心满满的打了一发连样例都没过)

如果按上面那样算的话,会有重复的,比如说\(A^2(x)\),会产生诸如\((x_i,x_i)\)之类的同一把斧头的贡献,所以定义\(B(x)\)为同一个斧头重复两次的方案数,那么\(A^2(x)-B(x)\)就是两把斧头时真正的贡献,又因为与顺序无关,所以还要除以\(2\)

然后\(A^3(x)\)的话,可能会有一把斧头重复两次或三次,如果重复两次,那么就是\((x_i,x_i,y_i),(x_i,y_i,x_i),(y_i,x_i,x_i)\),就是\(3A(x)B(x)\),但是减去这个的话又会把\((x_i,x_i,x_i)\)的情况多减去两次,所以定义\(C(x)\)为同一把斧头重复三次的生成函数,于是还要加上\(2C(x)\),然后无关顺序的话还要除掉\(3!=6\)

综上,最终的答案的生成函数为$$Ans(x)=A(x)+\frac{A2(x)-B(x)}{2}+\frac{A3(x)-3A(x)B(x)+2C(x)}{6}$$

//minamoto
#include<cstdio>
#include<cmath>
#include<algorithm>
#define R register
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
const int N=6e5+5;const double Pi=acos(-1.0);
struct complex{
	double x,y;
	complex(double xx=0,double yy=0){x=xx,y=yy;}
	inline complex operator +(const complex &b)const{return complex(x+b.x,y+b.y);}
	inline complex operator -(const complex &b)const{return complex(x-b.x,y-b.y);}
	inline complex operator *(const complex &b)const{return complex(x*b.x-y*b.y,x*b.y+y*b.x);}
	inline complex operator *(const int &b){return complex(x*b,y*b);}
	inline complex operator /(const int &b){return complex(x/b,y/b);}
}A[N],B[N],C[N],O[N],ans[N];
int r[N],lim,n,x,l,m;
void FFT(complex *A,int ty){
	fp(i,0,lim-1)if(i<r[i])swap(A[i],A[r[i]]);
	for(R int mid=1;mid<lim;mid<<=1){
		int I=(mid)<<1;
        complex Wn(cos(Pi/mid),ty*sin(Pi/mid));
		fp(i,1,mid-1)O[i]=O[i-1]*Wn;
		for(R int j=0;j<lim;j+=I)fp(k,0,mid-1){
			complex x=A[j+k],y=O[k]*A[j+k+mid];
			A[j+k]=x+y,A[j+k+mid]=x-y;
		}
	}if(ty==-1)fp(i,0,lim-1)A[i].x=(int)(A[i].x/lim+0.5);
}
int main(){
//	freopen("testdata.in","r",stdin);
	n=read();
	fp(i,1,n)x=read(),++A[x].x,++B[x<<1].x,++C[(x<<1)+x].x,cmax(m,x);
	m*=3,lim=1;while(lim<=m)lim<<=1,++l;O[0]=complex(1,0);
	fp(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	FFT(A,1),FFT(B,1),FFT(C,1);
	fp(i,0,lim-1)ans[i]=A[i]+(A[i]*A[i]-B[i])/2+(A[i]*A[i]*A[i]-A[i]*B[i]*3+C[i]*2)/6;
	FFT(ans,-1);
	fp(i,0,m)if(ans[i].x)printf("%d %.0lf\n",i,ans[i].x);
	return 0;
}
深深地明白自己的弱小
posted @ 2019-01-02 21:40  bztMinamoto  阅读(151)  评论(0编辑  收藏  举报
Live2D