zoukankan      html  css  js  c++  java
  • 转自 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) 咱们可以由 欧拉公式(这个公式怎么证?咱不知道....貌似有本叫费恩曼物理学讲义的有讲过一种推导过程...)
        欧拉公式:
        【算法】【快速FFT】 - z55250825 - z55250825

        咱们可以得到 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,然后用蝴蝶操作做即可。
          
         存储的时候使用了一个辅助数组,这里只记了公共部分的。方便查找。
  • 相关阅读:
    C# 实现 Aop [Emit动态生成代理类方式]
    分享一些最近在看的电子书
    Can't connect to your phone. Disconnect it, restart it, then try connecting again
    07,Windows Phone后台代理
    .NET 性能测试工具 性能计数器
    windows 8 metro 开发学习资源链接
    08,Windows Phone 本地存储
    06,Windows Phone 8程序的生命周期
    .NET 性能测试工具 事件跟踪器(ETW)
    LitJSONjson 和net 的完美组合用法
  • 原文地址:https://www.cnblogs.com/zyfzyf/p/3857982.html
Copyright © 2011-2022 走看看