BZOJ.3771.Triple(母函数 FFT 容斥)

题目链接

\(Description\)

　　有\(n\)个物品(斧头)，每个物品价值不同且只有一件，问取出一件、两件、三件物品，所有可能得到的价值和及其方案数。\((a,b),(b,a)\)算作一种方案。

\(Solution\)

　　尝试用母函数去表示。\(A\)表示取一个物品对应方案数和价值的母函数，即$$A=x^{v_1}+x+x^{v_3}\ldots$$
　　那么取两件就是\(A^2\)，取三件是\(A^3\)。因为是组合，而这么求的是排列，所以$$Ans=A+\frac{A^{2}{2!}+\frac{A}3}{3!}$$
　　但是有问题，\(A\times A\)会出现选了同一件物品的情况，所以要减掉。考虑构造选了同一个物品两次的母函数\(B\)，即\(B=x^{2v_1}+x^{2v_2}+x^{2v_3}\ldots\)
　　那么取两件对应的\(Ans'=\frac{A^2-B}{2!}\)
　　同理，构造取了同样物品三件的母函数\(C=x^{3v_1}+x^{3v_2}+x^{3v_3}+\ldots\)
　　对于三件物品选了两件同样物品，\(A\times A\times A\ ->\ A\times B\) 可能会出现\((x,x,y),(x,y,x),(y,x,x)\)三种情况，所以减去\(3AB\)。而\(A^3\)和\(AB\)中都包含选三件同样物品的方案，多减了两次，所以要加\(2C\)。对于取三件对应的母函数\(Ans''=\frac{A^3-3AB+2C}{3!}\)
　　所以

\[Ans=A+\frac{A^2-B}{2}+\frac{A^3-3AB+2C}{6} \]

　　多项式的点值表示下很多可以直接运算啊。。（并不晓得。。）

//7476kb	852ms
#include <cmath>
#include <cstdio>
#include <cctype>
#include <algorithm>
#define gc() getchar()
const int N=(1<<17)+5;//131072
const double PI=acos(-1);

int n,rev[N];
struct Complex
{
	double x,y;
	Complex() {}
	Complex(double x,double y):x(x),y(y) {}
	Complex operator +(const Complex &a) {return Complex(x+a.x, y+a.y);}//有点纠结这换不换行。。
	Complex operator -(const Complex &a) {return Complex(x-a.x, y-a.y);}
	Complex operator *(const Complex &a) {return Complex(x*a.x-y*a.y, x*a.y+y*a.x);}
	Complex operator *(const double &a) {return Complex(x*a, y*a);}
	Complex operator /(const double &a) {return Complex(x/a, y/a);}
}A[N],B[N],C[N];

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
void FFT(Complex *a,int lim,int type)
{
	for(int i=1; i<lim; ++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
	for(int i=2; i<=lim; i<<=1)
	{
		int mid=i>>1;
		Complex Wn(cos(PI/mid),type*sin(PI/mid)),t;//2.0*PI/i
//		W[0]=Complex(1,0);
//		for(int j=1; j<mid; ++j) W[j]=W[j-1]*Wn;//这两行出俩错误→_→ 
		for(int j=0; j<lim; j+=i)
		{
			Complex w(1,0);//不预处理更快。。是范围小的原因吗。。
			for(int k=0; k<mid; ++k,w=w*Wn)
				a[j+k+mid]=a[j+k]-(t=a[j+k+mid]*w),
				a[j+k]=a[j+k]+t;
		}
	}
	if(type==-1)
		for(int i=0; i<lim; ++i) a[i].x/=lim;
}

int main()
{
	n=read(); int mx=0;
	for(int v,i=1; i<=n; ++i)
		mx=std::max(mx,v=read()), A[v].x=B[2*v].x=C[3*v].x=1.0;
	int lim=1, L=0;
	while(lim<=3*mx) lim<<=1, ++L;
	for(int i=1; i<lim; ++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<L-1);

	FFT(A,lim,1), FFT(B,lim,1), FFT(C,lim,1);
	for(int i=0; i<lim; ++i)
		A[i]=A[i]+(A[i]*A[i]-B[i])*0.5+(A[i]*A[i]*A[i]-A[i]*B[i]*3.0+C[i]*2.0)/6.0;
	FFT(A,lim,-1);
	for(int i=0; i<lim; ++i)
		if((int)(A[i].x+0.5)) printf("%d %d\n",i,(int)(A[i].x+0.5));//这题还不需要longlong

	return 0;
}