2018.12.31 bzoj3771: Triple(生成函数+fft+容斥原理)

本文介绍了一种使用生成函数解决组合计数问题的方法,并通过快速傅立叶变换(FFT)优化多项式乘法,提高算法效率。具体展示了如何计算从n个数中选择1/2/3个数的所有可能和及其方案数。

传送门
生成函数经典题。
题意简述:给出nnn个数,可以从中选1/2/31/2/31/2/3个,问所有可能的和对应的方案数。


思路:
A(x),B(x),C(x)A(x),B(x),C(x)A(x),B(x),C(x)表示选111个,222个,333个的生成函数,ans1(x),ans2(x),ans3(x)ans1(x),ans2(x),ans3(x)ans1(x),ans2(x),ans3(x)表示选111个,222个,333个答案的生成函数。
那么ans1(x)=A(x)ans1(x)=A(x)ans1(x)=A(x)
ans2(x)=A2(x)−B(x)ans2(x)=A^2(x)-B(x)ans2(x)=A2(x)B(x)
ans3(x)=A3(x)−3A(x)B(x)+2C(x)ans3(x)=A^3(x)-3A(x)B(x)+2C(x)ans3(x)=A3(x)3A(x)B(x)+2C(x)
然后就可以算出答案了。
注意用fftfftfft优化多项式乘法。
代码:

#include<bits/stdc++.h>
#define ri register int
using namespace std;
const double pi=acos(-1.0);
const int N=5e5+5;
typedef long long ll;
struct Cp{
	double x,y;
	friend inline Cp operator+(const Cp&a,const Cp&b){return (Cp){a.x+b.x,a.y+b.y};}
	friend inline Cp operator-(const Cp&a,const Cp&b){return (Cp){a.x-b.x,a.y-b.y};}
	friend inline Cp operator*(const Cp&a,const Cp&b){return (Cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
	friend inline Cp operator*(const double&a,const Cp&b){return (Cp){a*b.x,a*b.y};}
	friend inline Cp operator/(const Cp&a,const double&b){return (Cp){a.x/b,a.y/b};}
}a[N],b[N],c[N];
int pos[N],lim=1,tim=0,mx=0,n;
inline void init(){
	while(lim<=mx*2)lim<<=1,++tim;
	for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}
inline int read(){
	int ans=0;
	char ch=getchar();
	while(!isdigit(ch))ch=getchar();
	while(isdigit(ch))ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar();
	return ans;
}
inline void fft(Cp a[],int type){
	for(ri i=0;i<lim;++i)if(i<pos[i])swap(a[i],a[pos[i]]);
	for(ri mid=1;mid<lim;mid<<=1){
		Cp wn=(Cp){cos(pi/mid),type*sin(pi/mid)};
		for(ri j=0,len=mid<<1;j<lim;j+=len){
			Cp w=(Cp){1,0};
			for(ri k=0;k<mid;++k,w=w*wn){
				Cp a0=a[j+k],a1=a[j+k+mid]*w;
				a[j+k]=a0+a1,a[j+k+mid]=a0-a1;
			}
		}
	}
	if(type==-1)for(ri i=0;i<lim;++i)a[i]=a[i]/lim;
}
int main(){
	n=read();
	for(ri i=1,x;i<=n;++i)x=read(),++a[x].x,++b[x*2].x,++c[x*3].x,mx=max(mx,x*3);
	init(),fft(a,1),fft(b,1),fft(c,1);
	for(ri i=0;i<lim;++i)a[i]=(a[i]*a[i]*a[i]-3*a[i]*b[i]+(Cp){2,0}*c[i])/6.0+(a[i]*a[i]-b[i])/2.0+a[i];
	fft(a,-1);
	for(ri i=0;i<=mx;++i){
		ll ans=(ll)(a[i].x+0.5);
		if(ans)cout<<i<<' '<<ans<<'\n';
	}
	return 0;
}

转载于:https://www.cnblogs.com/ldxcaicai/p/10367800.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值