bzoj3527: [Zjoi2014]力

当时就是从这题得知还有FFT这个东西的。

在上一篇博文里欠了一个FFT的高效实现。这里补上。

显然我们不会选择递归，但是我们要从递归出发。

递归中对A[0]和A[1]两个子问题的点值表达式分别求值。

就是把系数分成了两块。次数降为一半。

如何做到从递归到自底向上直接求？这用到一个很聪明的方法。

之前我们得到

        A(wn^(k+n/2))=A[0]((wn^k)^2)-wn^kA[1]((wn^k)^2)

        A(wn^k)=A[0]((wn^k)^2)+wn^kA[1]((wn^k)^2)
在递归实现中这样可行，但是自底向上如何寻找相应的A[0]和A[1]呢
我们对于(a0,a1,a2,a3,a4,a5,a6,a7)的系数输入尝试递归。

最后得到a0,a4,a2,a6,a1,a5,a3,a7的递归底层。

按照二进制，它们分别是000 100 010 110 001 101 011 111

没规律？没关系，把它们的码位取反。

就是000 001 010 011 100 101 110 111

十进制就是0,1,2,3,4,5,6,7

很神奇吧。所以我们对于相应的系数只要对其下标做码位翻转，就可以预测其在递归树底层的位置。

然而递归树是一层层排列的，且总是由前一层得到后一层。

于是我们就可以做到自底向上实现FFT了。

且这一层的相应位置的取值是由上一层相应位置求出的，不会出现覆盖问题导致的错误。

比如当前层的[0,4]由其递归层的[0,2]和[3,4]求出。

当前层0和3的位置可以由上一层的0和3的位置求出，2和4亦是如此，多么奇妙的巧合！

所以我们完全可以全程使用同一个数组。 注意，是复数数组。

至于复数怎么编。。自己打个结构体再重载一下运算符就可以了。

STL的complex太慢了。。

于是回到这道题，只要理解从多项式乘法到卷积的跨越，就可以了。

事实上多项式乘法本来就是个卷积的过程。

怎么理解？自己把多项式乘法拆开来试试。

#include<iostream>
#include<cstdio>
#include<cmath>

const int MN=1000010;
const double pi=acos(-1.0);

struct complex{
double r,i;
complex(double _r=0,double _i=0){
r=_r; i=_i;
}
}x,y,a[MN],b[MN],c[MN],A[MN];
int rev[MN],dig[MN],n;
double X[MN],ans[MN];
int N,L;

complex operator+(complex a,complex b){
return complex(a.r+b.r,a.i+b.i);
}
complex operator-(complex a,complex b){
return complex(a.r-b.r,a.i-b.i);
}
complex operator*(complex a,complex b){
return complex(a.r*b.r-a.i*b.i,a.r*b.i+b.r*a.i);
}

void FFT(complex a[],int f){
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){
complex wn(cos(2*pi/i),f*sin(2*pi/i));
for(int k=0;k<N;k+=i){
complex w(1);
for(int j=0;j<(i>>1);j++) x=a[k+j],y=w*a[k+j+(i>>1)],a[k+j]=x+y,a[k+j+(i>>1)]=x-y,w=w*wn;
}
}if (f==-1) for(int i=0;i<N;i++) a[i].r/=N;
}

int main(){
scanf("%d",&n);
for(N=1;N<n;N<<=1,L++); N<<=1,L++;
for(int i=0;i<N;i++){
for(int t=i,len=0;t;t>>=1) dig[len++]=t&1;
for(int j=0;j<L;j++) rev[i]=(rev[i]<<1)+dig[j];
}for(int i=0;i<n;i++) scanf("%lf",&X[i]);
for(int i=0;i<n;i++) a[i]=complex(X[i]);
for(int i=1;i<n;i++) b[i]=complex(1./i/i);
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++) ans[i]=c[i].r;
for(int i=0;i<N;i++) a[i]=b[i]=complex();
for(int i=0;i<n;i++) a[i]=complex(X[n-i-1]);
for(int i=1;i<n;i++) b[i]=complex(1./i/i);
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++) ans[i]-=c[n-i-1].r;
for (int i=0;i<n;i++) printf("%0.3lf\n",ans[i]);
return 0;
}