FFT入门
背景知识
FFT是快速求解两个多项式乘积的算法，时间复杂度为优秀的(O(nlogn)),主要想法是分治。
FFT = Fast离散傅里叶变换 = Fast DFT
多项式存在系数表示法和点值表示法，二者可以相互转换;
n次多项式的最高次为n,有(n+1)个系数，需要(n+1)个相异的点才能表示
FFT主要思路
- 预处理，填充高次为0,使得最终系数数目大于(N+M);(N次多项式和M次多项式)
- 系数表示转化成点值表示(代入的x坐标为(x = (w_n^0,w_n^1...w_n^{n-1}));这里的n是需要的点的个数;这个可以通过分治实现(O(nlogn))
- 点值表示法相乘,复杂度(O(n))
- 点值表示转化成系数表示,该问题可规约成步骤2
下面我们将重点放在2,3两点上
要点
多项式系数表示转化成点值表示
记多项式(A(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}), (A(x))是一个(n-1)次多项式，需要(n)个点才能表示。我们选取复数域上(w_n^k = e^{i*frac{2kpi}{n}} = cos(frac{2k pi}{n})+i sin(frac{2k pi}{n}) quad k=0,1,2,..n-1) 这样的i个点。
它们满足两条性质:
- (w_n^k = -w_n^{k+frac{n}{2}}quad) 逆时针旋转(pi)
- (w_n^k = w_{nd}^{kd}),假设(kd,nd)均为整数
上述两条性质在复数的单位圆上画一下即可知。
假设(n=2^t),实际上步骤0中预处理就会处理成(2^t);我们将(A(x)) 根据系数位置的就分成两部分,即:
(A1(x) = a_0+a_2x^2+..a_{n-2}^x{n-2} = sum_{i=0}^{i=frac{n}{2}-1}a_{2i}x^{2i})
(A2(x) = a_1+a_3x^3+...a_{n-1}x^{n-1} = sum_{i=0}^{i=frac{n}{2}-1}a_{2i+1}x^{2i+1} = xsum_{i=0}^{i=frac{n}{2}-1}a_{2i+1}x^{2i})
(A(x)=sum_{i=0}^{n-1}a_ix^i = A1(x)+A2(x))
注意到(A1(x), A2(x))虽然系数减少了一半，但需要带入的(x)依旧是(n)项;
假设我们将(x=w_n^k)代入(A1(x)) (代入(A2(x))是一样的),那么:
(A1(w_n^k)=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_n^{2ki}=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k=0,1,2....n-1),(根据性质2)
再一次变化,将k分为前半部分和后半部分处理:
(A1(w_n^k)=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k =0,1,2,frac{n}{2}-1), 前半部分
(A1(w_n^{k+frac{n}{2}})=-sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k =0,1,2,frac{n}{2}-1) 后半部分，(运用性质1)
由此可知我们在求得(A1(w_n^k))的同时已经把(A1(w_n^{k+frac{n}{2}}))求了出来，这样求解(A1(w_n^k))的规模相对于(A(x))将为了一半，多么nice的结果((A2(x))同理)。
所以 (T(n) = 2*(n/2)+O(n)),因为(x)有(n)种取值，所以通过(A1(x))和(A2(x))求得(A(x))需要(O(n)) 时间
递归写法
void FFT(Complex* a,int len,int opt){
if(len==1)
return;
Complex* a0=new Complex[len/2];
Complex* a1=new Complex[len/2];
for(int i=0;i<len;i+=2){
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}
FFT(a0,len/2);
FFT(a1,len/2);
Complex wn(cos(2*Pi/len),opt*sin(2*Pi/len));
Complex w(1,0);
for(int i=0;i<(len/2);i++){
a[i]=a0[i]+w*a1[i];
a[i+len/2]=a0[i]-w*a1[i];
w=w*wn;
}
delete[] a0;
delete[] a1;
}
点值转换成系数表示(有人称作IDFT)
思考一下系数表示到点值的矩阵形式:
我们记作 (W cdot A = Y); 那么从点值到系数的变换即为 (A = W^{-1}cdot Y)。
下面不加证明的给出 (W^{-1}=frac{1}{n}[w^{-ij}]), (w^{ij})表示(W) i行j列的元素。即W的逆的W原矩阵的每个复数取逆后的矩阵的(frac{1}{n}). 这样步骤3就规约到了步骤1
FFT 优化
从递归形式到迭代形式,直接给出详细注释的代码，可在luogu3803交题
迭代版本
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=1e7+10;
const double Pi=acos(-1.0);
struct complex{
double x,y;
complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[MAXN],b[MAXN];
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;
int l,r[MAXN];
int limit=1;
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)//待合并区间的中点,其实(mid<<1)是要求的DFT(某一层)的长度
{
complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根 本质是 2*PI/(2*mid)
for(int R=mid<<1,j=0;j<limit;j+=R)//R实际是步长，即同一层每组要处理的起始位置
{
complex w(1,0);//幂 单位1
for(int k=0;k<mid;k++,w=w*Wn)//枚举左半部分
{
complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应(只是装x),实际上什么有没有就是减少一次y的计算(复数乘法)
A[j+k]=x+y; // A[k]
A[j+mid+k]=x-y; //mid即位推导中的2/n;此处可以看作 A[k+2/n]
}
}
}
}
int main()
{
//n次多项式由(n+1)个系数,有n个复数解，最高次为n
scanf("%d",&N);scanf("%d",&M);
for(int i=0;i<=N;i++) scanf("%lf",&a[i].x);//复数域的实数部分
for(int i=0;i<=M;i++) scanf("%lf",&b[i].x);
while(limit<=N+M) limit<<=1,l++; //最高次为(N+M) 所以至少需要(N+M+1)个点值;为了fft的方便扩展成 2^l(2^l>(N+M)) 高次皆为0
for(int i=0;i<limit;i++)
r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ; //位逆置 r[i>>1]>>1 能反转前面的l-1 位； 然后将末尾的位或到首位即可
// 在原序列中 i 与 i/2 的关系是 ： i可以看做是i/2的二进制上的每一位左移一位得来
// 那么在反转后的数组中就需要右移一位，同时特殊处理一下复数
fft(a,1);// fft保留的是x= w_{limit}^1, w_{limit}^2 .... w_{limit}^{limit}是的取值
fft(b,1);
for(int i=0;i<=limit;i++) a[i]=a[i]*b[i]; //点值表示法的优秀乘法
fft(a,-1); //点值表示法到系数表示法，会涉及到逆矩阵，所以fft后还有一个1/limit (limit 实际上是矩阵的行列式(n*n的满秩矩阵才有行列式))
for(int i=0;i<=N+M;i++)
printf("%d ",(int)(a[i].x/limit+0.1));//直接用实数域
return 0;
}