FFT中的一个常见小问题
这里不细说FFT的内容,详细内容看这些就足以了解大概了
主要是对学习过程中一个容易困扰的小问题进行解释,以便于理解
-
用FFT将多项式的系数转换为点值时,原系数数组a最后存的是不同的点值,而不是只有第一个是点值
这一点最开始困扰了我很久
设 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n−1}x^{n−1} A(x)=a0+a1x+a2x2+...+an−1xn−1
则可将其移项 A ( x ) = ( a 0 + a 2 x 2 + . . . + a n − 2 x n − 2 ) + ( a 1 x + a 3 x 3 + . . . + a n − 1 x n − 1 ) A(x)=(a_0+a_2x^2+...+a_{n−2}x^{n−2})+(a_1x+a_{3}x^3+...+a_{n−1}x^{n−1}) A(x)=(a0+a2x2+...+an−2xn−2)+(a1x+a3x3+...+an−1xn−1)
a的下标为偶数的放在一起 A 1 ( x ) = a 0 + a 2 x + . . . + a n − 2 x n 2 − 1 A_1(x)=a_0+a_2x+...+a_{n−2}xn^{2−1} A1(x)=a0+a2x+...+an−2xn2−1
a的下标为奇数的放在一起 A 2 ( x ) = a 1 + a 3 x + . . . + a n − 1 x n 2 − 1 A_2(x)=a_1+a_3x+...+a_{n−1}xn^{2−1} A2(x)=a1+a3x+...+an−1xn2−1
则 A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 ) A(x)=A_1(x^2)+xA_2(x^2) A(x)=A1(x2)+xA2(x2)
注意此处为 x 2 x^2 x2所以有
A ( − x ) = A 1 ( x 2 ) − x A 2 ( x 2 ) A(-x)=A_1(x^2)-xA_2(x^2) A(−x)=A1(x2)−xA2(x2)
由于单位根的特殊性质,有
性质一 ω n k + n 2 = − ω n k ω_n^{k+\frac{n}{2}}=-ω_n^k ωnk+2n=−ωnk
性质二 ω n k = ω 2 n 2 k ω_n^k=ω_{2n}^{2k} ωnk=ω2n2k
所以才有了代码中的那两行for (int i=0;i<=mid-1;++i){ buf[i]=a[i]+w*a[i+mid]; buf[i+mid]=a[i]-w*a[i+mid]; w=w*wn; }
也就是说,我们可以由一个答案进而算出另外一个答案,这里可以理解为递推
所以当我们的递归递到最下面一层后往上走时每次都是将目前答案个数扩大两倍,而且这些答案是由不同的x算出来的,而且由于性质一,我们在计算过程中所用到的不同的 ω x ∗ k ω^{x*k} ωx∗k是没有问题的
最后附上板子#include <cstdio> #include <algorithm> #include <cmath> using namespace std; const int maxn = 4000006; const double pi = acos(-1.0); struct IO{ template<class T> IO operator >> (T &res){ res=0; char ch; bool flag=false; while ((ch=getchar())>'9'||ch<'0') flag|=ch=='-'; while (ch>='0'&&ch<='9') res=(res<<3)+(res<<1)+(ch^'0'),ch=getchar(); if (flag) res=~res+1; return *this; } }cin; struct complex { double x,y; complex (double xx=0,double yy=0) {x=xx,y=yy;} }; complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y);} complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y);} complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} int n,m,bit,len,val; int rev[maxn]; complex a[maxn],b[maxn],ans[maxn],buf[maxn]; //递归FFT void FFT (complex *a,int len,int on_off)//on_off=1 : FFT on_off=-1 : IFFT { if (len==1) return ; int mid=len/2; for (int i=0;i<=mid-1;++i) buf[i]=a[i*2],buf[i+mid]=a[i*2+1]; for (int i=0;i<=len;++i) a[i]=buf[i]; FFT(a,mid,on_off),FFT(a+mid,mid,on_off); complex wn=complex(cos(2*pi/len),on_off*sin(2*pi/len)),w(1,0); for (int i=0;i<=mid-1;++i){ buf[i]=a[i]+w*a[i+mid]; buf[i+mid]=a[i]-w*a[i+mid]; w=w*wn; } for (int i=0;i<=len;++i) a[i]=buf[i]; } //非递归FFT void FFT2 (complex *a,int len,int on_off)//on_off=1 : FFT on_off=-1 : IFFT { for (int i=0;i<=len-1;++i) if (i<rev[i]) swap(a[i],a[rev[i]]); for (int i=1;i<len;i<<=1){ complex wn=complex (cos(pi/i),on_off*sin(pi/i)); for (int j=0;j<len;j+=(i<<1)){ complex w(1,0); for (int k=0;k<i;++k){ complex u=a[j+k],t=w*a[i+j+k]; a[j+k]=u+t; a[i+j+k]=u-t; w=w*wn; } } } } int main () { cin>>n>>m; for (int i=0;i<=n;++i) cin>>val,a[i].x=val; for (int i=0;i<=m;++i) cin>>val,b[i].x=val; len=1; while (len<=n+m) ++bit,len<<=1; for (int i=0;i<=len-1;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); FFT2(a,len,1); FFT2(b,len,1); for (int i=0;i<=len;++i) ans[i]=a[i]*b[i]; FFT2(ans,len,-1); for (int i=0;i<=n+m;++i) printf("%d ",int(ans[i].x/len+0.5)); return 0; }