【OI向】快速傅里叶变换(Fast Fourier Transform)

转存于本人博客园地址

FFT的作用

在学习一项算法之前,我们总该关心这个算法究竟是为了干什么。

（以下应用只针对OI）

一句话:求多项式乘法(当然它的实际用处很多)

设多项式

A(x)=a0+a1x+a2x2+…+anxnA(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^nA(x)=a0+a1x+a2x2+…+anxn

B(x)=b0+b1x+b2x2+…+bmxmB(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^mB(x)=b0+b1x+b2x2+…+bmxm

我们想求 F(x)=A(x)B(x)=∑i=0n∑j=0maibjxi+jF(x)=A(x)B(x)=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx^{i+j}F(x)=A(x)B(x)=i=0∑nj=0∑maibjxi+j

肉眼可见的复杂度O(nm)O(nm)O(nm)

而使用一般的FFT则复杂度能够在O(nlogn)O(nlogn)O(nlogn) （实则带比较大的常数,但数据范围大时肯定比O(nm)O(nm)O(nm)好）

写这篇博客是为了让自己和看到的人能够梳理一遍FFT这个算法。

以下FFT的实现使用比较常见的复数单位根做法。

前置知识

首先要知道复数是什么,高中选修即有,网上资料也很多,不赘述,这里只讲FFT直接相关。

1 复数单位根

定义: zn=1z^n=1zn=1 那么 zzz 就是n次的单位根

对于一个nnn次单位根,就是将复数单位圆等分成nnn份,每个数 wnk(k∈[0,n−1))w_n^k (k\in[0,n-1) )wnk(k∈[0,n−1)) 就被称作n次的第k个单位根

如图

图中表示了所有8次单位根,其中w80=1,w81=1+i,w82=i,w83=−1+i…w_8^0=1,w_8^1=1+i,w_8^2=i,w_8^3=-1+i\ldotsw80=1,w81=1+i,w82=i,w83=−1+i… 所谓单位根就是这些复数,它们的八次方均为1(几次单位根就是几次方,见单位根的定义)。

虽然单位根的次数nnn在定义上可以取全体正整数,但是我们在FFT中只需要考虑n为2的整数次幂,以方便我们利用性质。

单位根有一些性质,我们先列出来看看（下边的n均为2的整次幂!）

0 单位根乘法

wni×wnj=wni+jw_n^i\times w_n^j=w_n^{i+j}wni×wnj=wni+j,这个性质实际上不能叫性质,从复数乘法,辐角的概念出发即能够得证。

1 消去引理

w2n2k=wnkw_{2n}^{2k}=w_n^kw2n2k=wnk

理解,可以看到 2k2k2k 和 kkk 在 2n2n2n 和 nnn 中占的比例是相同的,联想单位圆即可(w82w_8^2w82显然等于w41w_4^1w41,均为1)

2 折半引理

我看了好几篇博客,包括知乎,对折半引理的说法是: 对于 nnn 等于 222 的整次幂, wnkw_n^kwnk经过平方可以得到所有wn2kw_{\frac{n}{2}}^kw2nk(单纯来说只需要nnn是二的倍数就好)

两个小性质都可以证明

wn2k=wn2kw_n^{2k}=w_{\frac{n}{2}}^kwn2k=w2nk

由消去引理,这个小性质是显然的。从这个性质出发,如果 kkk 取遍[0,n−1][0,n-1][0,n−1],那么首先右边是全了,对于左边,2k∈[0,2n−2]2k\in[0,2n-2]2k∈[0,2n−2],经历了一次循环, (按道理这里需要指出一个循环性质,但我觉得没必要) wnk+n=wnk×wnnw_n^{k+n}=w_n^{k}\times w_n^{n}wnk+n=wnk×wnn, 而wnnw_n^nwnn即wn0=1w_n^0=1wn0=1,所以 2k2k2k 过了n−1n-1n−1之后会经过n2\frac{n}{2}2n个相同的点

wnk=−wnk+n2w_n^k=-w_n^{k+\frac{n}{2}}wnk=−wnk+2n

旋转180∘180^\circ180∘ 事情变得显然;从复数乘法的角度来说也可以：wnk+n2=wnk×wnn2w_n^{k+\frac{n}{2}}=w_n^k\times w_n^{\frac{n}{2}}wnk+2n=wnk×wn2n,而wnn2w_n^{\frac{n}{2}}wn2n恰好为−1-1−1。

那么这个性质告诉我们前一半k∈[0,n2−1]k\in[0,\frac{n}{2}-1]k∈[0,2n−1]和后一半 k∈[n2,n−1]k\in[\frac{n}{2},n-1]k∈[2n,n−1]的平方是一样的,所以同理,两半的平方值是一样的,也就是两部分平方只确定了n2\frac{n}{2}2n个值,那么利用前边第一个小性质,我们不需要循环性质就能够证明折半引理。

折半引理告诉我们,如果知道了n的所有单位根,那么n2\frac{n}{2}2n的单位根我们可以知道,反之,我们利用第二个性质也可以从n2\frac{n}{2}2n的所有单位根,加一个负号,就能推出所有的 nnn 次单位根。

3 求和引理

∑i=0n−1wnk=0\sum\limits_{i=0}^{n-1}w_n^k=0i=0∑n−1wnk=0

从折半引理看,我们把 wnkw_n^kwnk 分成两部分wna(a∈[0,n2−1]),wnb(b∈[n2,n−1])w_n^a (a\in[0,\frac{n}{2}-1] ),w_n^b(b\in[\frac{n}{2},n-1])wna(a∈[0,2n−1]),wnb(b∈[2n,n−1]),也就是圆的上半部分和下半部分,对应有每个wna+wnb=0w_n^a+w_n^b=0wna+wnb=0，实际上是wna+wna+n2=0w_n^a+w_n^{a+\frac{n}{2}}=0wna+wna+2n=0,所以两部分的累加和就是0(我们上边已经假设n是2的整次幂了,所以一定是相同大小的两部分)

单位根的表示

首先简单说一下辐角,在复数平面上,从实数轴正半轴逆时针旋转到复数zzz到原点的连线,中间经过的角度就叫做辐角(大小在2π2\pi2π之内)

从圆的知识类比一下,我们可以知道

在叙述 nnn 次单位根这件事情的时候,我们是将复数单位圆 nnn 等分之后得到单位根,这样来看的话,显然wn1w_n^1wn1的辐角是 θ=2πn\theta=\frac{2\pi}{n}θ=n2π, 那么wnk=(wn1)kw_n^k=(w_n^1)^kwnk=(wn1)k ,那么wnkw_n^kwnk的辐角就是 kθk\thetakθ,我们知道在普通二维平面上单位圆上的点表示为(cosθ,sinθ)(cos\theta,sin\theta)(cosθ,sinθ),类比到复平面,就会得到 wnk=cosθ+sinθi(θ=2πn)w_n^k=cos\theta+sin\theta\ i \ \ \ \ (\theta=\frac{2\pi}n)wnk=cosθ+sinθ i (θ=n2π)

至此,我们已经了解单位根的性质及其表示,然而我们还不知道为什么要介绍这么一个东西。下面我们会了解的。

2 多项式的点值表示法及其点值表示法的乘法

一开始,我们做介绍时设了两个多项式

A(x)=a0+a1x+a2x2+…+anxnA(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^nA(x)=a0+a1x+a2x2+…+anxn

B(x)=b0+b1x+b2x2+…+bmxmB(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^mB(x)=b0+b1x+b2x2+…+bmxm

仔细看的话,我们发现这两个多项式像什么?没错，像函数,实际上它完全可以当成一个函数表达式。

即y=A(x)y=A(x)y=A(x)

我们知道,求解一次函数只需要知道两个点的坐标,求解二次函数只需要知道三个点的坐标,那么我们可以换个说法,两个点确定了一个一次函数,三个点确定了一个二次函数,同样的,我们可以说 n+1n+1n+1 个点确定了一个 nnn 次函数,确定的方法很简单,我们带入联立即可(想想我们是怎么通过两个点确定一条直线的)。

现在我们把名词换一下,”函数“变成“多项式”,就可以改成说 n+1n+1n+1 个不同点确定一个 nnn 次多项式。

即点集((x0,y0),(x1,y1),…,(xn,yn))((x_0,y_0),(x_1,y_1),\ldots ,(x_n,y_n))((x0,y0),(x1,y1),…,(xn,yn)) 能够确定一个 nnn 次多项式。

代入y=A(x)y=A(x)y=A(x)看的形象一点，我们解下边的方程组(点是我们已知且在函数上的),就能够求出来所有的系数 aaa (就是和求函数解析式一样的概念),如果用高斯消元来解这样一个方程组,复杂度是 O(n3)O(n^3)O(n3)的，这个过程被称作离散傅里叶逆/反变换(IDFT) 。
{y0=a0+a1x0+a2x02+…+anx0ny1=a0+a1x1+a2x12+…+anx1n⋮yn=a0+a1xn+a2xn2+…+anxnn\begin{cases} y_0=a_0+a_1x_0+a_2x_0^2+\ldots+a_nx_0^n \\ y_1=a_0+a_1x_1+a_2x_1^2+\ldots+a_nx_1^n\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots \\ y_n=a_0+a_1x_n+a_2x_n^2+\ldots+a_nx_n^n \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧y0=a0+a1x0+a2x02+…+anx0ny1=a0+a1x1+a2x12+…+anx1n ⋮yn=a0+a1xn+a2xn2+…+anxnn
同时,如果我们知道一个多项式的系数表达式,我们只需要随便找n+1n+1n+1个点 xxx 代入到A(x)A(x)A(x)中,就可以得到所有的yyy。

这样我们就可以将一个多项式从系数表达式 O(n2)O(n^2)O(n2) 的转换成点值表达式。

把多项式的系数表示换成点值表示的变换被称作离散傅里叶变换(DFT)。(不小心先介绍了逆变换)

点值表达式之间想要确定新的点值表达式,只用将y相乘即可(xxx对应相等),比如说

A=((x0,y0),(x1,y1)…,(xn,yn))A=((x_0,y_0),(x_1,y_1)\ldots,(x_n,y_n))A=((x0,y0),(x1,y1)…,(xn,yn))

B=((x0,y1′),(x1,y2′)…,(xn,yn′))B=((x_0,y_1'),(x_1,y_2')\ldots,(x_n,y_n') )B=((x0,y1′),(x1,y2′)…,(xn,yn′))

F=AB=((x0,y1y1′),(x1,y2y2′),(x2,y2y2′)…,(xn,ynyn′))F=AB=((x_0,y_1y_1'),(x_1,y_2y_2'),(x_2,y_2y_2')\ldots,(x_n,y_ny_n') )F=AB=((x0,y1y1′),(x1,y2y2′),(x2,y2y2′)…,(xn,ynyn′))

因为我们的每个yiyi′y_iy_i'yiyi′都是一次多项式乘法的结果,按照定义,F(xk)=A(xk)B(xk)=ykyk′=∑i=0n∑j=0maibjxki+jF(x_k)=A(x_k)B(x_k)=y_ky_k'=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx_k^{i+j}F(xk)=A(xk)B(xk)=ykyk′=i=0∑nj=0∑maibjxki+j我们的多项式乘法本身是只关心系数的,和 xxx 无关(就像函数表达式y=f(x)y=f(x)y=f(x)只和各项系数有关而与具体的点无关一样),而根据我们点值表示法的定义 yk=F(xk)y_k=F(x_k)yk=F(xk)，所以(x0,y1y1′)(x_0,y_1y_1')(x0,y1y1′)就是F(x)F(x)F(x)上的一个点 (xk,F(xk))(x_k,F(x_k))(xk,F(xk))

一个小问题一开始我们知道,多项式A(x)B(x)A(x)B(x)A(x)B(x)的结果理应得到最高次项为n+mn+mn+m的F(x)F(x)F(x),而A，B的直接转化的点值表示法分别只有nnn个点和mmm个点，最后要真正确定F(x)F(x)F(x),需要 n+m+1n+m+1n+m+1个点怎么办？好办,我们一开始把A和B多算一些y,补若干个点补到n+m+1n+m+1n+m+1个点,多一些点并不会影响什么,而这样就能够创造出n+m+1n+m+1n+m+1个能够确定F(x)F(x)F(x)的点集了。

AB的运算得到的新点集,就能够确定新的多项式 F(x)F(x)F(x) 。而在已知 yiy_iyi和 yi′y_i'yi′ 的情况下，运算AB需要的复杂度是 O(n)O(n)O(n)的！

快速傅里叶变换(FFT)

从前置知识点值表示法的内容中,我们可以得到一种新的求多项式 F(x)=A(x)B(x)F(x)=A(x)B(x)F(x)=A(x)B(x) 的方法

DFT,将A(x)A(x)A(x) B(x)B(x)B(x) 通过离散傅里叶变换 (DFT)转化为点值表示法,为了能够真正确定F(x)F(x)F(x)，A和B都要转化为n+m+1n+m+1n+m+1个点以确定最高次项为n+mn+mn+m的F(x)F(x)F(x),复杂度 O(n2)O(n^2)O(n2)
将A(x)A(x)A(x)的所有点值和B(x)B(x)B(x)的所有点值对应相乘,得到F(x)F(x)F(x)的点值表示法。复杂度 O(n)O(n)O(n)
IDFT,以所有 xxx 的各个次项做系数矩阵(以xi1,xi2…x_i^1,x_i^2\ldotsxi1,xi2…做系数矩阵),进行高斯消元,得到F(x)F(x)F(x)的所有系数。复杂度 O(n3)O(n^3)O(n3)

非常漂亮啊!!我们一通操作,直接得到一个O(n³)的做法能够进行多项式乘法！！秒！

我们进行这么多操作,研究了点值表示法能够对应多项式的系数表示法的事实,当然不是为了花里胡哨得到一个O(n3)O(n^3)O(n3)的做法

这个时候,我们就要考虑优化了,优化什么呢？我们不忍心点值表示运算明明已经达到了极优秀的O(n)O(n)O(n),而DFT和IDFT却严重拖慢运行。

这个时候,我们就要用到复数单位根了。

平时我们做函数之类的事情,通常规定 x∈Rx\in Rx∈R 实际上我们可以扩展数域, 使得 x∈Cx\in Cx∈C 即让xxx从实数扩展到整个复数域,作为一个横坐标来使用。然后我们考虑将DFT的过程中任选的xxx,变成有规律的复数单位根,为了使得这些 xxx 有规律我们还规定,对于一个n−1n-1n−1次多项式,我们向多项式 A(x)A(x)A(x) 和B(x)B(x)B(x) 代入所有的单位根,即代入所有的 wnkw_n^kwnk ,而所有的单位根是不同的复数,所以点集(wnk,A(wnk))(w_n^k,A(w_n^k))(wnk,A(wnk)) 就能够作为A(x)A(x)A(x)的点值表示法,B(x)B(x)B(x)同理。

接下来我们将 nnn 次单位根代入到我们的多项式中,为了使用我们上边指出过的复数单位根的性质,以及为了能够唯一确定我们的目标多项式,我们选择的 nnn 应当是2的整次幂,并且它要大于等于n+m+1。

代入之后 A(x)=((wn0,A(wn0),),(wn1,A(wn1))…,(wnn−1,A(wnn−1))A(x)=((w_n^0,A(w_n^0),),(w_n^1,A(w_n^1))\ldots,(w_n^{n-1},A(w_n^{n-1}))A(x)=((wn0,A(wn0),),(wn1,A(wn1))…,(wnn−1,A(wnn−1)),考虑如何快速求出每个A(wnk)A(w_n^k)A(wnk)，我们从消去引理知道,只要知道偶数项的wnk=wn2k2w_n^k=w_\frac{n}2^{\frac{k}2}wnk=w2n2k,这提示我们需要将原来的表达式进行一波奇偶分组,从系数表达式入手,现在A(x)A(x)A(x)是一个n−1n-1n−1次多项式。
A(x)=(a0+a2x2+a4x4+…+an−2xn−2)+(a1+a3x3+…+an−1xn−1)=(a0+a2x2+a4x4+…+an−2xn−2)+x(a1+a3x2+a5x4+…+an−1xn−2)A(x)=(a_0+a_2x^2+a_4x^4+\ldots+a_{n-2}x^{n-2})+(a_1+a_3x^3+\ldots+a_{n-1}x^{n-1})\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =(a_0+a_2x^2+a_4x^4+\ldots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+\ldots+a_{n-1}x^{n-2})\\ A(x)=(a0+a2x2+a4x4+…+an−2xn−2)+(a1+a3x3+…+an−1xn−1) =(a0+a2x2+a4x4+…+an−2xn−2)+x(a1+a3x2+a5x4+…+an−1xn−2)
设A1(x)=(a0+a2x+a4x2+…+an−2xn2−1),A2(x)=(a1+a3x+…+an−1xn2−1)A_1(x)=(a_0+a_2x+a_4x^2+\ldots+a_{n-2}x^{\frac{n}2-1}),A_2(x)=(a_1+a_3x+\ldots+a_{n-1}x^{\frac{n}2-1})A1(x)=(a0+a2x+a4x2+…+an−2x2n−1),A2(x)=(a1+a3x+…+an−1x2n−1)

我们会得到A(x)=A1(x2)+xA2(x2)A(x)=A_1(x^2)+xA_2(x^2)A(x)=A1(x2)+xA2(x2)

将所有xxx 写成单位根的样子
A(wnk)=a0+a1wnk+a2wn2k+a33k+…+an−1wn(n−1)k=A1(wn2k)+wnkA2(wn2k)=(a0+a2wn2k+a4wn4k+…+an−2wn(n−2)k)+wnk(a1+a3wn2k+…+an−1(n−2)k)=(a0+a2wn2k+a4wn22k+…+an−2wn2(n2−1)k)+wnk(a1+a3wn2k+…+an−1(n2−1)k)=A1(wn2k)+wnkA2(wn2k)\begin{aligned} A(w_n^k)& =a_0+a_1w_n^k+a_2w_n^{2k}+a_3^{3k}+\ldots+a_{n-1}w_n^{(n-1)k} \\&=A_1(w_n^{2k} )+w_n^kA_2(w_n^{2k} )\\ &=(a_0+a_2w_n^{2k}+a_4w_n^{4k}+\ldots+a_{n-2}w_n^{(n-2)k})+w_n^k(a_1+a_3w_n^{2k}+\ldots+a_{n-1}^{(n-2)k})\\&=(a_0+a_2w_\frac{n}2^k+a_4w_\frac{n}2^{2k}+\ldots+a_{n-2}w_\frac{n}2^{(\frac{n}2-1)k})+w_n^k(a_1+a_3w_\frac{n}2^{k} +\ldots+a_{n-1}^{(\frac{n}2-1)k})\\ &=A_1(w_\frac{n}2^k)+w_n^kA_2(w_\frac{n}2^k) \end{aligned} A(wnk)=a0+a1wnk+a2wn2k+a33k+…+an−1wn(n−1)k=A1(wn2k)+wnkA2(wn2k)=(a0+a2wn2k+a4wn4k+…+an−2wn(n−2)k)+wnk(a1+a3wn2k+…+an−1(n−2)k)=(a0+a2w2nk+a4w2n2k+…+an−2w2n(2n−1)k)+wnk(a1+a3w2nk+…+an−1(2n−1)k)=A1(w2nk)+wnkA2(w2nk)
对于k∈[0,n2−1]k\in[0,\frac{n}2-1]k∈[0,2n−1]，A(wnk)=A1(wn2k)+wnkA2(wn2k)A(w_n^k)=A_1(w_\frac{n}2^k)+w_n^kA_2(w_\frac{n}2^k)A(wnk)=A1(w2nk)+wnkA2(w2nk)

对于k∈[n2,n−1]k\in[\frac{n}2,n-1]k∈[2n,n−1], 从wnk=−wnk+n2,wnk=wn2k2w_n^k=-w_n^{k+\frac{n}{2}},w_n^k=w_\frac{n}2^\frac{k}2wnk=−wnk+2n,wnk=w2n2k由于kkk已经大于n2\frac{n}22n,所以这里所有的A1(wn2k)=A1(wn2k−n2)A_1(w_\frac{n}2^k)=A_1(w_\frac{n}2^{k-\frac{n}2})A1(w2nk)=A1(w2nk−2n),A2A_2A2同理

所以A(wnk)=A1(wn2k−n2)−wnk−n2A2(wn2k−n2)A(w_n^k)=A_1(w_{\frac{n}2}^{k-\frac{n}2})-w_n^{k-\frac{n}2}A_2(w_{\frac{n}2}^{k-\frac{n}2})A(wnk)=A1(w2nk−2n)−wnk−2nA2(w2nk−2n) 也就是我们可以直接引用在前一半已经算出来过的A1A2A_1A_2A1A2,来计算出后一半的所有A(wnk)A(w_n^k)A(wnk) ,这样我们把计算的规模缩小了一半,同理,我们要知晓前一半A1A_1A1和A2A_2A2的值,也可以从一半的前一半推之,即递归便可以求解,最后需要进行log层计算,然后只需要在递归的底层对一个点进行 O(n)O(n)O(n) 的多项式运算然后每层O(n)O(n)O(n)的向上合并,最终的复杂度是 T(n)=2T(n2)+O(n)T(n)=2T(\frac{n}2)+O(n)T(n)=2T(2n)+O(n) 其最终复杂度是O(nlogn)O(nlogn)O(nlogn)的。这样我们通过将多项式拆分成两部分然后进行相似性的递归,使得我们的DFT能够在O(nlogn)O(nlogn)O(nlogn) 的复杂度下得到多项式的点值表达式,这个过程就被称作FFT。

快速傅里叶逆变换（IFFT）

如何快速的把得到的点值表示法的多项式变成系数表示法,先放一个结论吧
有多项式F(x)=∑i=0n−1fixi那么多项式的点值表示法就是((wn0,F(wn0)),(wn1,F(wn1)),…,(wnn−1,F(wnn−1)))则每个多项式的系数fk=1n∑i=0n−1yiwn−ki有多项式\ F(x)=\sum\limits_{i=0}^{n-1}f_ix^i\\ 那么多项式的点值表示法就是((w_n^0,F(w_n^0)),(w_n^1,F(w_n^1)),\ldots,(w_{n}^{n-1},F(w_n^{n-1})) ) \\则每个多项式的系数 f_k=\frac{1}n\sum\limits_{i=0}^{n-1}y_iw_n^{-ki} 有多项式 F(x)=i=0∑n−1fixi那么多项式的点值表示法就是((wn0,F(wn0)),(wn1,F(wn1)),…,(wnn−1,F(wnn−1)))则每个多项式的系数fk=n1i=0∑n−1yiwn−ki
也就是说,如果我们知道了一个多项式,可以用FFT化作点值表示,而如果我们知道了所有点值反过来求所有系数,我们发现求每个系数的形式和求每个点值的形式是相同的不同的是多了一个系数1n\frac{1}nn1以及在www上标多了一个负号,而这是不影响单位根的三项引理的,也就说,这个反向的过程依然可以用FFT解决。如此一来,我们就可以同样的在O(nlogn)O(nlogn)O(nlogn) 的复杂将点值表示转化为系数表示。

关于IFFT的证明:

我们把每个点展开成多项式,即把(wnk,A(wnk))(w_n^k,A(w_n^k) )(wnk,A(wnk))还原成多项式的模样然后以多项式的系数作为列向量相乘,就能够得到所有的点值yyy
[1,(wn0)1,(wn0)2,…,(wn0)n−11,(wn1)1,(wn1)2,…,(wn1)n−1⋮1,(wnn−1)1,(wnn−1)2,…,(wnn−1)n−1)][a0a1a2⋮an−1]=[y0y1y2⋮yn−1]\begin{bmatrix} &1,&(w_n^0)^1,&(w_n^0)^2,&\ldots,&(w_n^{0})^{n-1}\\ &1,&(w_n^1)^1,&(w_n^1)^2,&\ldots,&(w_n^1)^{n-1}\\ &&&\vdots\\ &1,&(w_n^{n-1})^1,&(w_n^{n-1})^2,&\ldots,&(w_n^{n-1})^{n-1}) \end{bmatrix} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_{n-1} \end{bmatrix}= \begin{bmatrix} y_0\\y_1\\y_2\\\vdots\\y_{n-1} \end{bmatrix} ⎣⎢⎢⎢⎡1,1,1,(wn0)1,(wn1)1,(wnn−1)1,(wn0)2,(wn1)2,⋮(wnn−1)2,…,…,…,(wn0)n−1(wn1)n−1(wnn−1)n−1)⎦⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡a0a1a2⋮an−1⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡y0y1y2⋮yn−1⎦⎥⎥⎥⎥⎥⎤
我们的点值表示法所得到的矩阵是一个范德蒙矩阵,而我们进行IFFT的过程就是已知向量 yyy 以及点矩阵WWW 求向量aaa

那么我们只需要计算 yW−1yW^{-1}yW−1 即可,对于这个范德蒙矩阵求逆(之所以要指出它是范德蒙矩阵就是因为范德蒙矩阵一定可逆)。

直接构造矩阵S
[1,(wn−0)1,(wn−0)2,…,(wn−0)n−11,(wn−1)1,(wn−1)2,…,(wn−1)n−1⋮1,(wn−(n−1))1,(wn−(n−1))2,…,(wn−(n−1))n−1)]\begin{bmatrix} &1,&(w_n^{-0})^1,&(w_n^{-0})^2,&\ldots,&(w_n^{-0})^{n-1}\\ &1,&(w_n^{-1})^1,&(w_n^{-1})^2,&\ldots,&(w_n^{-1})^{n-1}\\ &&&\vdots\\ &1,&(w_n^{-(n-1)})^1,&(w_n^{-(n-1)})^2,&\ldots,&(w_n^{-(n-1)})^{n-1}) \end{bmatrix} ⎣⎢⎢⎢⎡1,1,1,(wn−0)1,(wn−1)1,(wn−(n−1))1,(wn−0)2,(wn−1)2,⋮(wn−(n−1))2,…,…,…,(wn−0)n−1(wn−1)n−1(wn−(n−1))n−1)⎦⎥⎥⎥⎤

根据求和引理,ci,j=wi,ksk,jc_{i,j}=w_{i,k}s_{k,j}ci,j=wi,ksk,j两个矩阵相乘就会得到矩阵C
[n,0,0,…,00,n,0,…,0⋱0,0,0,…,n]\begin{bmatrix} &n,&0,&0,&\ldots,&0\\ &0,&n,&0,&\ldots,&0\\ &&&\ddots\\ &0,&0,&0,&\ldots,&n \end{bmatrix} ⎣⎢⎢⎡n,0,0,0,n,0,0,0,⋱0,…,…,…,00n⎦⎥⎥⎤
这个矩阵乘上一个1n\frac{1}nn1就是单位矩阵,所以以wn−kw_n^{-k}wn−k作为点的横坐标依次代入矩阵再乘一个1n\frac{1}nn1就是我们需要的W−1W^{-1}W−1,所以只需要我们再用wn−kw_n^{-k}wn−k做一遍FFT,将结果乘1n\frac{1}nn1即可。

这样一来,我们用点值表示法就能够O(nlogn)O(nlogn)O(nlogn)快速求出系数表示了。

蝴蝶变换

蝴蝶变换是FFT的一个优化,之所以要优化,是因为纯递归求解FFT效率并不高,所以考虑如何迭代求解。

每次奇偶分组的时候,我们都会挑一个当前所在位置编号iii 为偶数的分到左边,为奇数的分到右边。

偶数的偶数的偶数次项到最后,就会是0,多一个奇数就会在头多1…（感性理解）

总之从变换的规律来看,我们递归到最后,系数的分组是原来系数的二进制反转

比如我们有一个7次多项式,有8项系数a0−a7a_0-a_7a0−a7

分组得到的就是(根据指数的奇偶性分组)
0,1,2,3,4,5,6,70,2,4,6,1,3,5,70,4,2,6,1,5,3,70,4,2,6,1,5,3,7二进制表示下就是000,001,010,011,100,101,110,111000,010,100,110,001,011,101,111000,100,010,110,001,101,011,111000,100,010,110,001,101,011,111\begin{matrix} 0,1,2,3,4,5,6,7\\ 0,2,4,6,1,3,5,7\\ 0,4,2,6,1,5,3,7\\ 0,4,2,6,1,5,3,7\\ 二进制表示下就是\\ 000,001,010,011,100,101,110,111\\ 000,010,100,110,001,011,101,111\\ 000,100,010,110,001,101,011,111\\ 000,100,010,110,001,101,011,111 \end{matrix} 0,1,2,3,4,5,6,70,2,4,6,1,3,5,70,4,2,6,1,5,3,70,4,2,6,1,5,3,7二进制表示下就是000,001,010,011,100,101,110,111000,010,100,110,001,011,101,111000,100,010,110,001,101,011,111000,100,010,110,001,101,011,111
你看最后一组的分配,就是每个数二进制反转的结果。

换言之,我们一开始如果将系数二进制反转,然后从下往上推,每次都能推出一个长度更长的多项式,这样子问题就可以一步步向上推出全局的答案。

终于,我们得到了整个FFT的算法(typroa显示共4862词(不包括下边的代码),码了一天,淦)！下面会在例题里贴上一个模板,然后随着做题会更新一些题目总结吧。

####资料参考：
【学习笔记】超简单的快速傅里叶变换（FFT）（含全套证明）
单位根-百度百科
OIwiki(里边还有较详细的二进制反转证明,我懒了不想写了)
百度上的傅里叶变换(主要看严谨概念来的,然鹅看不大懂,然后题目就多了个【OI向】)
[学校的非公开资料]这个就莫得连接

例题

【模板】多项式乘法（FFT）

代码如下：

点击此处展开和收起代码

   #include<bits/stdc++.h>#define reg registertypedef long long ll;using namespace std;inline int qr(){int x=0,f=0;char ch=0;while(!isdigit(ch)){f|=ch=='-';ch=getchar();}while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}return f?-x:x;}const int maxn=3e6+199;const double PI=3.14159265358979;int n,m;struct Complex{double x,y;Complex operator+(const Complex &t)const {return (Complex){x+t.x,y+t.y};}Complex operator-(const Complex &t)const{return (Complex){x-t.x,y-t.y};} Complex operator* (const Complex &t)const {return (Complex){x*t.x-y*t.y,x*t.y+y*t.x};}}a[maxn],b[maxn];int rev[maxn],bit,tot;//反转的二进制数 void fft(Complex c[],int inv){//inv表示的是单位根转动的方向 for(reg int i=0;i<tot;i++){//蝴蝶变换！ if(i<rev[i]){//小于rev[i]时交换一次就好 swap(c[i],c[rev[i]]);}}for(reg int mid=1;mid<tot;mid<<=1){//从底层开始枚举,即从最小的子问题开始 Complex w1=(Complex){cos(PI/mid),inv*sin(PI/mid)};//求得单位根 for(reg int i=0;i<tot;i+=mid<<1){//枚举每一个子问题 Complex wk=(Complex){1,0};for(reg int j=0;j<mid;j++,wk=wk*w1){//扫左半边,此时大项目A的第i+j项就是原来处理出的左半边和右半边的结合,即之前的第i+j项和i+j+mid的组合,模拟递归 Complex x=c[i+j],y=wk*c[i+j+mid];//A1和A2 c[i+j]=x+y,c[i+j+mid]=x-y;//左半边和右半边的处理 }//c存储点值A(wk)=A1*(wk)+wk*A2(wk) }}}int main(){//  freopen(".in","r",stdin);// freopen(".out","w",stdout);n=qr();m=qr();for(reg int i=0;i<=n;i++) a[i].x=qr();for(reg int j=0;j<=m;j++) b[j].x=qr();while((1<<bit)<n+m+1) bit++;//找到第一个大于n+m+1的2的整次幂 tot=1<<bit;//扩展多项式长度 for(reg int i=0;i<tot;i++){rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//记录每个数的二进制反转 } //rev记录的是每个二进制数的二进制逆转,我们这个递推是求出来除去第0位的逆再把第0位放到最高位,就能够得到这个数的逆 fft(a,1),fft(b,1);//转化为点 for(reg int i=0;i<tot;i++) a[i]=a[i]*b[i];//求出目标多项式的点值表示 fft(a,-1);//求逆转化 for(reg int i=0;i<=n+m;i++){printf("%d ",(int)(a[i].x/tot+0.5));//除去tot,就像我们之前分析的 ,+0.5来以防精度出问题 } return 0;}