FFT中的一个常见小问题(递推式)
FFT中的一个常见小问题
这里不细说FFT的内容,详细内容看这些就足以了解大概了
小学生都能看懂的FFT!!!
FFT详解
补充——FFT中的二进制翻转问题
主要是对学习过程中一个容易困扰的小问题进行解释,以便于理解
用FFT将多项式的系数转换为点值时,原系数数组a最后存的是不同的点值,而不是只有第一个是点值
这一点最开始困扰了我很久
设A(x)=a0+a1x+a2x2+...+an−1xn−1
则可将其移项A(x)=(a0+a2x2+...+an−2xn−2)+(a1x+a3x3+...+an−1xn−1)
a的下标为偶数的放在一起A1(x)=a0+a2x+...+an−2xn−1
a的下标为奇数的放在一起A2(x)=a1+a3x+...+an−1xn−1
则A(x)=A1(x2)+xA2(x2)
注意此处为x2所以有
A(-x)=A1(x2)-xA2(x2)
由于单位根的特殊性质,有
性质一 ωnk+n/2+-ωnk
性质二 ωnk=ω2n2k
所以才有了代码中的那两行
1 for (int i=0;i<=mid-1;++i){ 2 buf[i]=a[i]+w*a[i+mid]; 3 buf[i+mid]=a[i]-w*a[i+mid]; 4 w=w*wn; 5 }
也就是说,我们可以由一个答案进而算出另外一个答案,这里可以理解为递推
所以当我们的递归递到最下面一层后往上走时每次都是将目前答案个数扩大两倍,而且这些答案是由不同的x算出来的,而且由于性质一,我们在计算过程中所用到的不同的$ω^{x*k}$是没有问题的
最后附上板子
原题 洛谷P3803 【模板】多项式乘法(FFT)
1 #include <cstdio> 2 #include <algorithm> 3 #include <cmath> 4 using namespace std; 5 const int maxn = 4000006; 6 const double pi = acos(-1.0); 7 struct IO{ 8 template<class T> 9 IO operator >> (T &res){ 10 res=0; 11 char ch; 12 bool flag=false; 13 while ((ch=getchar())>'9'||ch<'0') flag|=ch=='-'; 14 while (ch>='0'&&ch<='9') res=(res<<3)+(res<<1)+(ch^'0'),ch=getchar(); 15 if (flag) res=~res+1; 16 return *this; 17 } 18 }cin; 19 struct complex { 20 double x,y; 21 complex (double xx=0,double yy=0) {x=xx,y=yy;} 22 }; 23 complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y);} 24 complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y);} 25 complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} 26 int n,m,bit,len,val; 27 int rev[maxn]; 28 complex a[maxn],b[maxn],ans[maxn],buf[maxn]; 29 //递归FFT 30 void FFT (complex *a,int len,int on_off)//on_off=1 : FFT on_off=-1 : IFFT 31 { 32 if (len==1) return ; 33 int mid=len/2; 34 for (int i=0;i<=mid-1;++i) buf[i]=a[i*2],buf[i+mid]=a[i*2+1]; 35 for (int i=0;i<=len;++i) a[i]=buf[i]; 36 FFT(a,mid,on_off),FFT(a+mid,mid,on_off); 37 complex wn=complex(cos(2*pi/len),on_off*sin(2*pi/len)),w(1,0); 38 for (int i=0;i<=mid-1;++i){ 39 buf[i]=a[i]+w*a[i+mid]; 40 buf[i+mid]=a[i]-w*a[i+mid]; 41 w=w*wn; 42 } 43 for (int i=0;i<=len;++i) a[i]=buf[i]; 44 } 45 //非递归FFT 46 void FFT2 (complex *a,int len,int on_off)//on_off=1 : FFT on_off=-1 : IFFT 47 { 48 for (int i=0;i<=len-1;++i) 49 if (i<rev[i]) swap(a[i],a[rev[i]]); 50 for (int i=1;i<len;i<<=1){ 51 complex wn=complex (cos(pi/i),on_off*sin(pi/i)); 52 for (int j=0;j<len;j+=(i<<1)){ 53 complex w(1,0); 54 for (int k=0;k<i;++k){ 55 complex u=a[j+k],t=w*a[i+j+k]; 56 a[j+k]=u+t; 57 a[i+j+k]=u-t; 58 w=w*wn; 59 } 60 } 61 } 62 } 63 int main () 64 { 65 cin>>n>>m; 66 for (int i=0;i<=n;++i) cin>>val,a[i].x=val; 67 for (int i=0;i<=m;++i) cin>>val,b[i].x=val; 68 len=1; 69 while (len<=n+m) ++bit,len<<=1; 70 for (int i=0;i<=len-1;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); 71 FFT2(a,len,1); 72 FFT2(b,len,1); 73 for (int i=0;i<=len;++i) ans[i]=a[i]*b[i]; 74 FFT2(ans,len,-1); 75 for (int i=0;i<=n+m;++i) printf("%d ",int(ans[i].x/len+0.5)); 76 return 0; 77 }
如仍有问题或有其它问题可在下方指出,博主看到后会尽力解决,Thanks♪(・ω・)ノ