码迷,mamicode.com
首页 > 其他好文 > 详细

转自 z55250825 的几篇关于FFT的博文(一)

时间:2014-07-22 00:32:37      阅读:431      评论:0      收藏:0      [点我收藏+]

标签:style   blog   http   color   使用   os   

    关于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) 咱们可以由 欧拉公式(这个公式怎么证?咱不知道....貌似有本叫费恩曼物理学讲义的有讲过一种推导过程...)
    欧拉公式:
    bubuko.com,布布扣

    咱们可以得到 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 {nA的长度}
     if n==1 return A
     wn=e^(2πi/n)
     w=1
     B=(A0,A2,...,An-2) {B赋值Aa0,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 {Yn个点的答案}
   接下来考虑插值...也就是把 点值转化为 系数。
   注意咱们之前 系数 转化为 点权,实际上是求 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))

 

【七 模板】
   (以下排名不分先后,按字典序来)

 

【八 实现上的细节】
       感觉实现上还是需要理解一下的饿...
      首先要知道,咱们只在最后递归到尽头的时候使用了一下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,然后用蝴蝶操作做即可。
      
     存储的时候使用了一个辅助数组,这里只记了公共部分的。方便查找。

转自 z55250825 的几篇关于FFT的博文(一),布布扣,bubuko.com

转自 z55250825 的几篇关于FFT的博文(一)

标签:style   blog   http   color   使用   os   

原文地址:http://www.cnblogs.com/zyfzyf/p/3857982.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!