BZOJ3771:Triple FFT+容斥原理
Feb 04, 2015 09:09:19 AM
思路:
首先第一次看的时候认为是神题.
再一次看的时候卧槽这不是FFT裸题么.
然后写的时候发现有点坑...
(FFT都告诉你了难道还不会么)
#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef double f2;
static const f2 pi=acos(-1.0);
struct Complex{
f2 u,v;
Complex(){}
Complex(f2 _u,f2 _v):u(_u),v(_v){}
Complex operator+(const Complex&B)const{return Complex(u+B.u,v+B.v);}
Complex operator-(const Complex&B)const{return Complex(u-B.u,v-B.v);}
Complex operator*(const Complex&B)const{return Complex(u*B.u-v*B.v,u*B.v+v*B.u);}
Complex operator*(const f2&p)const{return Complex(p*u,p*v);}
Complex operator/(const f2&p)const{return Complex(u/p,v/p);}
inline void operator=(const f2&p){u=p,v=0;}
inline void operator+=(const Complex&B){u+=B.u,v+=B.v;}
inline void operator-=(const Complex&B){u-=B.u,v-=B.v;}
inline void operator*=(const Complex&B){f2 uu=u*B.u-v*B.v,vv=u*B.v+v*B.u;u=uu,v=vv;}
inline void operator*=(const f2&p){u*=p,v*=p;}
inline void operator/=(const f2&p){u/=p,v/=p;}
};
inline int revbit(int x,int b){
int res=0;for(int i=0;i<b;++i)if((x>>i)&1)res+=1<<(b-1-i);return res;
}
inline void fft(Complex A[],int n,int rev){
register int i,j,k;
int bit=0,tmp=n;for(;tmp^1;tmp>>=1)++bit;
int ano;for(i=0;i<n;++i){ano=revbit(i,bit);if(i<ano)swap(A[i],A[ano]);}
Complex wn,w,t;
for(k=2;k<=n;k<<=1)
for(wn=Complex(cos(2*pi/k),rev*sin(2*pi/k)),i=0;i<n;i+=k)
for(w=1,j=0;j<k/2;++j,w*=wn)
t=w*A[i+j+k/2],A[i+j+k/2]=A[i+j]-t,A[i+j]+=t;
if(rev==-1)for(i=0;i<n;++i)A[i]/=n;
}
#define N 131072
inline void mul(Complex A[],Complex B[],Complex C[]){
static Complex sa[N],sb[N];
memcpy(sa,A,N*sizeof(Complex)),memcpy(sb,B,N*sizeof(Complex));
fft(sa,N,1),fft(sb,N,1);
for(int i=0;i<N;++i)C[i]=sa[i]*sb[i];
fft(C,N,-1);
}
Complex A[N],A2[N],A3[N];
Complex mul2[N],mul3[N],mul3_2[N];
long long res[N];
long long up(const f2&p){
return(long long)(p+.5);
}
int main(){
int n;scanf("%d",&n);register int i,j;
int x;while(n--)scanf("%d",&x),A[x]=1,A2[2*x]=1,A3[3*x]=1;
mul(A,A,mul2);
mul(mul2,A,mul3);
mul(A,A2,mul3_2);
for(int i=0;i<N;++i)mul2[i]-=A2[i],mul2[i]*=.5;
for(int i=0;i<N;++i)mul3_2[i]-=A3[i];
for(int i=0;i<N;++i)mul3[i]-=A3[i]+mul3_2[i]*3,mul3[i]/=6;
for(int i=0;i<N;++i)res[i]+=up(A[i].u)+up(mul2[i].u)+up(mul3[i].u);
for(int i=0;i<N;++i)if(res[i]!=0)printf("%d %lld\n",i,res[i]);
return 0;
}