FFT—快速傅里叶变换学习小记

我们要干什么

FFT能够快速地解决多项式乘法。额···好像也只能干这个事。
一个特殊的例子：高精度乘法。

问题简述

我们记一个多项式A(x)的次数界为n，则
A(x)=∑n−1i=0ai∗xi，其中a为系数，x为变量。
注意最高次系数为n-1，不是n。实际上次数界就是有多少项系数。
两个多项式相乘，我们一般记为C(x)=A(x)B(x)，C的次数界显然是A,B的和-1。一般来说，乘起来要O(n2)

FFT的中心思想

转换思路

普通的相乘方法：cj=∑ji=0ai∗bj−i
提出概念： 点值，插值。

点值

一个次数界为n的多项式A的点值表达为n个点值对所组成的集合。形如
{(x0,y0),(x1,y1)...(xn−1,yn−1)}
其中yi=A(xi)，x各不相同。
比如多项式A(x)=x2−3x+2的点值表达可以为，{(-1,6),(0,2),(5,12)}。

插值

插值运算为点值运算的逆运算，如果给定一个n个点值对的点值表达，我们可以确定唯一一个次数界为n的多项式。就比如通过上面的三个点，可以求出A(x)的各项系数。

综合两者进行多项式乘法

确定一组x，可以得到A,B的点值表达。
A:{(x0,y0),(x1,y1)...(xn−1,yn−1)}
B:{(x0,y′0),(x1,y′1)...(xn−1,y′n−1)}
那么可以知道C的点值表达为
C:{(x0,y0∗y′0),(x1,y1∗y′1)...(xn−1,yn−1∗y′n−1)}
不太理解的话，可以把y0和y′0的表达式写出来相乘，发现y0∗y′0的表达式和前者确实相同
通过插值运算，就可以得到多项式C的系数了。
注意：A,B次数界均为n的话，C的次数界为2n-1(为了方便直接弄成2n)，所以x集合大小至少为2n，不然插值失败。

流程

对于次数界均为n的两个多项式A,B
1.点值运算，构造长度为2n的点值表达；
2.每个点的y相乘，算出C的点值表达；
3.插值运算，通过点值表达求出系数。
然而暴力来做时间复杂度仍为O(n2)，没有区别。
其实我们可以通过选特殊的x集合，来优化时间，我们选的数叫做n次单位复数根。

N次单位复数根及其性质

n次单位复数根是满足wn=1的复数w。那么我们知道给定了n，n次单位复数根恰好有n个。对于根k=0,1…n-1，这些根是e2πik/n(i为虚数单位)，为了解释，我们需要用到复数的指数形式的定义。
eiu=cos(u)+isin(u)。
这一条是欧拉定理(之一)，实际上，我们发现eiu的泰勒展开和cos(u)的泰勒展开加上isin(u)的泰勒展开是一样的。所以两边相等。
我们来看个图，直观地发现规律。

可以发现，他们均匀分布在复平面(以实数域为x轴，虚数域i为y轴的平面)的单位圆上。跟学三角函数的时候那些东西是一样的。

性质1，主n次单位根

我们称wn=w1n为主n次单位根，注意到每个n次单位复数根都是由旋转得到的，每次旋转有一个固定的角度，可见win=wi−1n∗wn。
n次单位复数根可视为公比为wn的等比数列。

性质2，群的性质

由于单位根们可表达为一对三角函数，那么他也是具有周期的
w0n=wnn=1,wkn=wk%nn
而由性质1得，一个wkn实际上可以看成(wn)k，那么wkn∗wjn=wj+kn
又可以得到(wkn)2=w2kn

两个引理

消去引理

wdkdn=wkn，可以从图象角度理解，他们所对应的点是相同的。
推论：
wn/2n=w2=−1
w2kn=wkn/2

折半引理

若n为>0的偶数，那么n次单位复数根的平方的集合，等于n/2次单位复数根的集合。
这个也好理解，可以知道(wkn)2=(wk+n/2n)2，这个由上面的一些性质可以推出来。图象上理解：n次单位复数根的平方，就对应着三角函数里的角度乘了2。

求和引理

对于任何整数n>=1和不能被n整除的非负整数k（负数可以通过群的性质变为正数），有
∑n−1j=0(wkn)j=0
等比数列求和，化出来发现末项和首项都是1，所以是0。k为n的时候分母为0，所以不成立。

FFT及其关键算法

DFT——离散傅里叶变换

我们的目标：计算次数界为n的多项式A(x)=∑jaj∗xj在n次单位复数根上的值。
定义结果y
yk=A(wkn)=∑n−1j=0aj∗wjkn
y即为a的离散傅里叶变换。
可记为y=DFTn(a)

FFT——快速傅里叶变换

为了充分利用n次单位复数根，我们令n为2的幂，这样可以在O(nlog2n)时间内求出DFTn(a)。为了这样做，我们用分治的策略。

分治策略

如何求出单个数x的函数值A(x)？我们定义两个多项式
A0(x)=a0+a2x+a4x2⋅⋅⋅an−2xn/2
A1(x)=a1+a1x+a3x2⋅⋅⋅an−1xn/2
即一个全是偶数编号的系数，一个全是奇数的。我们发现这两个东西的次数是n/2，变小了一半。
可得A(x)=A0(x2)+xA1(x2)，那么再使用点值相乘。

阶段性总结

原问题：求A(x)在n次单位复数根上的函数值。
转化成：求A0(x)和A1(x)在n/2次单位复数根上的值。
可以发现不同的自变量少了一半，而每个自变量出现2次。
于是，我们递归地对n/2的多项式A0(x)与A1(x)在n/2个n/2次单位复数根进行求值。

伪代码

看到最后的处理，我们对于(wkn)2相同的一起处理，而平方相同的两个n次单位复数根恰好互为相反数。
运行时间显然是O(nlogn)的。
这样我们解决了点值运算。可以处理出C的点值表达了！

插值运算

这个其实跟点值运算差不多。
我们如果把点值运算写成矩阵方程的形式，可以得到表达式y=Vna。

只是把∑换个形式写而已。
那么可以知道a=y∗V−1n。也就是说，只要得到逆矩阵，再像上面做类似的点值运算，就能得到系数了！

逆矩阵

定理，V−1n[j,k]=w−jkn/n。
证明很麻烦，考虑证明Vn∗V−1n=单位矩阵，即(i,i)处为1，其他为0的矩阵。
证明：
[Vn∗V−1n]j,j′=∑n−1k=0wk∗(j′−j)n/n，如果j’=j，那后面为1，否则根据求和引理，为0。
那么，只需要把原来fft里面的a和y互换，wn变成w−1n，再进行fft，结果再除n，就能得到系数了！我们称这个为逆DFT。

实现的优化

我们不可能像伪代码一样赋值数组嘛···考虑给数组重新分配位置，使得最后算出来又是原来的顺序。
观察这幅图

最下面一层为我们数组应该保存的顺序。可以发现，他们的位置，刚好就是原下标在二进制下的翻转(开头为0，结尾为次数界n在二进制下的最高位的翻转)。
然后对于上面几层，我们每个组里面有多个值。

发现如果按照刚刚那个顺序存，每一层转移到下一层的时候，一个A(x0)要提取的A0(x0)和A1(x0)的间隔都是一样的。
设一行分为cnt组，即相同的x0出现cnt次，共有m个不同的x0，可以发现转移间隔为m/2。而平方相等的n次单位复数根的间隔也是m/2，这就很优美了。
这个大概理解一下就好，反正代码很简洁明了。

代码

#include<cstdio>#include<algorithm>#include<cstring>#include<bitset>#include<cmath>#define fo(i,j,k) for(i=j;i<=k;i++)#define fd(i,j,k) for(i=j;i>=k;i--)typedef long long ll;typedef double db;using namespace std;const int N=400005;const db pi=acos(-1),eps=1e-5;struct vec{    db a,b;    vec (db ax=0,db bx=0) {a=ax,b=bx;}};vec operator +(vec a,vec b){    return vec(a.a+b.a,a.b+b.b);}vec operator -(vec a,vec b){    return vec(a.a-b.a,a.b-b.b);}vec operator *(vec a,vec b){return vec(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a);}int Log,m,len,i,n,ans[N];vec t[N],a[N],b[N];void DFT(vec *a,int n,int sig){    int i,j,k,half,siz;    vec u,v,w;    fo(i,0,n-1)    {        int p=0;        for(int pos=0,tmp=i;pos<Log;pos++,tmp/=2) p=(p<<1)+(tmp&1);        t[p]=a[i];    }    for(int m=2; m<=n; m*=2)    {         int half = m/2;         for(int i=0; i<half; i++)        {             vec w( cos(i*sig*pi/half), sin(i*sig*pi/half));            for(int j=i; j<n; j+=m)            {                 k=j+half;                u=t[j];v=w*t[k];                t[j]=u+v;                t[k]=u-v;            }        }    }    fo(i,0,n-1) a[i]=t[i];}void FFT(vec *a,vec *b){    DFT(a,len,1);    DFT(b,len,1);    fo(i,0,len-1) a[i]=a[i]*b[i];    DFT(a,len,-1);    fo(i,0,len-1) a[i].a/=(db)len;}int max2(int x){    int ret;    for(ret=1;ret<x;ret<<=1);    return ret<<1;}int ci(int x) {return (x==1)?0:ci(x/2)+1;}int main(){    freopen("taxi.in","r",stdin);    scanf("%d %d",&n,&m);    n++;m++;    len=max2(max(n,m));    Log=ci(len);    fo(i,0,n-1) scanf("%lf",&a[i].a);    fo(i,0,m-1) scanf("%lf",&b[i].a);    FFT(a,b);    fo(i,0,n+m-1) ans[i]=round(a[i].a+eps);    fo(i,0,n+m-2) printf("%d ",ans[i]);}//注意，由于是实数运算，所以给出系数在10^5以上时很容易出现误差

总结

中心思想：点值运算与插值运算
优化时间复杂度：选取n次单位复数根
难点：n次单位复数根的各种性质
重点：FFT的分治思想，DFT与逆DFT
要背的地方：FFT实现优化