首页 > 代码库 > 转自 z55250825 的几篇关于FFT的博文（一）
转自 z55250825 的几篇关于FFT的博文（一）
关于FFT，咱们都会迫不及待地 @ .....(大雾)(貌似被玩坏了...)
.....0.0学习FFT前先orz FFT君。
首先先是更详细的链接（手写版题解点赞0v0） FFT的资料
其实众所周知的最详细的算法解释在《算法导论》上...然后咱就是边看着那个边码理解的...
首先来看看多项式乘法和快速FFT的关系，然后咱们再来看能否聊到卷积什么的东西...
其实觉得还是去看算法导论最好。
【一.多项式及其表达方式。】
首先什么是多项式额.....实际上就是一个类似这个的东西。
A(x)=∑aj*x^j (0<=j<=n-1)
对于 这样一个 式子咱们就称它是多项式，这里的 aj 便是系数，然后x 是一个变量，而满足 aj<>0的最大的j是这个多项式的次数，比如 f(x)=1+x+2x^x的次数就是3（它最大的不等于0的ak是k=3的时候a3=2）
然后考虑多项式乘法，咱们先来看看多项式乘法有什么表示方法。
如果没学快速fft之前，咱只会用这样一个式子 A(x)=∑aj*x^j 来表示 一个多项式，然后这种方法叫做系数表达法。咱们还有另外一种表达方法，叫做 点值表达法。即 咱们在次数为N的 A(x) 任意取N个点xi，求出 yi = A(xi)，然后咱们就可以用N个数对 （xi,yi）来表示 这个多项式，比如 A(x)=1+x，咱们可以用两个数 (0,1) (1,2) 来表示这个多项式。
为什么可以这样呢？只要咱们可以通过这个N个点表示出这个多项式就可以了。然后最方便的方法自然是设 N个系数为N个未知数，然后可以对N个点列出N个一次的方程，然后解方程即可，由于N个未知数，N个方程，这是显然可以解出来的。所以这种表示法就是正确的。另外实际上这个就是拉格朗日插值定理可以直接构造出来的...最后还有一种更严谨的证明方法见算导的那个矩阵的证明...咱线性代数弱...各种看不懂...但是上面两个比较不严谨的证明的方法应该可以了吧？
上面咱们证明了点值表达式和系数表达式实际上是一样的，由系数表达可以通过枚举N个点算值得到点值表达，然后由点值表达式又可以通过解方程或者插值得到系数表达。它们之间的相互转化实际上就是优化多项式乘法的关键。
实际上多项式乘法，即 计算两个n次多项式乘积 A(x)*B(x) 的结果，如果从系数表达考虑的话，实际上是无论如何也得不到一个优化的，复杂度就是O(n^2) 但是假设咱们把角度换到 点值表达式 的话，会发现只需要O(n)的，首先注意咱们这里的红字的结果指的是，得到最终乘出来的多项式的系数表达或者点值表达即可。这个用点值表示的话，咱们会发现，假设咱们知道了 A(x)的点值表达 (x0,y0),.....(xn-1,yn-1)和 B(x)的点值表达 (x0,y0‘),(x1,y1‘)....(xn-1,yn-1‘)的话，咱们可以迅速知道 结果多项式 C(x)的点值表达式 是 (x0,y0*y0‘) (x1,y1*y1‘)... (xn-1,yn-1*yn-1‘)。这个就是O(n)的，也就是说咱们对于多项式乘法，可以在O(n)的时间从两个次数为n的多项式的点值表达得到 相乘后的多项式的点值表达，那么咱们的任务实际上就是，考虑如何在O(n^2)以下的时间由 点值表达 变成 系数表达 或者由 系数表达 变成 点值表达。
【二.一些复数的证明什么的...】
首先考虑 咱们知道 一个多项式的系数表达 ，如何迅速变成 点值表达，如果直接枚举 N个点死算的话，咱们会发现，算一个需要O(n)的时间，算n个则是O(n^2)的时间，这个不优化的话显然还是比不上直接暴力。然而咱们的确有O(nlogn)的方法来做这件事情。
先得知道，这里的n个点任意的都是可以的，所以咱们的任务就是找到n个好算的点以便可以O(nlogn)做出来。然后还真有这样的点，那就是 e^(2kπi/n) (0<=k<=n-1)，这个的表示不是特别舒服（尤其是在一大堆这个东西的时候，目测很壮观0.0...），所以咱们就 用 w(k,n)表示 e^(2kπi/n)
这里先假设 n 是 2的幂，实际上不是的也可以多算一点让他变成2的幂。
咱们对于 次数为n的方程，咱们就求 e^(2kπi/n)（0<=k<=n-1）这n个点的值，其实就是求 w(k,n)的值。为什么用这个求会使复杂度变成O(nlogn)呢？
咱们来看看一些性质：
1）消去定理 ： 对于任何证书 n>=0,k>=0,以及 d>0，咱们有以下的性质：
w(dk,dn)=w(k,n)
证明：
根据定义 w(k,n)=e^(2kπi/n)
w(dk,dn)=e^(2dkπi/dn)=e^(2kπi/n)=w(k,n)
然后由该定理可以导出一个较为关键的东西：
(设n为偶数)w(n/2,n)=w(1,2)=e^(2πi/2)=e^(πi)
e^(πi) 咱们可以由 欧拉公式（这个公式怎么证？咱不知道....貌似有本叫费恩曼物理学讲义的有讲过一种推导过程...）
欧拉公式：
咱们可以得到 e^(iπ)=cosπ+isinπ=-1
即：
w(n/2,n)=-1
2)折半定理：如果n>0为偶数，那么n个n次单位的复数根（这个就是w(k,n)们(0<=k<=n-1)）的平方的集合就是n/2个n/2单位的复数根（即w(k,n/2)(0<=k<=n/2-1)）的集合。
证明： 首先咱们证明 w(k,n)^2=w(k+n/2,n)^2
首先证明 w(i+j,n)=w(i,n)*w(j,n)，这个为什么对的自己化成 e^blabla 的形式就可以了。
这个实际上就是 w(k+n/2,n)^2=(w(k,n)*w(n/2,n))^2=(w(k,n)*(-1))^2 (w(n/2,n)=-1是上面1)的结果)
=w(k,n)^2
这个证明完了，然后咱们再证明w(k,n)^2=w(k,n/2)即可，这个其实也是消去定理的直接运用。
咱们有 w(k,n)^2=w(2k,n)=w(k,n/2)
即 w(k,n)和w(k+n/2,n)的平方都对应 w(k,n/2)
所以折半定理是正确的。
其实还有一个求和引理的..貌似没看到用所以就先跳过，待会儿如果有需要用了就补上。
【三.DFT和FFT】
实际上咱们这个问题就转化成了求 A（x）在 w(k,n) (0<=k<=n-1)的值，总共n个，咱们把结果 yi 称作 系数 ai 的离散傅里叶变换，而离散傅里叶变换又叫DFT。（以前一直以为是FFT的简化版本的...现在看起来应该不是的）
而FFT即在O(nlogn)求出DFT的算法。它是利用了复数根的特殊性来求的。
首先确定咱们现在的问题，是求 A（x）取 w(k,n) 的n个值。
这里先考虑特殊情况，即n为2的幂的情况。
即求 A(x)=a0+a1*x+a2*x^2...+an-1*x^(n-1)在 w(k,n)的n个值。
咱们考虑把 A(x)拆开。
咱们按照 ai的i的 奇偶性将他拆成下面两个：
B(x)=a0+a2*x^2+a4*x^4...+a(n-2)*x^(n-2)
C(x)=a1*x+a3*x^3+a5*x^5...+a(n-1)*x^(n-1)
即 A(x)=B(x)+C(x)
然后咱们看到C(x)实际上还有化简的余地额...即
C(x)=x(a1+a3*x^2+a5*x^4...+a(n-1)*x^(n-2))=xD(x)
即 A(x)=B(x)+xD(x)
要求A(x)，实际上求B(x)和D(x)即可。
先考虑求B(x)，如果咱们把 B(x)中的x^2看成 y会发生什么？
咱们会发现问题变成了
求 B(y)=a0+a2*y+a4*y^2....+a(n-2)*y^(n/2-1) 这个式子代入 y=w(k,n)^2 的n个答案。
咱们发现了什么？ 咱们把w(k,n)^2变成w(k,n/2)，咱们发现，要算的实际上变成了一半了!(折半定理)，即咱们把原问题转化成了这样一个问题，求 B(y)=a0+a2*y+a4*y^2...在 y=w(k,n/2) 的 n/2个答案。
这个的规模貌似和原来的问题相比小了一半额....但是还是这个问题额...
然后咱们递归下去做...一直到 B(y)次数为1的时候（就是a0）的时候，咱们直接返回答案a0即可。
这个就是分治的!复杂度O（nlogn）。
考虑怎么把返回来的答案组合成这个区间的答案。
之前咱们考虑求 A(x)=B(x)+xD(x)，然后B(x)又可以递归求，咱们就可以得到B(x)的n/2个解答，然后咱们把它们复制一遍就是B(x)的n个解，然后咱们的D(x)也是类似的，咱们也复制一遍，再每一个乘以对应的x，在分治的时候复杂度是O(n)的，就可以得到自己的n个解答，由主定理可以得到这样的复杂度就是O(nlogn)。然后 由于 w(k+n/2,n)=w(k,n)*w(n/2,n)=-w(k,n)，这个可以只算一半的w即可。剩余一半加负数。
再附上 算导上的 伪代码把：
Recursive-FFT(A)
n=A.length {n为A的长度}
if n==1 return A
wn=e^(2πi/n)
w=1
B=(A0,A2,...,An-2) {B赋值A的a0,a2...a(n-2)}
D=(A1,A3,...,An-1) {D类似}
yB=Recursive-FFT(B) {递归求B}
yD=Recursive-FFT(D) {递归求D}
for k=0 to n/2-1
Y[k]=yB[k]+wyD[k]
Y[k]+n/2=yB[k]-wyD[k]
w=w*wn
return Y {Y是n个点的答案}
接下来考虑插值...也就是把 点值转化为 系数。
注意咱们之前 系数 转化为 点权，实际上是求 n个 ∑A(w(k,n))的和，在O(nlogn)的时间内。
而咱们可以证明（这个证明咱还不会...）
aj=∑ykw(-kj,n)/n
然后这个的话实际上 和 yj=∑ak*w(kj,n) 是相互逆运算的（乃把这个式子的a看成y,y看成a就可以发现了）
所以这个也可以用FFT转换，咱们先把求出来的y，用O(n)求出乘后的点值表达式，然后把表达式的y看做a,再做一遍转换的FFT,转换的方法就是 对于 伪代码 第四行的那个 wn 一开始初始化为指数为负数的，然后求一遍FFT之后再给求出来的每一个y 除以n，就可以得到多项式乘法后的式子的 a了。这个的复杂度还是O(nlogn)的，所以总的复杂度就是O(nlogn)
【四 FFT的优化及迭代实现啥的】
首先有一个优化叫蝴蝶操作，名字很漂亮但是还是比较简单的...注意到12~13行的 wyD[k] 算了两次，实际上咱们可以只算一次的，咱们用个变量存储它，然后再算。
0.0...然后 vfk 的代码是递归的貌似...不知道该不该看迭代的饿...
vfk的fft
然后迭代的算法其实感觉...比较类似zkw自底向上的做法...递归实际上可以看做一棵完全二叉树（注意咱们已经假设了n是2次幂了）,然后咱们可以类似zkw直接自底向上的做。然后这个咱们可以用那种蝴蝶操作做，问题是咱们得确定叶子节点的顺序。实际上这个的确定的顺序是按某个点是0还是1来的，第一位是0的点全去了左子树，1的去了右子树，咱们按照这个方法按普通DFS遍历一遍数O(n)就可以求出来了。
再给一个算导上的伪代码：
Iterative-FFT(a)
Bit-Reverse-Copy(a,A){算出叶子的顺序存在A中}
n=a.length
for s=1 to lgn
m=2^s
wm=e^2π/m
for k=0 to n-1
w=1
for j=0 to m/2-1
t=wA[k+j+m/2]
u=A[k+j]
A[k+j]=u+t
A[k+j+m/2]=u-t
w=w*wm
Return A
[五 入门题目】
哪一个OJ的高精度乘法其实都可以额...
[六 卷积咱也不知道是什么...]
只是发现了一个比较好的讲卷积的博文...总算不用纠结连续函数的卷积啥的了0v0
好的卷积
然后貌似看懂了额...
咱们OI一般都是讲离散化的东西...所以这里讲的卷积也是 离散化的。
咱们设有两个序列 f[n],g[n]，那么卷积 s[k]=∑f[i]*g[k-i]，很容易证明多项式的每一位的系数实际上就是那一位的卷积。卷积实际上也是一个函数，对于某一个数 k，它的值为 s[k]=∑f[i]*g[k-i]，而 多项式乘法实际上就是它们的系数向量的卷积（向量也可以看做一个函数把），然后咱们的快速傅里叶变换（FFT）算法，实际上就是可以在O(nlogn)求出两个函数的卷积。令 A x B 表示 A与B的卷积，即
A x B = DFT^-1(DFT(A)*DFT(B))
【七 模板】
（以下排名不分先后，按字典序来）
Evil君的模板
Error君的模板
vfk君的模板
【八 实现上的细节】
感觉实现上还是需要理解一下的饿...
首先要知道，咱们只在最后递归到尽头的时候使用了一下a[]来计算，在递归返回的时候实际上它们已经变成了y[]了。
首先是看w[]的变化，每一次往下除以2，实际上可以看做另一部分乘以2，所以咱们预处理出 0~n-1的w是机智的。除以2的实际上就是 1<<t，然后咱们枚举当前能够枚举的数量，乘以这个1<<t就是咱们的w[]，也就是说咱们可以用 i <<t来得到w[]。i的范围为当前能够循环的w的一半，也就是 0~n>>(t+1)-1。即咱们使用w[i shl t]即可。（还有另一半，为负数）
咱们在递归的时候记录一下咱们所得到的a[]的具体信息，然后再来算，首先就是s，s记录的实际上是当前所涉及的所有a[]的公共二进制部分，而以上的可以靠自己枚举了，也就是 0~n>>(t+1)-1（这个是默认剩余部分最低位为0的）它们所对应的就是剩余部分为1的，即再加上个 1<<t即可，分到左边的为 p=i<<(t+1)+s，而其对应的系数应该就是 p+1<<t，然后用蝴蝶操作做即可。
存储的时候使用了一个辅助数组，这里只记了公共部分的。方便查找。
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。