转至:http://www.cnblogs.com/zwfymqz/p/8244902.html
快速傅里叶变换(FFT)详解
本文只讨论FFT在信息学奥赛中的应用
文中内容均为个人理解,如有错误请指出,不胜感激
前言
先解释几个比较容易混淆的缩写吧
DFT:离散傅里叶变换—>O(n2)O(n2)计算多项式乘法
FFT:快速傅里叶变换—>O(n∗log(n)O(n∗log(n)计算多项式乘法
FNTT/NTT:快速傅里叶变换的优化版—>优化常数及误差
FWT:快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题
MTT:毛爷爷的FFT—>非常nb/任意模数
FMT 快速莫比乌斯变化—>感谢stump提供
多项式
系数表示法
设A(x)A(x)表示一个n−1n−1次多项式
则A(x)=∑ni=0ai∗xiA(x)=∑i=0nai∗xi
例如:A(3)=2+3∗x+x2A(3)=2+3∗x+x2
利用这种方法计算多项式乘法复杂度为O(n2)O(n2)
(第一个多项式中每个系数都需要与第二个多项式的每个系数相乘)
点值表示法
将nn互不相同的xx带入多项式,会得到nn个不同的取值yy
则该多项式被这nn个点(x1,y1),(x2,y2),…,(xn,yn)(x1,y1),(x2,y2),…,(xn,yn)唯一确定
其中yi=∑n−1j=0aj∗xjiyi=∑j=0n−1aj∗xij
例如:上面的例子用点值表示法可以为(0,2),(1,6),(2,12)(0,2),(1,6),(2,12)
利用这种方法计算多项式乘法的时间复杂度仍然为O(n2)O(n2)
(选点O(n)O(n),每次计算O(n)O(n))
我们可以看到,两种方法的时间复杂度都为O(n2)O(n2),我们考虑对其进行优化
对于第一种方法,由于每个点的系数都是固定的,想要优化比较困难
对于第二种方法,貌似也没有什么好的优化方法,不过当你看完下面的知识,或许就不这么想了
复数
在介绍复数之前,首先介绍一些可能会用到的东西
向量
同时具有大小和方向的量
在几何中通常用带有箭头的线段表示
圆的弧度制
等于半径长的圆弧所对的圆心角叫做1弧度的角,用符号rad表示,读作弧度。用弧度作单位来度量角的制度叫做弧度制
公式:
1∘=π180rad1∘=π180rad
180∘=πrad180∘=πrad
平行四边形定则
(好像画的不是很标准。。)
平行四边形定则:AB+AD=AC
复数
定义
设a,ba,b为实数,i2=−1i2=−1,形如a+bia+bi的数叫复数,其中ii被称为虚数单位,复数域是目前已知最大的域
在复平面中,xx代表实数,yy轴(除原点外的点)代表虚数,从原点(0,0)(0,0)到(a,b)(a,b)的向量表示复数a+bia+bi
模长:从原点(0,0)(0,0)到点(a,b)(a,b)的距离,即a2+b2−−−−−−√a2+b2
幅角:假设以逆时针为正方向,从xx轴正半轴到已知向量的转角的有向角叫做幅角
运算法则
加法:
因为在复平面中,复数可以被表示为向量,因此复数的加法与向量的加法相同,都满足平行四边形定则(就是上面那个)
乘法:
几何定义:复数相乘,模长相乘,幅角相加
代数定义:
单位根
下文中,默认nn为22的正整数次幂
在复平面上,以原点为圆心,11为半径作圆,所得的圆叫单位圆。以圆点为起点,圆的nn等分点为终点,做nn个向量,设幅角为正且最小的向量对应的复数为ωnωn,称为nn次单位根。
根据复数乘法的运算法则,其余n−1n−1个复数为ω2n,ω3n,…,ωnnωn2,ωn3,…,ωnn
注意ω0n=ωnn=1ωn0=ωnn=1(对应复平面上以xx轴为正方向的向量)
那么如何计算它们的值呢?这个问题可以由欧拉公式解决
例如
图中向量ABAB表示的复数为88次单位根
单位根的幅角为周角的1n1n
在代数中,若zn=1zn=1,我们把zz称为nn次单位根
单位根的性质
- ωkn=cosk2πn+isink2πnωnk=cosk2πn+isink2πn(即上面的公式)
- ω2k2n=ωknω2n2k=ωnk
证明:
- ωk+n2n=−ωknωnk+n2=−ωnk
- ω0n=ωnn=1ωn0=ωnn=1
讲了这么多,貌似跟我们的正题没啥关系啊。。
OK!各位坐稳了,前方高能!
快速傅里叶变换
我们前面提到过,一个nn次多项式可以被nn个点唯一确定。
那么我们可以把单位根的00到n−1n−1次幂带入,这样也可以把这个多项式确定出来。但是这样仍然是O(n2)O(n2)的呀!
我们设多项式A(x)A(x)的系数为(ao,a1,a2,…,an−1)(ao,a1,a2,…,an−1)
那么
将其下标按照奇偶性分类
设
那么不难得到
我们将ωkn(k<n2)ωnk(k<n2)代入得
同理,将ωk+n2nωnk+n2代入得
大家有没有发现什么规律?
没错!这两个式子只有一个常数项不同!
那么当我们在枚举第一个式子的时候,我们可以O(1)O(1)的得到第二个式子的值
又因为第一个式子的kk在取遍[0,n2−1][0,n2−1]时,k+n2k+n2取遍了[n2,n−1][n2,n−1]
所以我们将原来的问题缩小了一半!
而缩小后的问题仍然满足原问题的性质,所以我们可以递归的去搞这件事情!
直到多项式仅剩一个常数项,这时候我们直接返回就好啦
时间复杂度:
不难看出FFT是类似于线段树一样的分治算法。
因此它的时间复杂度为O(nlogn)O(nlogn)
快速傅里叶逆变换
不要以为FFT到这里就结束了。
我们上面的讨论是基于点值表示法的。
但是在平常的学习和研究中很少用点值表示法来表示一个多项式。
所以我们要考虑如何把点值表示法转换为系数表示法,这个过程叫做傅里叶逆变换
(y0,y1,y2,…,yn−1)(y0,y1,y2,…,yn−1)为(a0,a1,a2,…,an−1)(a0,a1,a2,…,an−1)的傅里叶变换(即点值表示)
设有另一个向量(c0,c1,c2,…,cn−1)(c0,c1,c2,…,cn−1)满足
即多项式B(x)=y0,y1x,y2x2,…,yn−1xn−1B(x)=y0,y1x,y2x2,…,yn−1xn−1在ω0n,ω−1n,ω−2n,…,ω−(n−1)n−1ωn0,ωn−1,ωn−2,…,ωn−1−(n−1)处的点值表示
emmmm又到推公式时间啦
(c0,c1,c2,…,cn−1)(c0,c1,c2,…,cn−1)满足
设S(x)=∑n−1i=0xiS(x)=∑i=0n−1xi
将ωknωnk代入得
当k!=0k!=0时
等式两边同乘ωknωnk得
两式相减得
观察这个式子,不难看出它分母不为0,但是分子为0
因此,当K!=0K!=0时,S(ωkn)=0S(ωnk)=0
那当k=0k=0时呢?
很显然,S(ω0n)=nS(ωn0)=n
继续考虑刚刚的式子
当j≠kj≠k时,值为00
当j=kj=k时,值为nn
因此,
这样我们就得到点值与系数之间的表示啦
理论总结
至此,FFT的基础理论部分就结束了。
我们来小结一下FFT是怎么成功实现的
首先,人们在用系数表示法研究多项式的时候遇阻
于是开始考虑能否用点值表示法优化这个东西。
然后根据复数的两条性质(这个思维跨度比较大)得到了一种分治算法。
最后又推了一波公式,找到了点值表示法与系数表示法之间转换关系。
emmmm
其实FFT的实现思路大概就是
系数表示法—>点值表示法—>系数表示法
引用一下远航之曲大佬的图
当然,再实现的过程中还有很多技巧
我们根据代码来理解一下
递归实现
递归实现的方法比较简单。
就是按找我们上面说的过程,不断把要求的序列分成两部分,再进行合并
在c++的STL中提供了现成的complex类,但是我不建议大家用,毕竟手写也就那么几行,而且万一某个毒瘤卡STL那岂不是很GG?
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=2*1e6+10;
inline int read()
{
char c=getchar();int x=0,f=1;
while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
return x*f;
}
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