Live2D

FFT算法详解

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

好的，那我们就开始我们的飞飞兔FFTFFTFFT算法吧！

偷偷说一句，FFTFFTFFT的代码十分的短哦~并且如果你不喜欢看算法，你可以翻到最下面看心得哟！

写在前面

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

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

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

多项式乘法（卷积

对于多项式AAA和多项式BBB,两者是关于xxx的nnn次多项式，算上常数项均有N+1N+1N+1项，那么他们的简单版表示就是A(x)=∑i=0naixi=a0+a1x+a2x2.....+anxnA(x) = \sum\limits_{i = 0}^{n}{a_ix^i} = a_0 + a_1x + a_2x^2.....+a_nx^nA(x)=i=0∑n​ai​xi=a0​+a1​x+a2​x2.....+an​xnB(x)=∑i=0nbixi=b0+b1x+b2x2.....+bnxnB(x) = \sum\limits_{i = 0}^{n}{b_ix^i} = b_0 + b_1x + b_2x^2.....+b_nx^nB(x)=i=0∑n​bi​xi=b0​+b1​x+b2​x2.....+bn​xn

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

不难看出C(x)C(x)C(x)一共由2n+12n+12n+1项。如果你单纯地用FFTFFTFFT来做那些友好的题的话，那这其实也无妨理解为CCC就是AAA和BBB的卷积。

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

First_step\mathcal{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个点是随意选取的——只要求它们相异即可。

前置技能一·点乘法：

假设我们有两个关于xxx的NNN次多项式A(x)和B(x)，他们的点乘表达式均已知：A(x)=(x0,y0),(x1,y1)....(xn,yn)A(x) = {(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​,y0​y0′​),(x1​,y1​y1′​)(x2n​,y2n​y2n′​) 所以我们考虑补项，之后再做点乘，得到C(x)C(x)C(x)的点乘表达式。

我们观察点乘法，它的时间复杂度达到了惊人的O(n)O(n)O(n)，完全可以接受。这便是为什么化系数乘法为点乘法的原因。那么现在我们已经解决了它的主要部分，算法的大体思路也有了：

得到A(x)A(x)A(x)和B(x)B(x)B(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))

这就是FFTFFTFFT的大体流程。转化之后怎么没多快常数还大了

虽然其余部分的时间复杂度还是很麻烦的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)求出所有点。那么现在我们来思索第三步：

前置技能二·插值法

看过选修4−64-64−6的小伙伴们对此一定不陌生。那么对于插值法的简要定义如下：

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

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

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

∣1x0x02⋯x0n1x1x12⋯x1n1x2x22⋯x2n⋮⋮⋮⋱⋮1xnxn2⋯xnn∣\begin{vmatrix} 1 & 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⋮1​x0​x1​x2​⋮xn​​x02​x12​x22​⋮xn2​​⋯⋯⋯⋱⋯​x0n​x1n​x2n​⋮xnn​​∣∣∣∣∣∣∣∣∣∣​∣a0a1a2⋮an∣\begin{vmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{vmatrix}∣∣∣∣∣∣∣∣∣∣​a0​a1​a2​⋮an​​∣∣∣∣∣∣∣∣∣∣​ = ∣y0y1y2⋮yn∣\begin{vmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \end{vmatrix}∣∣∣∣∣∣∣∣∣∣​y0​y1​y2​⋮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) = \prod\limits_{0 \leq j<k \leq n}^{}{(x_k - x_j)}det(V)=0≤j<k≤n∏​(xk​−xj​)而因为xxx们两两互异，所以det(v)̸=0det(v) \not= 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} \times \vec{y}A=V(x0​,x1​⋯xn​)−1×y​

证毕。

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

好吧，成功负优化了。

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

F(x)=∑k=0nyk∏j̸=k(x−xj)∏j̸=k(xk−xj)F(x) = \sum\limits_{k=0}^{n}{y_k\frac{\prod\limits_{j \not= k}^{}{(x-x_j)}}{\prod\limits_{j \not= k}^{}{(x_k-x_j)}}}F(x)=k=0∑n​yk​j̸=k∏​(xk​−xj​)j̸=k∏​(x−xj​)​

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

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

嗯，读到这儿，你已经会了DFTDFTDFT（离散傅立叶变换）和IDFTIDFTIDFT（逆离散傅立叶变换）了！实际上，DFTDFTDFT对应着的就是我们把系数表达式映射到点值表达式的过程，IDFTIDFTIDFT对应着的就是我们把点值表达式映射到系数表达式的过程！

那么对这个东西进行加速的就是FFT      Fast  Fourier  TransformationFFT \ \ \ \ \ \ Fast \ \ Fourier \ \ TransformationFFT      Fast  Fourier  Transformation啦！

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

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

因为实际上，我们的第一步——求值（系数转点值）和我们的第三步（点值转系数）都是可以做到nlognnlognnlogn的，那么总的时间复杂度，渐进意义下就是O(nlogn)O(nlogn)O(nlogn)的！大胜利！

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

前置技能三·单位复根

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

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

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

诶，这个有啥用啊QAQQAQQAQ?

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

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

消去引理

引理：对任何整数n≥0,k≥0,d>0n \geq 0,k \geq 0,d >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{2\pi dk}{dn}} = e^{\frac{2\pi ik}{n}} = w_n^k,wdndk​=edn2πdk​=en2πik​=wnk​,那么由此我们可以推导出下面的引理：

折半引理

引理：对于任何大于000的偶数nnn，都有nnn个nnn次单位复根的平方的集合，等于n2\frac{n}{2}2n​个n2\frac{n}{2}2n​次单位复根的集合。

我们可以由消去引理得到(wnk)2=wn/2k(w_n^k)^2 = w_{n/2}^k(wnk​)2=wn/2k​这个……也可以尝试意会。那么接下来，如果对所有的nnn次单位跟平方一下，我们会发现n2\frac{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} \times 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+1\frac{n}{2}+12n​+1号单位根只能和111号对应。

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

求和引理

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

由几何级数的求和公式（等比数列求和公式)∑j=0nxj=xj+1−1x−1\sum\limits_{j = 0}^{n}{x^j} = \frac{x^{j +1} -1}{x -1}j=0∑n​xj=x−1xj+1−1​可得∑j=0n−1(wnk)j=(wnk)n−1wnk−1=(wnn)k−1wnk−1=(1)k−1wnk−1\sum\limits_{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.

真正地认识DFTDFTDFT与FFTFFTFFT

那么我们在了解完单位复数根之后，便可以正式地对DFTDFTDFT给出定义与操作方案了。

DFTDFTDFT

对于我们已知的一个多项式A(x)=∑i=0n−1aixiA(x) = \sum\limits_{i =0}^{n - 1}{a_ix^i}A(x)=i=0∑n−1​ai​xi在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分治的过程中，使之可以一直分治下去且每次分治得出的两半可以进行均衡运算。

那我们现在会有一个AAA的向量组a⃗=a1,a2,a3 cdotsan−1\vec{a} = {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)=\sum\limits_{j =0}^{n -1}{a_j \times w_n^{kj}}yk​=A(wnk​)=j=0∑n−1​aj​×wnkj​，那么向量y⃗=y0,y1,y2⋯yn−1\vec{y} = {y_0, y_1, y_2 \cdots y_{n-1}}y​=y0​,y1​,y2​⋯yn−1​就称作系数向量a⃗=a1,a2,a3⋯an−1\vec{a} = {a_1, a_2, a_3 \cdots a_{n-1}}a=a1​,a2​,a3​⋯an−1​的离散型傅立叶变换（Discrete  Fourier  TransformationDiscrete \ \ Fourier \ \ TransformationDiscrete  Fourier  Transformation）。

嗯，这个离散型我们可以由点乘法联想意会一下：本来A(x)A(x)A(x)是一个优美的多项式，转变成函数之后是一条优美的曲线（优美只是定语……不是重要内容qwqqwqqwq），结果你突然把它拆成了一堆离散的点，把它用点值法表示，故称之曰：“离散型” 。

FFTFFTFFT是如何求值的

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

我们考虑将原来的多项式A(x)=a0+a1x+a2x2⋯+an−1xn−1A(x) = a_0+a_1x+ a_2x^2 \cdots +a_{n-1}x^{n-1}A(x)=a0​+a1​x+a2​x2⋯+an−1​xn−1重定义成两个次数域为n2\frac{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​+a2​x+a4​x2⋯+an−2​x2n​−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​+a3​x+a5​x2⋯+an−1​x2n​−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)的值，然后根据上式进行合并即可。

而显然的是，根据折半引理，我们根本不需要O(n)O(n)O(n)求，而是通过数据规模不断减小使之成为O(logn)O(logn)O(logn)。于是，我们成功通过FFTFFTFFT优化了求值的复杂度qwqqwqqwq

那么同时对于另一边，我们也有公式可以计算出来

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 
}

以上是基于pkspkspks标准下的伪代码你可以试试在c++标准下运行，其中forforfor循环部分，grow表示当前循环变量的单次增量，之后带有withwithwith表示每次循环结束都会进行的运算（下同

嗯，这就是求值的方法，好像很nicenicenice地达到了O(nlogn)O(nlogn)O(nlogn)

FFTFFTFFT是如何插值的

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

∣111⋯11wnwn2⋯wnn−11wn2wn4⋯wn2(n−1)⋮⋮⋮⋱⋮1wnn−1wn2(n−1)⋯wn(n−1)(n−1)∣\begin{vmatrix} 1 & 1 & 1 & \cdots & 1 \\ 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⋮1​1wn​wn2​⋮wnn−1​​1wn2​wn4​⋮wn2(n−1)​​⋯⋯⋯⋱⋯​1wnn−1​wn2(n−1)​⋮wn(n−1)(n−1)​​∣∣∣∣∣∣∣∣∣∣​∣a0a1a2⋮an−1∣\begin{vmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{vmatrix}∣∣∣∣∣∣∣∣∣∣​a0​a1​a2​⋮an−1​​∣∣∣∣∣∣∣∣∣∣​ = ∣y0y1y2⋮yn−1∣\begin{vmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \end{vmatrix}∣∣∣∣∣∣∣∣∣∣​y0​y1​y2​⋮yn−1​​∣∣∣∣∣∣∣∣∣∣​

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

定理：对于j,k=0,1,2⋯n−1,Vn−1j,k = 0,1, 2 \cdots 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' \times 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' \times V = \sum\limits^ {n-1}_{k=0}{(w_n^{-ki}/n)} \times {w_n^{kj}} = \sum\limits^ {n-1}_{k=0}{w_n^{k(j-i)}/n}V′×V=k=0∑n−1​(wn−ki​/n)×wnkj​=k=0∑n−1​wnk(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}\sum\limits_{k = 0}^{n-1}{y_kw_n^{-kj}}IDFTn​(y)=n1​k=0∑n−1​yk​wn−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)的！

那么就亮出我们的cppcppcpp吧！

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) ;
......................
}

好的，现在嘛……可以考虑撒花花啦！因为我们的FFTFFTFFT实际上已经结束了！~~~~ButButBut，这个递归版本的FFTFFTFFT的常数，由于牵扯到sin/cossin/cossin/cos的运算、doubledoubledouble、递归时的入栈出栈（底层），所以特别的大emmmmmemmmmmemmmmm，那我们该怎么办呢——

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

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

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

emmmemmmemmm先上一个不是特别卡常数的优化。我们观察之前的代码中，有这么一步：

   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];
}

我们会发现……w×Ai[i]w \times Ai[i]w×Ai[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 ;
    }

嗯，这就是全部的优化啦！那么，FFTFFTFFT，完！

qwqqwqqwq这分明是骗小孩子的啦……如果单单这一步就可以卡常数的话，那这个世界该多么美好QAQ\mathcal{QAQ}QAQ。好吧，说这个的原因，只是为了引出我们关于蝴蝶操作的定义：

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

在了解完蝴蝶操作之后，我们首先考虑按照递归的思路，将FFTFFTFFT的分治流程刻画一下：

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

step  1step \ \ 1step  1成对地取出儿子节点，用蝴蝶操作计算出其DFTDFTDFT。 step  2step \ \ 2step  2用这一步的DFTDFTDFT替换之前的； step  3step \ \ 3step  3直到我们迭代到根节点为止，否则返回step  1step \ \ 1step  1

而反向迭代似乎有规律可循。我们发现只要我们用迭代的过程模拟出回溯的过程即可。那么思路便有了：三层forforfor，先枚举区间长度（1，2，4，8……），第二层枚举每个区间的起点（其实只有两段），第三层负责遍历每段区间（其实也只有两段ORZORZORZ），然后cppcppcpp是这样子的：

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 ;
            }
        }
    }

嗯，好像……海星？哈，思维不严谨的我们漏了一个地方：我们并不可以得到迄今为止的结构体参数JJJ（即这一段多项式）每一位的值。因为现在编号为111的是原序列中编号为444的，以此类推，我们需要重新排序……但实际上并不需要这么麻烦——

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

这个嘛……我们可以选择打个表观察： 原来的序号 0  1  2  3  4  5  6  70 \ \ 1 \ \ 2 \ \ 3 \ \ 4 \ \ 5 \ \ 6 \ \ 70  1  2  3  4  5  6  7 现在的序号 0  4  2  6  1  5  3  70 \ \ 4 \ \ 2 \ \ 6 \ \ 1 \ \ 5 \ \ 3 \ \ 70  4  2  6  1  5  3  7 原来的二进制表示 000  001  010  011  100  101  110  111000 \ \ 001 \ \ 010 \ \ 011 \ \ 100 \ \ 101 \ \ 110 \ \ 111000  001  010  011  100  101  110  111 现在的二进制表示 000  100  010  110  100  101  011  111000 \ \ 100 \ \ 010 \ \ 110 \ \ 100 \ \ 101 \ \ 011 \ \ 111000  100  010  110  100  101  011  111

诶，二进制好像是反序的嗷~~这便是我们的最后一个tricktricktrick，蝴蝶定理。

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

嗯，然后我们可以将这个反序存在一个数组里面。类似这样求出来：

    for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;

那么我们可以看到，这就简化了很多冗余的步骤，并让我们脱离递归的大常数。真开森啊

最后附迭代版的代码(我写的常数好像有点儿大QAQQAQQAQ)

#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_step\mathcal{Last\_step}Last_step·感悟

嗯……怎么说呢，现在看这个算法，真是简单的一匹啊……代码这么短这么容易背过。但是当时理解起来却花了很大心思呢！这篇博客我写了整整三天qwqqwqqwq，由于要培训和考试，所以拖拖沓沓地写了三天，一边写一边感叹自己理解的简直太浅显了，每一个证明、每一个引理、甚至每一个符号，我都需要去和其他DALAODALAODALAO比对审核、或者缠着rqyrqyrqy问个没完；每次一讲到原理，我都发现自己原来并不理解那些，于是不得不推倒重来。我希望你们可以经历这些，但同时更希望你们不经历这些。这篇博客会持续更新，补充语意不明、证明难以理解的地方。欢迎hackhackhack哟！

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

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

·rvaluervaluervalue大佬的博客[LinkLinkLink]【与《算法导论》结合起来食用效果更佳

·rqyrqyrqy（我缠着他问了好多关键点，大多是原理部分的）

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

posted @ 2019-04-01 20:33  γひん0ΖΖƦ  阅读(306)  评论(0)  编辑  收藏  举报