[FFT] 快速傅里叶变换学习笔记

-1、为什么学FFT

退役（很早）之前听说FFT很神（e）奇（xin），Po姐来讲的时候也是膜（sha）了（ye）一（bu）发（dong），于是就放那里了。退役之后有（xian）了（de）时（mei）间（shi），并且在篮球赛之前立了赢一场的flag否则学FFT结果又双叒叕全输了，于是下面就是成果了……网上FFT讲解多得是不看也罢……

0、什么是FFT？为什么FFT？

0.1、为什么FFT？

从一个简单问题说起：大整数乘法。在做VijosP2000的时候，看数据范围O(n2)有点卡常？在做HDOJ1402的时候O(n2)接着卡常？大整数乘法的时间复杂度只能是O(n2)了？不是的。从算法优化上可以实现为O(nlog23)的方法（log23≈1.585，此方法为Karatsuba乘法，详解看这里），还是不太友好。那么就会用到O(nlog2n)的FFT了。

0.2、什么是FFT？

FFT（Fast Fourier Transform），全名快速傅里叶变换，这与傅里叶变换只有半毛钱的关系，傅里叶变换是分析波的成分的方法，通过推广傅里叶变换得到了优化的快速傅里叶变换，为了展示区别和联系，下面是傅里叶变换和傅里叶逆变换的公式：
傅里叶变换：

F(ω)=F[f(t)]=∫−∞∞f(t)e−iωtdt

傅里叶逆变换：

f(t)=F−1[F(ω)]=12π∫−∞∞F(ω)eiωtdω

其中f(t)为以t为周期的周期函数，t满足狄利赫里条件，F(ω)称为f(t)的像函数，f(t)称为F(ω)的原像函数。具体问题请参考《高等数学》，博主还是高中生……（高中高等数学233333）
然后，我们终于可以开始讲FFT了……

1 、（初中知识）多项式

1.1、定义

我们把形如a0+a1x+a2x2+⋯+an−1xn−1的式子叫做多项式，其中a0,a1,⋯,an−1称为多项式的系数，多项式中单项式个数为多项式的项数。因为多项式由单项式组成，规定最高次项的次数叫做多项式的次数，这里A(x)的次数为n−1。x为未知数。
以上内容均为初二（我记得是）的内容，不明白的自觉面壁思过（初中数学老师：这个定义抄10遍给我）……
下面就是没学过的了：
为方便表示，我们记多项式a0+a1x+a2x2+⋯+an−1xn−1为A(x)，即：

A(x)=∑i=0n−1aixi

我们规定：严格大于n的任意整数为这个多项式的次数界，即次数的范围。多项式中任意单项式次数均严格小于次数界。
哦，行……
下面无特殊说明，多项式均用A(x)等类似形式表示。

1.2、多项式的两种表示法

1.2.1、系数表示法

我们知道，向量是个好东西（蛤？），矩阵也是个好东西（蛤？），于是我们把三者结合起来看（蛤玩意？）。
我们将n−1次多项式的系数看成n维向量，即a⃗ =(a0,a1,⋯,an−1)，则称a⃗ 为多项式A(x)的系数表示。
在选修4-2中，我们知道向量与矩阵有着密切联系，这就是为什么我们把多项式，向量和矩阵联系起来。（然而我们并没选4-2蛤蛤蛤）
于是这个n维向量可以表示为一个n维列向量，为：

⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥

1.2.2、点值表示法

我们把多项式A(x)看做一个函数，选取不同的n个值x0,x1,⋯,xn−1代入A(x)中，可以得到A(x)的图像上n个不同的点，那么称点集α={(xi,A(xi))∣i∈[0,n−1],i∈Z}为多项式A(x)的点值表示。
到这里，是不是差不多忘了要干什么了……

1.3、多项式的运算

1.3.1、多项式求值

这个比较简单，拿题说。这是我还没出也不想出的一道题，拿出来娱乐一下……
题目描述 Description

都知道FZ酱数学不好，数学老师很着急，于是让她求个多项式的值。但是数学老师一不小心数据就出大了却并没有注意到这一点。老师共出了T道题，FZ酱拿到了这T道题就崩溃了，请帮她尽快完成作业。

输入 Input

第一行一个整数T，表示数据组数。
对于每组数据：
第一行有一个数n表示这个多项式的次数界；
第二行有n个整数a0,a1,⋯ ,an−1，表示多项式∑i=0n−1aixi；
第三行，为未知数的值x。

输出 Output

对于每组数据，输出一行，为答案。
样例输入 Sample Input
2
2
3 2
100
3
1 2 3
-9

样例输出 Sample Output

203
226

样例解释 Explanation

对于第一个多项式为A(x)=2x+3，代入后得A(100)=2×100+3=203；
对于第二个多项式为A(x)=3x2+2x+1，代入后得A(−9)=226。

限制 Limits

对于50%的数据，答案在[0,2128−1]范围内；
对于100%的数据，T≤2，n≤105，保证输入均为绝对值不大于231−1的整数并且答案小于3∗106位。

那么很好我们必须写高精度了……
2128−1≈3.4∗1038，这里的50%数据是给不想写高精度的童鞋写__int128准备的（__int128的范围是[0,2128−1]）。于是我们就有这样三种的方法。
记数字长度为p。

1.3.1.1、方法一：O(pn2) ——>50pts

暴力不会吗？
核心代码如下：

1.3.1.2、方法二：O(p2nlog2n)——>[50,80]pts

这里挨个乘太慢了，快速幂处理会快一些。
（其实不一定会快多少）
核心代码如下：

1.3.1.3、方法三：O(p2n)——>100pts

学过必修三吗？
学过必修三还不会？
拿出数学必修三翻到37页你看到了什么？
秦九韶算法可以大量减少乘法和加法次数，并且把运算量简化为O(n)的。
核心代码如下：

当然，还可以处理前缀积，也是O(n)的。
核心代码如下：

不过没上面的优雅不是吗……（呵呵）
回到正题上，我们要谈的其实是……

1.3.2、多项式的加法和乘法

1.3.2.1、多项式加法

没啥难度，直接加即可：

C(x)=A(x)+B(x)=∑i=0n−1aixi+∑j=0n−1bjxj=∑i=0n−1(ai+bi)xi=∑i=0n−1cixi

这个操作可以O(n)完成。

1.3.2.2、多项式乘法

这就是个开括号的问题……
两个次数界为na,nb的多项式A(x),B(x)相乘，得到的多项式C(x)的次数界为na+nb−1。
我们把不足的次数用系数0补充，变成两个齐次多项式。
考虑怎样的两项乘完后，结果的幂指数为同一个。即aixi×bjxj=ckxk，i,j怎样取值才能使k为定值。显然k=i+j。
那么乘完之后，C(x)就可以知道了。

C(x)=∑i=02n−2(∑j=0iajbi−j)xi=∑i=02n−2cixi

我们把中间的那个加法，即ci=∑j=0iajbi−j叫做卷积。
这样，朴素算法时间复杂度为O(n2)，显然不是太好。
但是点值表达式的表示就很好了，直接乘起来就好了嘛！即：
设A(x)的点值表示为α={(xi,A(xi))∣i∈[0,n−1],i∈Z}，设B(x)的点值表示为β={(xi,B(xi))∣i∈[0,n−1],i∈Z}，然后就能得到C(x)的点值表示为：γ={(xi,A(xi)B(xi))∣i∈[0,n−1],i∈Z}，但是这里只是两多项式的积。并不是真正的点值表示，不过也够用了。

1.4、插值-点值转换

就要接近重点了同志们加油！
先定义两个运算：
1、定义一种运算，可以将系数表示转化为点值表示，记为DFT(A(x))=α。
2、定义其逆运算，即从一个多项式的点值表示确定其系数表示为插值。记为DFT−1(α)=IDFT(α)=A(x)
第一个运算在xi确定时唯一确定，那么第二个在xi确定时可不可以唯一确定A(x)呢？
答案是肯定的。我们把这些点代入，可以得出一个n元一次方程组，当xi均不相同时，可以唯一确定a0,a1,⋯,an。（如果有一个xi相同立刻变成不定方程）
当然可以拿矩阵证明网上多得是。
那么就很尴尬，DFT(A(x))的复杂度是O(n2)的，它的逆运算是O(n3)的（高斯消元）。还不如不优化了。
所以我们还需要一些姿势。

2、（高中知识）复数

一些奇奇怪怪的问题无法解决，可以把它从实数域扩展到复数域上解决，加入了复数这种类似向量的东西，问题就会变得简单。
噫，向量？点值表示不就是个向量吗！

2.1、定义

一堆定义就不说了，选修2-2内容，这里简单的强调一些运算：
i2=−1
设z1,z2∈C，z1=a+bi,z2=c+di，则：
z1±z2=(a±c)+(b±d)i
z1×z2=(ac−bd)+(ad+bc)i
z1z2=ac+bdc2+d2+bc−adc2+d2i
设z∈C，z=a+bi，则z的共轭复数为z¯=a−bi
这是学过的，下面就说说没学过的。

2.2、单位根

根？方程的根？什么方程的根？
定义：xn=1的所有根叫做n次单位根。
当我们把问题引入复数域，x就不能简单计算为x=1√n=1了，而是在复平面内解决问题。

2.3、欧拉公式

一个很6的公式。内容是：

eix=cosx+isinx

式子很简单，证明如下：
由Taylor展开：

excosxsinx=1+x+x22!+⋯+xnn!+⋯=1−x22!+x44!−x66!+⋯=x−x33!+x55!−x77!+⋯

将Taylor展开推广到复数域上，有复数运算±i=±i,(±i)2=−1,(±i)3=∓i,(±i)4=1，可得：

eix=1+ix−x22!−ix33!+x44!+ix55!+⋯=(1−x22!+x44!−⋯)+i(x−x33!+x55!−⋯)=cosx+isinx

证毕。
将x换成π，就出现了eπi=−1，即eπi+1=0，被誉为上帝创造的公式。
将x换成2kπ  (k∈Z)，就出现了e2kπi=1。
我们不是要解xn=1吗！
那么x=e2kπin  (k∈Z)，记这些单位根为ωkn=e2kπin。
下面就是一些数论问题了……
记k=nr+p，其中p∈[0,n−1]，则ωkn=e2kπin

ωkn=e2kπin=e2(nr+p)πin=e2nrπin+2pπin=e2nrπin×e2pπin=e2πir×e2pπin=1r×e2pπin=e2pπin

这不一样吗……
所以可以证明k∈[0,n−1]。也就是n次单位根有n个，不难发现，这n个单位根就是单位圆上的n个n等分点。
于是记这些根为ω0n,ω1n,⋯,ωn−1n，这些根各不相同，并且在乘法意义下构成一个群。即：ωinωjn=ωi+jn=ω(i+j)modnn，且ω−1n=ωn−1n，证明如下：

ω−1n=1×ω−1n=e2πi×e−2πin=e2πi−−2πin=e2(n−1)πin=ωn−1n

问题得证。

2.4、三个引理

2.4.1、消去引理

内容：当d≠0时，ωdkdn=ωkn。
证明：

ωdkdn=e2dkπidn=e2kπin=ωkn

很简单……也很智障

2.4.2、折半引理

内容：如果n>0且n为偶数，(ωkn)2=ωkn2
证明：

(ωkn)2=(e2kπin)2=e2kπin2=ωkn2

还行吧……可以看出，n个n次单位根平方的集合就是n2个n2次单位根的集合，并且每个元素出现两次。可以通过单位根的对称性来理解。
更一般的结论：
ωk+n2n=ωn2nωkn=−ωkn
ωn2n=eπi=−1

2.4.3、求和引理

内容：对于∀k∈(0,n)，有∑j=0n−1(ωkn)j=0。
证明：

∑j=0n−1(ωkn)j=1−(ωkn)n1−ωkn（等比数列求和）=1−(ωnn)k1−ωkn=0

注意：ωnn=1，化简省略了部分很多步骤。
现在全都就绪了，终于要到目标了……

3、离散傅里叶变换（DFT）

DFT，即Discrete Fourier Transform，离散傅里叶变换，就是我们之前定义的那个，现在有了复数，我们就可以在O(nlog2n)的时间内做完DFT运算。运用到了Cooley-Tukey算法。
我们将多项式按照指数的奇偶分类，记原多项式为A(x)，那么构造两个多项式：

A0(x)A1(x)=a0+a2x+a4x2+⋯+an−2xn2−1=∑i=0n2−1a2ixi=a1+a3x+a5x3+⋯+an−1xn2−1=∑i=0n2−1a2i+1xi

通过观察可知，A(x)=A0(x2)+xA1(x2)，根据折半引理，我们考虑将n次单位根作为x代入计算。又因为折半引理，我们即可分治操作。
那么，对于k<n2，有：

A(ωkn)A(ωk+n2n)=A0((ωkn)2)+ωknA1((ωkn)2)=A0(ωkn2)+ωknA1(ωkn2)=∑i=0n2−1a2iωkin2+ωkn∑i=0n2−1a2i+1ωkin2=A0((ωkn)2)+ωk+n2nA1((ωkn)2)=A0(ωkn2)−ωknA1(ωkn2)=∑i=0n2−1a2iωkin2−ωkn∑i=0n2−1a2i+1ωkin2

这样，由奇偶性使得需要代入的值范围减半，递归做下去就行了，时间复杂度为：

T(n)=2T(n2)+O(n)=O(nlog2n)

因为永远是奇偶合并，如果想要更方便地合并必然n取2的整数次幂，否则合并到一定时候左右不一定都是相同次数（可以考虑一个完全二叉树）。所以做的时候取2的整数次幂作为n。

4、逆离散傅里叶变换（IDFT）

IDFT，即Inverse Discrete Fourier Transform，逆离散傅里叶变换，还是我们之前定义的那个。我们在O(nlog2n)的时间里搞定了DFT，可是IDFT还是个O(n3)的，所以我们继续观（tui）察（da）性（shi）质（zi）。
我们说过，IDFT相当于解个n元一次方程组，这个方程组写成矩阵形式是这样的：

⎡⎣⎢⎢⎢⎢⎢⎢(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1(ω0n)2(ω1n)2⋮(ωn−1n)2⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1⎤⎦⎥⎥⎥⎥⎥⎥×⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢A(ω0n)A(ω1n)A(ω2n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥⎥⎥⎥

这同时也是做DFT时我们乘的矩阵。只不过这次ai变成未知数了。
记上面的系数矩阵为D，再考虑以下矩阵：

⎡⎣⎢⎢⎢⎢⎢⎢(ω0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω0n)1(ω−1n)1⋮(ω−(n−1)n)1(ω0n)2(ω−1n)2⋮(ω−(n−1)n)2⋯⋯⋱⋯(ω0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎥

记此矩阵为W，则记E=D\times WE=D×W。可知：
\begin{equation}\begin{aligned}E_{ij}&=\sum\limits_{k=0}^{n-1}D_{ik}W_{kj}\\&=\sum\limits_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}\\&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\end{aligned}\end{equation}

Eij=∑k=0n−1DikWkj=∑k=0n−1ω−iknωkjn=∑k=0n−1ωk(j−i)n

当i=ji=j时：
\begin{equation}\begin{aligned}E_{ij}&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\\&=\sum\limits_{k=0}^{n-1}\omega_n^0=n\\\end{aligned}\end{equation}

Eij=∑k=0n−1ωk(j−i)n=∑k=0n−1ω0n=n

当i\not =ji≠j时，由求和引理：
\begin{equation}\begin{aligned}E_{ij}&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\\&=\sum\limits_{k=0}^{n-1}(\omega_n^{j-i})^k=0\\\end{aligned}\end{equation}

Eij=∑k=0n−1ωk(j−i)n=∑k=0n−1(ωj−in)k=0

记单位矩阵为II，所以I=\frac{1}{n}EI=1nE，即\frac{1}{n}D=W^{-1}1nD=W−1。
将由方程组推得那个矩阵左乘1nD，可得：

⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢(ω0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω0n)1(ω−1n)1⋮(ω−(n−1)n)1(ω0n)2(ω−1n)2⋮(ω−(n−1)n)2⋯⋯⋱⋯(ω0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎥×⎡⎣⎢⎢⎢⎢⎢⎢⎢A(ω0n)A(ω1n)A(ω2n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥⎥⎥⎥

这和DFT的操作很像，所以用DFT的过程实现IDFT即可。
负号？把ω1n改成ω−1n即可；
1n？做完之后再除个n即可。现在IDFT也是O(nlog2n)的了。
至此，所有的时间都是O(nlog2n)的了！即FFT的复杂度为O(nlog2n)！

5、实现问题

说是一回事，写是另一回事。
奇偶分组这个玩意就很坑，容易发现，设ri为将i的高位翻转到低位得到的数字（如0001就变成1000），那么初始位置为i的数就会被移动到ri，所以我们可以预先完成移动，然后合并计算。时间就会加快很多了。
不过FFT本身常数巨大，主要因为double类型的计算时间，所以常数问题一定要注意。
我们看到，DFT做完后，得到的是将n次单位根代入后表达式的点值表示，IDFT做完后，得到的是两多项式相乘的系数表示。
由于FFT涉及复数和除法，容易出现精度问题，所以考虑一种整数在模的前提下的FFT，就出现了FNT，也叫NTT，即快速数论变换（Fast Number-theory Transform或Number Theory Transform）。这里有学习笔记……
如果我们计算卷积的时候a,b的下标不是相加为定值，而是一或二元逻辑运算（与，或，异或等等）后为定值，则可应用快速沃尔什变换（Fast Walsh Hadamard Transform）。学习笔记在这里
FFT真是大毒瘤……
FFT和FWT都在图像处理，音频处理方面有较大应用。
FFT的板子在这里，以HDOJ1402为例……Code
VijosP2000那道题除了写FFT还得压位……或者用一种神奇的O(n2)乘法：Code
终于把我的flag填完了我好开心啊！！！