FFT快速傅里叶
快速傅里叶
（这个博客主要帮助自己记着FFT这个算法，并不是讲解用的QAQ）
定义：
现在有两个多项式：
\(f(x)=a_1+a_2x+a_3x^2+...+a_nx^{n-1}\)
\(g(x)=b_1+b_2x+g_3x^2+...+g_mx^{m-1}\)
加入我们计算 \(f(x)*g(x)\) 的表达式时，暴力算法的复杂度为 \(O(n^2)\) ,我们这时候就需要引入快速傅里叶算法来解决这个问题。
思想：
1.初始想法
枚举出 \(n+1\) 个点，计算出其分别在 \(f(x),g(x)\) 上对应的坐标，利用拉格朗日插值去求解。
2.减少时间
虽然这样看上去简单了，但是时间复杂度还是 \(O(n^2)\) ，因此我们可以观察原来每个函数的性质：如果这个函数是偶函数，那么枚举 \(n+1\) 个点，根据偶函数性质我们就可以得到 \(2n+2\) 个点的坐标，这样就可以减少枚举的次数了
3.转化多项式
例如
\(x^5+x^4+3x^3+2x^2+5x\)
可以转化成：
\((x^4+2x^2)+x(x^4+3x^2+5)\)
这样，每一个括号内就变成了偶函数。
将前后两个式子分别定义为 \(a,b\) ，那么 \(x>0\) 时，\(f(x)=a+b\)，反之，则是 \(f(x)=a-b\)
4.利用复数进行枚举
根据就当是我已经会了的单位圆的知识可以得到：\(n\) 等分的单位元可以利用复数得到不同的值，从而很好的减少了枚举的时间复杂度。
同时，单位元还有个性质，就是前后可以同时算（QAQ懒得证明了）。
后面：
\(\large{F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2}})\)
前面：
\(\large{F(ω^{k}_n)=FL(ω^{k}_{n/2})+ω^k_nFR(ω^{k}_{n/2}})\)
前后除了正负号其他都一样。
这样进行分治就可以啦啦啦啦啦啦。
5.保证limit是2的整数次幂
为什么要这样呢？我们可以想一想，如果分的时候不是2的整数次幂，就会分乱，因此要定义一个 \(limit\) 进行左右划分。
代码：
while(limit<=n+m) limit<<=1,l++;
6.复数计算：
利用 \(STL\) 中的 $complex $ 代码，重定义三维运算。
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);
}
7.反演:
就是传的参数不一样。
8.蝴蝶优化：
我们在进行多项式拆解的时候还是需要去一个一个算，那么怎么快速拆解呢?
根据网上的理论：
一个数的二进制表达反过来就是其应该在的位置的二进制表达
例如5的二进制是100,其应该在的位置就是001,也就是1。
就有接下来一串代码：
fft中：
for(int i=0;i<limit;i++)//每对数只反转一次
if(i<r[i]) swap(A[i],A[r[i]]);
main中：
for(int i=0;i<limit;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
模板：
#include<cstdio>
#include<cmath>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const double PI=acos(-1.0);
const int N=2e6+5;
int n,m;
int limit,l,r[N];
struct complex{
double x,y;
complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
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);
}
void fft(complex *A,int type){
for(int i=0;i<limit;i++)//每对数只反转一次
if(i<r[i]) swap(A[i],A[r[i]]);
for(int mid=1;mid<limit;mid<<=1){//等待合并区间的中点
complex Wn( cos(Pi/mid),type*sin(Pi/mid));//单位根
for(int R=mid<<1,j=0;j<limit;j+=R){
complex w(1,0);
for(int k=0;k<mid;k++,w=w*Wn){//枚举左半部分
complex x=A[j+k],y=w*A[j+mid+k];
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
int main()
{
cin>>n>>m;
for(int i=1;i<=n;i++) cin>>a[i].x;
for(int i=1;i<=m;i++) cin>>b[i].x;
while(limit<=n+m) limit<<=1,l++;
for(int i=0;i<limit;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
fft(a,1),fft(b,1);
for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;i++)
.......//此时求出来的就是多项式每一位的系数
return 0;
}
注意事项：
1.不能用万能头......
2.暂时没发现QAQ
