【BZOJ】【P3527】【ZJOI2014】【力】【题解】【FFT】

来源：互联网发布：淘宝客服晚班可在家上编辑：程序博客网时间：2024/05/20 03:38

传送门 www.lydsy.com/JudgeOnline/problem.php?id=3527

好久没写FFT了……全忘干净了……推完式子不会写FFT了……

============================================================

更正：倒数第三个式子应为 E[i]=sigma( f(j)g(j-i) ) - sigma( f(n-j-1)g(j-i) )

下面两个式子也有类似错误.

============================================================

大概就是这样了，经过各种转换成了卷积的形式……然后FFT，重新学了一遍FFT，YM&学习了 pyc神犇的FFT，

Code：

#include<cmath>#include<cstdio>#include<cstring>#include<iostream>#include<algorithm>using namespace std;const int maxn=1000010;int n,N,L;int rev[maxn];int dig[maxn];double p[maxn];struct cp{double r,i;cp(double _r=0,double _i=0):r(_r),i(_i){}cp operator+(cp x){return cp(r+x.r,i+x.i);}cp operator-(cp x){return cp(r-x.r,i-x.i);}cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}};cp a[maxn],b[maxn],c[maxn],A[maxn],x,y;void FFT(cp a[],int flag){for(int i=0;i<N;i++)A[i]=a[rev[i]];for(int i=0;i<N;i++)a[i]=A[i];for(int i=2;i<=N;i<<=1){cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));for(int k=0;k<N;k+=i){cp w(1,0);for(int j=0;j<i/2;j++){x=a[k+j];y=w*a[k+j+i/2];a[k+j]=x+y;a[k+j+i/2]=x-y;w=w*wn;  }}}if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;}double anss[maxn];int main(){scanf("%d",&n);for(int i=0;i<n;i++)scanf("%lf",&p[i]);for(L=0,N=1;N<n;N<<=1,L++);L++;N<<=1;for(int i=0;i<N;i++){int len=0;for(int t=i;t;t>>=1)dig[len++]=t&1;for(int j=0;j<L;j++)rev[i]=rev[i]*2+dig[j];}for(int i=0;i<n;i++)a[i]=cp(p[i],0);for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);FFT(a,1);FFT(b,1);for(int i=0;i<N;i++)c[i]=a[i]*b[i];FFT(c,-1);for(int i=0;i<n;i++)anss[i]=c[i].r;memset(a,0,sizeof(a));memset(b,0,sizeof(b));for(int i=0;i<n;i++)a[i]=cp(p[n-i-1],0);for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);FFT(a,1);FFT(b,1);for(int i=0;i<N;i++)c[i]=a[i]*b[i];FFT(c,-1);for(int i=0;i<n;i++)anss[i]-=c[n-i-1].r;for(int i=0;i<n;i++)printf("%.9f\n",anss[i]);return 0;}

0 0