FFT算法详解

zoukankan html css js c++ java

FFT算法详解
啊…本来觉得这是个比较良心的算法~~没想到这么抽搐~~这个算法真是将一个人的自学能力锻炼到了极致

好的，那我们就开始我们的~~飞飞兔~~

偷偷说一句，

写在前面

·好多你不理解的地方在代码里就只有半行。

·三个引理中，只有消去引理跟算法的实现没有关系——消去引理主要是用来证明的

·不是特别熟悉矩阵的同学，看到“矩阵”字样时不要慌，矩阵只是辅助工具，跟算法的实现关系不大，但是却是理论的基石。

多项式乘法（卷积

对于多项式 $sumlimits_{i = 0}^{n}{a_ix^i} = a_0 + a_1x + a_2x^2.....+a_nx^nA(x)=i=0∑naixi=a0+a1x+a2x2.....+anxnB(x)=∑i=0nbixi=b0+b1x+b2x2.....+bnxnB(x) = sumlimits_{i = 0}^{n}{b_ix^i} = b_0 + b_1x + b_2x^2.....+b_nx^nB(x)=i=0∑nbixi=b0+b1x+b2x2.....+bnxn$

那么显然的话，两者相乘之后得到 $CCC，他应该有ci=∑j=0iajbi−jc_i = sumlimits_{j=0}^{i}{a_jb_{i-j}} ci=j=0∑iajbi−j此处隐含着多项式是按次数的单调性给出的qwq.qwq.qwq.$

不难看出

$emmmmemmmmemmmm那么很显然的是，朴素的多项式乘法的时间复杂度是O(n2)O(n^2)O(n2)，不是特别可以接受，于是我们思考优化，便引入了DFTDFTDFT和IDFTIDFTIDFT这两个毒瘤可爱的操作qwqqwqqwq$

$First_stepmathcal{First\_step}First_step·换个思路思考多项式$

好的，我们现在的问题是，对于这个 $O(n2)O(n^2)O(n2)的算法进行优化。我们考虑系数乘法的复杂度，一定会是O(n2)O(n^2)O(n2)的，因为其中的每一步都不可忽略，所以我们期望换个思路来搞一搞。$

我们知道，原本的多项式是系数表示法，现在我们换个思路，将其转化为点值表示法。即我们可以把多项式 $f(x)f(x)f(x)看作在平面直角坐标系上的函数f(x)f(x)f(x)，那么这个nnn阶函数就可以由n+1n+1n+1个点唯一确定$（$不了解的的右转选修4−64-64−6)，即f(x)<=>(x0,y0),(x1,y1),(x2,y2)....(xn,yn)f(x)<=>{(x_0,y_0),(x_1,y_1),(x_2,y_2)....(x_n,y_n)}f(x)<=>(x0,y0),(x1,y1),(x2,y2)....(xn,yn)那么∀k,yk=f(xk)forall k,y_k = f(x_k)∀k,yk=f(xk)这是很显然的，并且学过函数的同学们应该知道，这n+1n+1n+1个点是随意选取的——只要求它们相异即可。$

前置技能一·点乘法：

假设我们有两个关于 ${(x_0,y_0),(x_1,y_1)....(x_n,y_n)}A(x)=(x0,y0),(x1,y1)....(xn,yn)B(x)=(x0,y0′),(x1,y1′)....(xn,yn)B(x) = {(x_0,y_0'),(x_1,y_1')....(x_n,y_n)}B(x)=(x0,y0′),(x1,y1′)....(xn,yn)那么他们相加的话，O(N)O(N)O(N)即可qwqqwqqwq A(x)+B(x)=(x0,y0+y0′),(x1,y1+y1′)(xn,yn+yn′)A(x)+B(x) = {(x_0,y_0+y_0'),(x_1,y_1+y_1')(x_n,y_n+y_n')}A(x)+B(x)=(x0,y0+y0′),(x1,y1+y1′)(xn,yn+yn′)对于多项式乘法，我们不能传统地按位相加，因为上文中我们提到过C(x)C(x)C(x)的项数为2n+12n+12n+1，所以我们考虑补上一堆项（并不是空项），让他们可以像多项式加法一样直接按位运算，类似这样A(x)=(x0,y0),(x1,y1)....(x2n,y2n)A(x) = {(x_0,y_0),(x_1,y_1)....(x_{2n},y_{2n})}A(x)=(x0,y0),(x1,y1)....(x2n,y2n)B(x)=(x0,y0′),(x1,y1′)....(x2n,y2n)B(x) = {(x_0,y_0'),(x_1,y_1')....(x_{2n},y_{2n})}B(x)=(x0,y0′),(x1,y1′)....(x2n,y2n) 做乘法可得 A(x)B(x)=(x0,y0y0′),(x1,y1y1′)(x2n,y2ny2n′)A(x)B(x) = {(x_0,y_0y_0'),(x_1,y_1y_1')(x_{2n},y_{2n}y_{2n}')}A(x)B(x)=(x0,y0y0′),(x1,y1y1′)(x2n,y2ny2n′) 所以我们考虑补项，之后再做点乘，得到C(x)C(x)C(x)的点乘表达式。$

我们观察点乘法，它的时间复杂度达到了惊人的

得到 ——>对于每个多项式，选取 $n+1n+1n+1个点，得出点值表达式(O(n2))(O(n^2))(O(n2)) ——>点乘法(O(n))(O(n))(O(n))——>将得出来的C(x)C(x)C(x)的点值表达式再转换成系数表达式(O(n2))(O(n^2))(O(n2))$

这就是

虽然其余部分的时间复杂度还是很麻烦的 $O(n2)O(n^2)O(n2)，但是都是可以优化成O(nlogn)O(nlogn)O(nlogn)的，我们先来看比较朴素的剩余步骤。因为对于第一步，我们可以用秦九韶算法，O(n)O(n)O(n)求出一个点，O(n2)O(n^2)O(n2)求出所有点。那么现在我们来思索第三步：$

前置技能二·插值法

看过选修

求值计算的逆（从一个多项式的点值表示确定其系数表示中的系数）称为插值(

首先我们要明确的一点是，多项式插值具有唯一性，~~你可以考虑使用“显然意会”法~~我们接下来证明一下：

我们证明这玩意儿类似于证明矩阵可逆，比如我们借助增广矩阵的相关知识可以得出有如下：

$x_0 & x_0^2 & cdots & x_0^n \ 1 & x_1 & x_1^2 & cdots & x_1^n \ 1 & x_2 & x_2^2 & cdots & x_2^n \ vdots & vdots & vdots & ddots & vdots \ 1 & x_n & x_n^2 & cdots & x_n^n \ end{vmatrix}∣∣∣∣∣∣∣∣∣∣111⋮1x0x1x2⋮xnx02x12x22⋮xn2⋯⋯⋯⋱⋯x0nx1nx2n⋮xnn∣∣∣∣∣∣∣∣∣∣∣a0a1a2⋮an∣egin{vmatrix} a_0 \ a_1 \ a_2 \ vdots \ a_n end{vmatrix}∣∣∣∣∣∣∣∣∣∣a0a1a2⋮an∣∣∣∣∣∣∣∣∣∣ = ∣y0y1y2⋮yn∣egin{vmatrix} y_0 \ y_1 \ y_2 \ vdots \ y_n end{vmatrix}∣∣∣∣∣∣∣∣∣∣y0y1y2⋮yn∣∣∣∣∣∣∣∣∣∣$

这可以很显然地表达出我们的第一步——将系数表达式转化为点值表达式。

我们观察最左边的矩阵，没错，他就是范德蒙德矩阵，记作 $V(x0,x1⋯xn)V(x_0,x_1 cdots x_n)V(x0,x1⋯xn)。而对于范德蒙德矩阵，他的行列式计算公式为det(V)=∏0≤j<k≤n(xk−xj)det(V) = prodlimits_{0 leq j<k leq n}^{}{(x_k - x_j)}det(V)=0≤j<k≤n∏(xk−xj)而因为xxx们两两互异，所以det(v)̸=0det(v) ot= 0det(v)̸=0 那我们考虑在弦形袋鼠中有讲过，这个矩阵的行列式非000即可逆，所以我们可以得到唯一的系数表达A⃗vec{A}A，使得A⃗=V(x0,x1⋯xn)−1×y⃗vec{A} = V(x_0,x_1 cdots x_n)^{-1} imes vec{y}A=V(x0,x1⋯xn)−1×y$

证毕。

然而这样的话，我们在插值的时候可以考虑直接求出范德蒙德矩阵的逆矩阵。 $emmmmmemmmmmemmmmm但是我们惊奇的发现对矩阵求逆是O(n3)O(n^3)O(n3)的！$

好吧，成功负优化了。

嗯，兄弟莫慌，我们有拉格朗日插值公式呢！所以我们可以 $O(n2)O(n^2)O(n2)求出来这个多项式$

$sumlimits_{k=0}^{n}{y_kfrac{prodlimits_{j ot= k}^{}{(x-x_j)}}{prodlimits_{j ot= k}^{}{(x_k-x_j)}}}F(x)=k=0∑nykj̸=k∏(xk−xj)j̸=k∏(x−xj)$

哈，我们成功地追平了没有优化时的、渐进意义下的时间复杂度！

~~那这优化好像除了增添一个贼大的常数啥都不能干~~

嗯，读到这儿，你已经会了

那么对这个东西进行加速的就是啦！

$Second_stepmathcal{Second\_step}Second_step·FFTFFTFFT究竟是如何优化的$

下面来解释真相——为何这种变换策略可以称之为优化？

因为实际上，我们的第一步——求值（系数转点值）和我们的第三步（点值转系数）都是可以做到

下面就让我们来看看如何优化：

前置技能三·单位复根

$nnn次单位复根是满足wn=1w^n = 1wn=1 的复数www,其中我们可以由复数的运算法则（辐角相乘，模相加），我们可以很简单地得出nnn次单位根有nnn个这个结论，而又因为复数wnw^nwn在复数平面上的模都是一，所以相乘之后还会是一，那么所有的wi,1≤i≤<nw_i,1 leq i leq <nwi,1≤i≤<n就会均匀分布在单位圆上，类似当n=8n = 8n=8时它是这样的：$

图源自于《算法导论》。

我们考虑欧拉公式：

不了解的戳这里，文章的第四个标题 $eix=cosx+isinxe^{ix} = cosx + isinxeix=cosx+isinx$

我们取 $=2pix=2π，可以得到如下关系式：e2πi=1=wne^{2 pi i} = 1 = w^ne2πi=1=wn从而w=e2πinw = e^{frac{2pi i}{n}}w=en2πi$

我们把此时的单位根称之为主次单位根，记作 $wn=e2πinw_n = e^{frac{2pi i}{n}} wn=en2πi那么对于其他的单位根，记作wnk=e2πikn,0≤k<nw_n^k=e^{frac{2pi ik}{n}},0 leq k < nwnk=en2πik,0≤k<n都是主次单位根的整次幂。大家可以考虑借助上面那张图图意会一下qwq.qwq.qwq.$

诶，这个有啥用啊

那是因为单位根们有一堆特别好用的性质，让我们可以将数据规模不断折半，使得其达到

那么我们先来看其支持其规模减半的引理：

消去引理

引理：对任何整数 $>0n≥0,k≥0,d>0,有wdndk=wnkw_{dn}^{dk} = w_n^kwdndk=wnk这个好像很好证的样子……代入定义可以获得wdndk=e2πdkdn=e2πikn=wnk,w_{dn}^{dk} = e^{frac{2pi dk}{dn}} = e^{frac{2pi ik}{n}} = w_n^k,wdndk=edn2πdk=en2πik=wnk,那么由此我们可以推导出下面的引理：$

折半引理

引理：对于任何大于

我们可以由消去引理得到 $(wnk)2=wn/2k(w_n^k)^2 = w_{n/2}^k(wnk)2=wn/2k这个……也可以尝试意会。那么接下来，如果对所有的nnn次单位跟平方一下，我们会发现n2frac{n}{2}2n次单位根每个都恰好出现了两次，证明如下：(wnk+n2)2=wn2k+n=wn2k×wnn=wn2k=(wnk)2(w_n^{k + frac{n}{2}})^2 = w_n^{2k + n} = w_n^{2k} imes w_n^n = w_n^{2k} = (w_n^k)^2(wnk+2n)2=wn2k+n=wn2k×wnn=wn2k=(wnk)2我们参考上面那张图，可以发现它的意义就在于方向相对的两组向量，其平方相等。基于这两个式子便可以保证我们时间O(nlogn)O(nlogn)O(nlogn)。$

嗯，好的，猜都能猜到你现在一定认为这玩意儿……让人云里雾里的。那简单总结一下吧，我们现在根据这个特别易证~~但根本无法理解~~的引理得到了如下两个公式： $(wnk+n2)2=(wnk)2(w_n^{k + frac{n}{2}})^2 = (w_n^k)^2(wnk+2n)2=(wnk)2(wnk)2=wn/2k(w_n^k)^2 = w_{n/2}^k(wnk)2=wn/2k这两个公式被我们所证，一定有它独特的理由，那么证明一的原因是，据此我们可以发现有wn0=wnn2w_n^0 = w_n^{frac{n}{2}}wn0=wn2n，wn1=wnn2+1w_n^1 = w_n^{frac{n}{2} + 1}wn1=wn2n+1⋯cdots⋯，把所有的单位根画到一个数列上就是这样$

但是事实上，只有第一个式子我们并不能说明太多，只能说明 $LrangeLrangeLrange和RrangeRrangeRrange具有某种映射关系，但我们不知道到底是否是单射关系。即迄今为止，我们在上图中还不可以确定那些红色的带箭头的连线是否恰好正确。这个时候我们就需要公式二(wnk)2=wn/2k(w_n^k)^2 = w_{n/2}^k(wnk)2=wn/2k这才可以保证，序列左边与右边具有一一对应的映射关系，并且映射序列满足单调性，即n2+1frac{n}{2}+12n+1号单位根只能和111号对应。$

而我们在代码实现中，不能直接得到，所以这个时候求单位根的任务就交给了我们上文中提到过的欧拉公式。

求和引理

引理：对于任意 $n>0n>0n>0且kkk不能整除nnn，我们都有∑j=0n−1(wnk)j=0sumlimits_{j =0}^{n-1}{(w_n^k)^j} = 0j=0∑n−1(wnk)j=0$

由几何级数的求和公式（等比数列求和公式) $∑j=0nxj=xj+1−1x−1sumlimits_{j = 0}^{n}{x^j} = frac{x^{j +1} -1}{x -1}j=0∑nxj=x−1xj+1−1可得∑j=0n−1(wnk)j=(wnk)n−1wnk−1=(wnn)k−1wnk−1=(1)k−1wnk−1sumlimits_{j =0}^{n-1}{(w_n^k)^j} = frac{(w_n^k)^n -1}{w_n^k -1} = frac{(w_n^n)^k -1}{w_n^k -1} = frac{(1)^k -1}{w_n^k -1}j=0∑n−1(wnk)j=wnk−1(wnk)n−1=wnk−1(wnn)k−1=wnk−1(1)k−1由于保证了kkk不可整除nnn所以分母一定不为0.0.0.$

真正地认识

那么我们在了解完单位复数根之后，便可以正式地对

对于我们已知的一个多项式 $sumlimits_{i =0}^{n - 1}{a_ix^i}A(x)=i=0∑n−1aixi在wn0,wn1,wn2⋯wnn−1w_n^0,w_n^1,w_n^2 cdots w_n^{n-1}wn0,wn1,wn2⋯wnn−1处的取值，我们可以假定nnn是222的幂，因为即使它本身不是222的幂，我们也可以通过向高次幂补值为000的项来解决这个问题。而补足222的幂的目的，就是为了在FFTFFTFFT分治的过程中，使之可以一直分治下去且每次分治得出的两半可以进行均衡运算。$

那我们现在会有一个 ${a_1, a_2, a_3 cdots a_{n-1}}a=a1,a2,a3 cdotsan−1，对于k=0,1,2,⋯n−1k = 0, 1, 2, cdots n -1k=0,1,2,⋯n−1，定义yky_kyk如下：yk=A(wnk)=∑j=0n−1aj×wnkjy_k = A(w_n^k)=sumlimits_{j =0}^{n -1}{a_j imes w_n^{kj}}yk=A(wnk)=j=0∑n−1aj×wnkj，那么向量y⃗=y0,y1,y2⋯yn−1vec{y} = {y_0, y_1, y_2 cdots y_{n-1}}y=y0,y1,y2⋯yn−1就称作系数向量a⃗=a1,a2,a3⋯an−1vec{a} = {a_1, a_2, a_3 cdots a_{n-1}}a=a1,a2,a3⋯an−1的离散型傅立叶变换（Discrete Fourier TransformationDiscrete Fourier TransformationDiscrete Fourier Transformation）。$

嗯，这个离散型我们可以由点乘法联想意会一下：本来

在上文中我们分析过，将系数表达式转化为点值表达式需要的时间复杂度为 $O(n2)O(n^2)O(n2)，这是朴素算法。而我们只需要用一种被称作快速傅立叶变换（Fast Fourier TransformationFast Fourier TransformationFast Fourier Transformation）的方式，就可以将其时间复杂度压缩成O(nlogn)O(nlogn)O(nlogn)。而在这里我们就用到了刚才证明的引理——折半引理。$

我们考虑将原来的多项式 $a_0+a_1x+ a_2x^2 cdots +a_{n-1}x^{n-1}A(x)=a0+a1x+a2x2⋯+an−1xn−1重定义成两个次数域为n2frac{n}{2}2n的小多项式A[0](x)A^{[0]}(x)A[0](x)和A[1](x)A^{[1]}(x)A[1](x)：$

$A[0](x)=a0+a2x+a4x2⋯+an−2xn2−1A^{[0]}(x) = a_0 + a_2x+a_4x^2 cdots +a_{n-2}x^{frac{n}{2} - 1}A[0](x)=a0+a2x+a4x2⋯+an−2x2n−1A[1](x)=a1+a3x+a5x2⋯+an−1xn2−1 A^{[1]}(x) = a_1 + a_3x+a_5x^2 cdots +a_{n-1}x^{frac{n}{2} - 1}A[1](x)=a1+a3x+a5x2⋯+an−1x2n−1 那么也就是说，A[0](x)A^{[0]}(x)A[0](x)存储的是所有偶数位（二进制位最后一位是000），而A[1](x)A^{[1]}(x)A[1](x)存储的是所有的奇数位（二进制位最后一位是111），那么有下式：A(x)=A[0](x2)+xA[1](x2)A(x) = A^{[0]}(x^2)+xA^{[1]}(x^2)A(x)=A[0](x2)+xA[1](x2)那我们求A(x)A(x)A(x)在单位根们wn0,wn1,wn2⋯,wnn−1w_n^0,w_n^1,w_n^2 cdots ,w_n^{n-1}wn0,wn1,wn2⋯,wnn−1处的值，就变成了先求出A[0](x2)A^{[0]}(x^2)A[0](x2)和A[1](x2)A^{[1]}(x^2)A[1](x2)的值，然后根据上式进行合并即可。$

而显然的是，根据折半引理，我们根本不需要

那么同时对于另一边，我们也有公式可以计算出来
int Lim = 1, N, M ; function FFT(int lenth, complex *A, int flag){ IF (Lim == 1) return ; complex A0[lenth >> 1], A1[lenth >> 1] ;//分成两部分 for(int j : 0 to lenth by_grow 2) A0[j >> 1] = A[j], A1[j >> 1] = A[j + 1] ; FFT(lenth >> 1, A0, flag) ; FFT(lenth >> 1, A1, flag) ; complex Wn = unit(,) , w = (1, 0) ;//Wn是单位根，w用来枚举幂，即我们令主次单位根不变，由于其余单位根都是其整次幂，所以可以用一个w来记录到第几次幂 /*此处求单位根的时候会用到我们的参数flag……嗯没错就用这一次，并且flag的值域为(-1, 1)……是的，只会有两个值*/ for(int j : 0 to (lenth >> 1) by_grow 1 with w = w * Wn){ A[i] = A0[i] + A1[i] * w ;//应用公式，下同 A[i + (lenth >> 1)] = A0[i] - A1[i] * w ; //顺便求出另一半，由折半引理可显然。 } } function Main{ input(N), input(M) ; for(i : 0 to N by_grow 1) => input(A) ; for(i : 0 to M by_grow 1) => input(B) ; while(Lim < N + M) Lim <<= 1 ;//Lim为结果多项式的长度（暂时），化为2的幂的便于分治（二分） FFT(Lim, A, 1) ;//两遍FFT表示从系数化为点值 FFT(Lim, B, 1) ; for(i : 0 to Lim by_grow 2) => A[i] *= B[i] ;//点乘法，此处需要重定义乘号，因为每一项现在表示的是一个点，有x和y两个属性qwq }
以上是基于

嗯，这就是求值的方法，好像很

上文中我们曾经提及过的范德蒙德矩阵可以放到这儿用：

$\ 1 & w_n & w_n^2 & cdots & w_n^{n-1} \ 1 & w_n^2 & w_n^4 & cdots & w_n^{2(n-1)} \ vdots & vdots & vdots & ddots & vdots \ 1 & w_n^{n-1} & w_n^{2(n-1)} & cdots & w_n^{(n-1)(n-1)} \ end{vmatrix}∣∣∣∣∣∣∣∣∣∣111⋮11wnwn2⋮wnn−11wn2wn4⋮wn2(n−1)⋯⋯⋯⋱⋯1wnn−1wn2(n−1)⋮wn(n−1)(n−1)∣∣∣∣∣∣∣∣∣∣∣a0a1a2⋮an−1∣egin{vmatrix} a_0 \ a_1 \ a_2 \ vdots \ a_{n-1} end{vmatrix}∣∣∣∣∣∣∣∣∣∣a0a1a2⋮an−1∣∣∣∣∣∣∣∣∣∣ = ∣y0y1y2⋮yn−1∣egin{vmatrix} y_0 \ y_1 \ y_2 \ vdots \ y_{n-1} end{vmatrix}∣∣∣∣∣∣∣∣∣∣y0y1y2⋮yn−1∣∣∣∣∣∣∣∣∣∣$

那为了求出我们的 ${a_0, a_1 cdots ,a_{n-1}}a=a0,a1⋯,an−1我们应该让刚刚求值算出的y⃗vec{y}y乘上我们V⃗−1vec{V}^{-1}V−1（V⃗vec{V}V的逆矩阵）即可。但是桥豆麻袋~~~不需要用什么高消啊…余子式啊…我们只需要以下：$

定理：对于 $n-1,V_n^{-1}j,k=0,1,2⋯n−1,Vn−1的(j,k)(j, k)(j,k)处的值为wn−kj/nw_n^{-kj}/nwn−kj/n 证明：我们考虑反向证明，已知Vn′V_n'Vn′是一个(j,k)(j,k)(j,k)处值为wn−kj/nw_n^{-kj}/nwn−kj/n，那我们只需要证明V′×V=InV' imes V = I_nV′×V=In即可，其中InI_nIn是单位矩阵，即主对角线都是111，其余位置上是000的矩阵。qwqqwqqwq这一步转化应该挺简单的吧，至于不知道单位矩阵是啥的…百度去QAQQAQQAQ。那么我们考察V′VV' VV′V中的元素(i,j)(i, j)(i,j)有如下的式子V′×V=∑k=0n−1(wn−ki/n)×wnkj=∑k=0n−1wnk(j−i)/nV' imes V = sumlimits^ {n-1}_{k=0}{(w_n^{-ki}/n)} imes {w_n^{kj}} = sumlimits^ {n-1}_{k=0}{w_n^{k(j-i)}/n}V′×V=k=0∑n−1(wn−ki/n)×wnkj=k=0∑n−1wnk(j−i)/n由求和引理当且仅当i=ji=ji=j时其值为一，其余的时刻均为零，正好和我们单位矩阵的样子一毛一样！我们可以想到单位矩阵InI_nIn里面也是当i=ji =ji=j时，在主对角线上其值为000. 于是乎，证毕。$

那么我们把我们刚刚求出来的逆矩阵 $V−1V^{-1}V−1美化一下，提出每一项所除的nnn，可以得到IDFTIDFTIDFT可以如此求：IDFTn(y)=1n∑k=0n−1ykwn−kjIDFT_n(y) = frac{1}{n}sumlimits_{k = 0}^{n-1}{y_kw_n^{-kj}}IDFTn(y)=n1k=0∑n−1ykwn−kj诶，这个好像……跟我们求值时的公式差不多？没错，除了带个负号，其余的都差不多,所以这就可以解释我们刚才的伪代码中读进去的参数flagflagflag看似没有什么用，现在知道了？当flag=1flag=1flag=1时，他是正向DFTDFTDFT；当它等于−1-1−1时，它是逆向DFTDFTDFT。这可以让我们通过这一个函数解决两个过程。我们只需要用yyy替换aaa，用wn−1w_n^{-1}wn−1替换wnw_nwn，其余的没什么差别，于是……时间复杂度还是O(nlogn)O(nlogn)O(nlogn)的！$

那么就亮出我们的
void FFT(int Lim,complex *A,int flag){ if(Lim == 1) return ; complex A0[Lim >> 1], A1[Lim >> 1] ; for(int i = 0; i <= Lim ; i += 2) A0[i >> 1] = A[i], A1[i >> 1] = A[i+1] ; FFT(Lim >> 1, A0, flag) ; FFT(Lim >> 1, A1, flag) ; complex unit = complex(cos(2.0 * Pi / Lim) , flag * sin(2.0 * Pi / Lim)), w = complex(1, 0) ;//欧拉公式 for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) { A[i] = A0[i] + w * A1[i] ; A[i + (Lim>>1)] = A0[i] - w*A1[i]; } } int main(){ ...................... FFT(A, 1), FFT(B, 1) ; for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ; FFT(A, -1) ; ...................... }
好的，现在嘛……可以考虑撒花花啦！因为我们的

$Third_stepmathcal{Third\_step}Third_step·最后来一些优化吧！$

我们现在要引出的就是迭代版的

前置技能四·小
for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) { a[i] = A0[i] + w * A1[i] ; a[i + (Lim>>1)] = A0[i] - w*A1[i]; }
我们会发现……
for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) { int temp = w * A1[i] ; a[i] = A0[i] + t ; a[i + (Lim>>1)] = A0[i] - t ; }
嗯，这就是全部的优化啦！那么，

我们定义 $wnkw_n^kwnk为旋转因子，那么每一次我们先将yk[1]y_k^{[1]}yk[1]与旋转因子的乘积存储在一个变量ttt里，并在yk[0]y_k^{[0]}yk[0]增加、减去ttt的操作称为一次蝴蝶操作。$

在了解完蝴蝶操作之后，我们首先考虑按照递归的思路，将

我们会发现，其实我们是可以对它进行反向迭代的。以上面的迭代树为例，我们的步骤大体如下：

而反向迭代似乎有规律可循。我们发现只要我们用迭代的过程模拟出回溯的过程即可。那么思路便有了：三层
for(j = 1; j < Lim; j <<= 1){ node T(cos(Pi / j), flag * sin(Pi / j)) ; for(k = 0; k < Lim; k += (j << 1) ){ node t(1, 0) ; for(l = 0 ; l < j; l ++, t = t * T){ node Nx = J[k + l], Ny = t * J[k + j + l] ; J[k + l] = Nx + Ny ; J[k + j + l] = Nx - Ny ; } } }
嗯，好像……海星？哈，思维不严谨的我们漏了一个地方：我们并不可以得到迄今为止的结构体参数

前置（最后的）技能·蝴蝶定理

这个嘛……我们可以选择打个表观察：原来的序号

诶，二进制好像是反序的嗷~~这便是我们的最后一个

因为既然是重新排序，所以必然有 $ji<j但AiA_iAi反序之后到了AjA_jAj的后面。因为我们观察到的蝴蝶定理是满足一一对应性的，所以我们在FFTFFTFFT之前swapswapswap一遍即可。$

嗯，然后我们可以将这个反序存在一个数组里面。类似这样求出来：
```
    for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;
```
那么我们可以看到，这就简化了很多冗余的步骤，并让我们脱离递归的大常数。~~真开森啊~~

最后附迭代版的代码(我写的常数好像有点儿大
#include <cmath> #include <cstdio> #include <iostream> #define il inline using namespace std ; int N, M, K ; const int MAXN = 3000100 ; const double Pi = acos(-1.0) ; int i, j, k, l, Lim = 1, L, R[MAXN] ; struct node{ double x, y ; node (double xx = 0, double yy = 0){ x = xx, y = yy ; } }A[MAXN], B[MAXN] ; node operator * (node J, node Q){ return node(J.x * Q.x - J.y * Q.y , J.x * Q.y + J.y * Q.x); } node operator + (node J, node Q){ return node(J.x + Q.x , J.y + Q.y); } node operator - (node J, node Q){ return node(J.x - Q.x , J.y - Q.y ); } il int qr(){ int k = 0, f = 1 ; char c = getchar() ; while(!isdigit(c)){if(c == '-') f = -1 ;c = getchar() ;} while(isdigit(c)) k = (k << 1) + (k << 3) + c - 48 ,c = getchar() ; return k * f ; } void FFT(node *J, int flag){ for(i = 0; i < Lim; i ++) if(i < R[i]) swap(J[i], J[R[i]]) ; for(j = 1; j < Lim; j <<= 1){ node T(cos(Pi / j), flag * sin(Pi / j)) ; for(k = 0; k < Lim; k += (j << 1) ){ node t(1, 0) ; for(l = 0 ; l < j; l ++, t = t * T){ node Nx = J[k + l], Ny = t * J[k + j + l] ; J[k + l] = Nx + Ny ; J[k + j + l] = Nx - Ny ; } } } } int main(){ N = qr(), M = qr() ; for(i = 0; i <= N; i ++) A[i].x = qr() ; for(i = 0; i <= M; i ++) B[i].x = qr() ; while(Lim <= N + M) Lim <<= 1, L ++ ; for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ; FFT(A, 1), FFT(B, 1) ; for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ; FFT(A, -1) ; for(i = 0; i <= N + M; i ++) printf("%d ", (int)(A[i].x / Lim + 0.5)) ;//我们推过的公式里面有一个1/n这一项，最后输出的时候添上即可 return 0 ; }
啊……那就撒花花吧！！

$Last_stepmathcal{Last\_step}Last_step·感悟$

嗯……怎么说呢，现在看这个算法，真是简单的一匹啊……代码这么短~~这么容易背过~~。但是当时理解起来却花了很大心思呢！这篇博客我写了整整三天

鸣谢以下对我有帮助的典籍或者

·《算法导论》（讲的……还可以吧）

·

·

$Flower_pksby mathcal{color{cyan}{Flower}color{purple}{\_}color{pink}{pks}}by Flower_pks$
查看全文

相关阅读:
独立人格，让人生更美好
 版本控制入门简介
 Android蓝牙操作
 Android Rect和RectF的区别
 做个环保主义的程序员
 java/android下JNI编程总结
 别太单纯，也别太不单纯
 C++ namespace的用法
 JNI编程
 企业架构－架构原则

原文地址：https://www.cnblogs.com/zzrblogs/p/10638747.html

写在前面

多项式乘法（卷积

First_stepmathcal{First\_step}First_step·换个思路思考多项式

前置技能一·点乘法：

前置技能二·插值法

Second_stepmathcal{Second\_step}Second_step·FFTFFTFFT究竟是如何优化的

前置技能三·单位复根

消去引理

折半引理

求和引理

真正地认识DFTDFTDFT与FFTFFTFFT

DFTDFTDFT

FFTFFTFFT是如何求值的

FFTFFTFFT是如何插值的

Third_stepmathcal{Third\_step}Third_step·最后来一些优化吧！

前置技能四·小tricktricktrick与“蝴蝶操作”

前置（最后的）技能·蝴蝶定理

Last_stepmathcal{Last\_step}Last_step·感悟

鸣谢以下对我有帮助的典籍或者dalaodalaodalao

by Flower_pksby mathcal{color{cyan}{Flower}color{purple}{\_}color{pink}{pks}}by Flower_pks

$First_stepmathcal{First\_step}First_step·换个思路思考多项式$

$Second_stepmathcal{Second\_step}Second_step·FFTFFTFFT究竟是如何优化的$

真正地认识

$Third_stepmathcal{Third\_step}Third_step·最后来一些优化吧！$

前置技能四·小

$Last_stepmathcal{Last\_step}Last_step·感悟$

鸣谢以下对我有帮助的典籍或者

$Flower_pksby mathcal{color{cyan}{Flower}color{purple}{\_}color{pink}{pks}}by Flower_pks$