1. 程式人生 > >2018.11.18 spoj Triple Sums(容斥原理+fft)

2018.11.18 spoj Triple Sums(容斥原理+fft)

傳送門
這次 f f t fft 亂搞居然沒有被卡常?
題目簡述:給你 n n 個數,每三個數 a

i , a j , a k ( i
< j < k ) a_i,a_j,a_k(i<j<k) 組成的所有和以及這些和出現的次數。
讀完題直接讓我聯想到了昨天寫過的一道 f
f t fft
優化點分治合併的題

,這不是差不多嘛?
只是這一次的容斥要麻煩一些。
我們令原數列轉化成的係數序列為 { a n } \{a_n\}
那麼如果允許重複答案就應該是 a n 3 a_n^3
然後展開式子。
我們需要容斥掉的就是有3個數相同的和有2個數相同的。
這個時候已經可以用 f f t fft 了。
但是還可以進一步優化。
如何優化?
觀察到有兩個數相同的其實可以再用一次容斥求出。
因為兩個數相同的個數是等於至少兩個數相同的個數扣去有三個數相同的個數的。
這樣可以優化更多常數。
剩下就是 f f t fft 的事啦。
程式碼:

#include<bits/stdc++.h>
#define ri register int
using namespace std;
inline int read(){
	int ans=0,w=1;
	char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
	while(isdigit(ch))ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar();
	return ans*w;
}
typedef long long ll;
const int N=(1<<17)+5,delta=20000;
const double pi=acos(-1.0);
struct Complex{
	double x,y;
	inline Complex operator+(const Complex&b){return (Complex){x+b.x,y+b.y};}
	inline Complex operator-(const Complex&b){return (Complex){x-b.x,y-b.y};}
	inline Complex operator*(const Complex&b){return (Complex){x*b.x-y*b.y,y*b.x+x*b.y};}
	inline Complex operator/(const double&b){return (Complex){x/b,y/b};}
}a[N],b[N],c[N];
int n,lim,tim,A[N],B[N],C[N],pos[N];
inline void fft(Complex *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){
		Complex wn=(Complex){cos(pi/mid),type*sin(pi/mid)};
		for(ri j=0,len=mid<<1;j<lim;j+=len){
			Complex w=(Complex){1,0};
			for(ri k=0;k<mid;++k,w=w*wn){
				Complex 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;
}
inline void init(){
	lim=1<<17,tim=17;
	for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}
int main(){
	freopen("lx.in","r",stdin);
	n=read(),init();
	for(ri i=1,val;i<=n;++i)val=read()+delta,++A[val],++B[val*2],++C[val*3];
	for(ri i=0;i<lim;++i)a[i].x=A[i],b[i].x=B[i],c[i].x=C[i];
	fft(a,1),fft(b,1);
	for(ri i=0;i<lim;++i)c[i]=a[i]*(a[i]*a[i]-(Complex){3,0}*b[i]);
	fft(c,-1);
	for(ri i=0;i<lim;++i){
		ll cnt=((ll)(c[i].x+0.5)+2*C[i])/6;
		if(cnt)cout<<i-3*delta<<" : "<<cnt<<'\n';
	}
	return 0;
}