[bzoj3771]Triple——生成函数+容斥原理+FFT

题目大意：

“这把斧头，是不是你的？”

“这把斧头，是不是你的？”

“这把斧头，是不是你的？”

“你看看你现在的样子，真是丑陋！”

思路：

/*=======================================
* Author : ylsoi
* Time : 2019.1.30
* Problem : bzoj3771
* E-mail : ylsoi@foxmail.com
* ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
freopen("bzoj3771.in","r",stdin);
freopen("bzoj3771.out","w",stdout);
}

_=0; T fl=1; char ch=getchar();
for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
_*=fl;
}

const int maxn=240000+10;
const double pi=acos(-1);

struct cp{
double x,y;
cp(double xx=0,double yy=0){
x=xx,y=yy;
}
};
cp operator + (cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator - (cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator * (cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp operator / (cp a,double b){return cp(a.x/b,a.y/b);}
cp operator * (cp a,double b){return cp(a.x*b,a.y*b);}

cp g[maxn],ig[maxn];
int lim,cnt,dn[maxn];

void fft(cp *A,int ty){
REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
for(int len=1;len<lim;len<<=1){
cp w= ty==1 ? g[len<<1] : ig[len<<1];
for(int L=0;L<lim;L+=len<<1){
cp wk=cp(1,0);
REP(i,L,L+len-1){
cp u=A[i],v=A[i+len]*wk;
A[i]=u+v;
A[i+len]=u-v;
wk=wk*w;
}
}
}
}

int n,ans[maxn];
cp a[maxn],b[maxn],c[maxn],y[maxn],z[maxn];

void init(){
lim=1,cnt=0;
while(lim<=12e4)lim<<=1,++cnt;
if(!cnt)cnt=1;
REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
g[lim]=cp(cos(pi*2.0/lim),sin(pi*2.0/lim));
ig[lim]=cp(cos(pi*2.0/lim),-sin(pi*2.0/lim));
for(int i=lim>>1;i;i>>=1){
g[i]=g[i<<1]*g[i<<1];
ig[i]=ig[i<<1]*ig[i<<1];
}
}

void work(){
int t;
REP(i,0,lim-1)if(a[i].x!=0){
b[i*2].x+=a[i].x*a[i].x;
c[i*3].x+=a[i].x*a[i].x*a[i].x;
}

fft(a,1),fft(b,1),fft(c,1);
REP(i,0,lim-1){
y[i]=(a[i]*a[i]-b[i])/2;
z[i]=(a[i]*a[i]*a[i]-a[i]*b[i]*3+c[i]*2)/6;
}
fft(a,-1),fft(y,-1),fft(z,-1);
REP(i,0,lim-1){
a[i].x/=lim;
y[i].x/=lim;
z[i].x/=lim;
ans[i]+=(int)(a[i].x+0.5)+(int)(y[i].x+0.5)+(int)(z[i].x+0.5);
}
REP(i,0,lim-1)if(ans[i])
printf("%d %d\n",i,ans[i]);
}

int main(){
File();
init();
work();
return 0;
}


posted @ 2019-01-31 14:37  ylsoi  阅读(138)  评论(0编辑  收藏  举报