[BZOJ3771]Triple

Posted 2020-10-31 jefflyy

真是个愉悦的故事>_<

$\newcommand{\pretty}[1]{\begin{align*}#1\end{align*}}$设$\pretty{A(x)=\sum\limits_{i=1}^nx^{a_i},B(x)=\sum\limits_{i=1}^nx^{2a_i},C(x)=\sum\limits_{i=1}^nx^{3a_i}}$，那么只拿一个斧头，价值为$k$的的方案数是$\left[x^k\right]A(x)$，拿两个斧头，价值为$k$的方案数是$\left[x^k\right]\dfrac{A^2(x)-B(x)}2$（把“一把斧头拿两次”的情况容斥掉再消序），拿三把斧头，价值为$k$的方案数是$\left[x^k\right]\dfrac{A^3(x)-3A(x)B(x)+2C(x)}6$（容斥掉“一个斧头选两次，另一个斧头选一次”$\text{abb,bab,bba!}$和“一个斧头选三次”的情况，最后消序）

#include<stdio.h>
#include<math.h>
typedef double du;
typedef long long ll;
struct complex{
	du x,y;
	complex(du a=0,du b=0){x=a;y=b;}
};
complex operator+(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator-(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator*(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int rev[262144],N,iN;
void pre(int n){
	int i,j,k;
	for(N=1,k=0;N<n;N<<=1)k++;
	for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
void swap(complex&a,complex&b){
	complex c=a;
	a=b;
	b=c;
}
void fft(complex*a,int on){
	int i,j,k;
	complex t;
	for(i=0;i<N;i++){
		if(i<rev[i])swap(a[i],a[rev[i]]);
	}
	for(i=2;i<=N;i<<=1){
		for(j=0;j<N;j+=i){
			for(k=0;k<i>>1;k++){
				t=complex(cos(k*M_PI/(i/2)),on*sin(k*M_PI/(i/2)))*a[i/2+j+k];
				a[i/2+j+k]=a[j+k]-t;
				a[j+k]=a[j+k]+t;
			}
		}
	}
	if(on==-1){
		for(i=0;i<N;i++)a[i].x/=N;
	}
}
int max(int a,int b){return a>b?a:b;}
complex a[262144],b[262144],c[262144];
int main(){
	int n,m,i,x;
	ll s;
	scanf("%d",&n);
	m=0;
	for(i=0;i<n;i++){
		scanf("%d",&x);
		a[x].x+=1;
		b[x*2].x+=1;
		c[x*3].x+=1	;
		m=max(m,x*3);
	}
	pre(m<<1|1);
	fft(a,1);
	fft(b,1);
	fft(c,1);
	for(i=0;i<N;i++)a[i]=(a[i]*a[i]*a[i]-3*a[i]*b[i]+2*c[i])*(1/6.)+(a[i]*a[i]-b[i])*.5+a[i];
	fft(a,-1);
	for(i=1;i<=m;i++){
		s=llround(a[i].x);
		if(s)printf("%d %lld\n",i,s);
	}
}