zoukankan      html  css  js  c++  java
  • 模板:快速傅里叶变换(FFT)

    参考:http://blog.csdn.net/f_zyj/article/details/76037583
    如果公式炸了请去我的csdn博客:http://blog.csdn.net/luyouqi233/article/details/79323568
    原文即是一篇很好的FFT入门博客,但是笔者打算为了日后的学习,则将原篇章的结构删改增添一下,如有思路上的雷同十分正常。
    “是时候打开FFT的大门了!”

    预备知识:

    1.至少知道基础数论与一定解三角形知识(大概是高中水平)。
    2.定义(i=sqrt{-1})
    3.引入复数(即形如(a+bi)(a,b均为实数)的数的集合)
    4.((cos heta+i imes sin heta)^k=cos(k heta)+i imes sin(k heta))
    5.显然我们对多项式FFT之后得到的答案不是我们想要的,那么这时候就需要反着用FFT把式子再变回去(本文记做IFFT)。

    这里证明一下第四条,用归纳法。
    显然当(k=1)时成立。
    (k)成立时,我们有:
    ((cos heta+i imes sin heta)^{k+1})
    (=(cos heta+i imes sin heta)^k imes (cos heta+i imes sin heta))
    (=(cos(k heta)+i imes sin(k heta)) imes (cos heta+i imes sin heta))
    (=cos(k heta)cos heta+i imes sin(k heta)cos heta+i imes cos(k heta)sin heta+i^2 imes sin(k heta) sin heta)
    (=cos(k heta)cos heta-sin(k heta) sin heta+i imes (sin(k heta)cos heta+cos(k heta)sin heta))
    (=cos((k+1) heta)+i imes sin((k+1) heta))
    得证。

    问题引入:

    (A(x)=sum_{i=0}^{n-1}a_ix^i,B(x)=sum_{i=0}^{n-1}b_ix^i),求(A(x) imes B(x))后的多项式系数。

    初探:

    显然我们有一个(O(n^2))的解法,但是实在是太慢了。
    考虑到一个(n-1)次多项式可以看做是定义在复数域上的函数,则我们一定可以找到n个点来唯一确定这个函数。
    当然我们也可以通过这些点来表示这个多项式。
    假设:
    (A(x))被表示为:(<(x_0,y_{a_0}),(x_1,y_{a_1}),ldots,(x_{2n-2},y_{a_{2n-2}})>)
    (B(x))被表示为:(<(x_0,y_{b_0}),(x_1,y_{b_1}),ldots,(x_{2n-2},y_{b_{2n-2}})>)
    显然(A(x) imes B(x))被表示为:(<(x_0,y_{a_0}y_{b_0}),(x_1,y_{a_1}y_{b_1}),ldots,(x_{2n-2},y_{a_{2n-2}}y_{b_{2n-2}})>)

    这里多取了点的原因在于(A(x) imes B(x))是一个(2n-2)次多项式,则至少要取(2n-1)个点才能保证正确。

    但是显然还是(O(n^2))的。

    再试:

    考虑设(A(x_i)=A_0(x_i^2)+x_iA_1(x_i^2)),其中:
    (A_0(x)=a_0+a_2x+a_4x^2+ldots+a_{n-2}x^{frac{n}{2}+1})
    (A_1(x)=a_1+a_3x+a_5x^2+ldots+a_{n-1}x^{frac{n}{2}+1})

    其实就是按照系数下标的奇偶性分类了一下。

    此时我们再令取点的(x)值为(<x_0,x_1,ldots,x_{frac{n}{2}-1},-x_0,-x_1,ldots,-x_{frac{n}{2}-1}>)
    我们发现把(x)平方后我们的取值瞬间缩小了一半,而原式唯一变化的就是(A_1(x))前的符号。
    看起来我们似乎找到了(O(nlogn))的可行方案。
    但是很可惜,这样优秀的(x)取值的性质只会保留一次,也就是说我们只是得到了一个(O(frac{n^2}{2}))
    如何才能每次将问题的规模缩小一半是我们的目标。

    插曲:

    有个人告诉你:不如试试(X_n=cosfrac{2pi}{n}+i imes sinfrac{2pi}{n})(0ldots n-1)次方作为(x)的取值。

    这块大家一直有个疑惑:这是怎么构造出来的啊?
    事实上傅里叶变换最早是应用于信号处理上的,傅里叶提出:任何连续周期信号可以由一组适当的正弦曲线组合而成。
    多项式可以看做非连续周期信号,然后通过各种奇妙的姿势让它逼近正弦曲线的组合形,详情可以看松松松WC2018的课件。
    “逼近”显然用到了微积分,不适合初学者,所以就直接跳过了。(其实我也不会……)
    (再多说一点吧,其实上面和下面的数学推理完全可以从物理层面理解,还是可以参考松松松WC2018的课件)

    继续:

    那么令取点的(x)值为(<X_n^0,X_n^1,ldots,X_n^{n-1}>)
    我们可知:
    ((X_n^{k})^2)

    (=X_n^{2k})

    (=cosfrac{2k imes 2pi}{n}+i imes sinfrac{2k imes 2pi}{n})

    (=cosfrac{2kpi}{frac{n}{2}}+i imes sinfrac{2kpi}{frac{n}{2}})

    (=X_{frac{n}{2}}^k)


    (X_n^{k})

    (=cosfrac{k imes 2pi}{n}+i imes sinfrac{k imes 2pi}{n})

    根据三角函数的周期性可知,(k)(n)取模显然不会对答案造成影响。
    于是我们有(X_n^{k}=X_n^{k\%n})

    那么显然对于(<(X_n^0)^2,(X_n^1)^2,ldots,(X_n^{n-1})^2>)

    它等效于(<X_{frac{n}{2}}^0,X_{frac{n}{2}}^1,ldots,X_{frac{n}{2}}^{frac{n}{2}-1},X_{frac{n}{2}}^0,X_{frac{n}{2}}^1,ldots,X_{frac{n}{2}}^{frac{n}{2}-1}>)

    我们好像看到了(O(nlogn))的曙光了。

    尾声:

    显然我们可以对(x)的取值折半,然后对于左右区间的(x)值递归下去即可。
    Q1:诶等等,“再试”里面的内容好像没有应用上啊……

    A1:那就转化一下,其实我们只需要求一个区间的(A_0(x))(A_1(x))值递归下去求(A(x))即可。
    也就是说其实我们是得到了:
    (<(A_0)_0,(A_0)_1,ldots,(A_0)_{frac{n}{2}-1},(A_1)_0,(A_1)_1,ldots,(A_1)_{frac{n}{2}-1}>)

    Q2:这好像是画蛇添足……

    A2:emmm……我说这个可以用于常数优化你信吗……
    显然(A(X_n^k)=(A_0)_{k\%frac{n}{2}}+X_n^k(A_1)_{k\%frac{n}{2}})

    取模是因为,不要忘了我们的取值是由两个一样的左右区间合并在一起的。

    那么我们得到了(<A_0,A_1,ldots,A_{n-1}>)

    (其中(A_k=A(X_n^k))

    我们好像把这个序列的长度减少了一半诶!那自然是快了二倍啊。

    不要忘了n要满足始终是2的倍数,所以n要取2的整数次幂,同时将没用的次幂的系数填成0。

    Q3:IFFT怎么做啊?

    A3:继续看下去……?

    补遗:

    略讲一下IFFT。
    显然我们可以把FFT的最初算法(也就是DFT)看做两个矩阵相乘。

    两个矩阵分别一个填((X_n^k)^m),一个填系数,可以上参考处原博客看矩阵。

    那么我们把第一个矩阵变成逆矩阵岂不是为IFFT?
    其实就是这样,并且事实上就是填(((X_n^{-k})^m)/n),具体证明过程看参考处原博客。
    剩下的做法就和FFT一样啦。

    谢幕:

    事实上我上述讲的内容其实没有多少用(滑稽。
    因为你理解半天也不如不理解知道怎么用然后默写下来。
    但是理解了更好背啊。

    例题:

    模板:
    HDU1402:A * B Problem Plus:http://www.cnblogs.com/luyouqi233/p/8448969.html

    应用:
    BZOJ3527:[ZJOI2014]力:http://www.cnblogs.com/luyouqi233/p/8452117.html

    +++++++++++++++++++++++++++++++++++++++++++

    +本文作者:luyouqi233。               +

    +欢迎访问我的博客:http://www.cnblogs.com/luyouqi233/+

    +++++++++++++++++++++++++++++++++++++++++++

  • 相关阅读:
    gc buffer busy/gcs log flush sync与log file sync
    给Oracle年轻的初学者的几点建议
    Android 编程下帧动画在 Activity 启动时自动运行的几种方式
    Android 编程下 Touch 事件的分发和消费机制
    Java 编程下 static 关键字
    Java 编程下 final 关键字
    Android 编程下模拟 HOME 键效果
    Why Are Thread.stop, Thread.suspend, Thread.resume and Runtime.runFinalizersOnExit Deprecated ?
    Extjs4 大型项目目录结构重构
    [转]SQLServer 2008 允许远程连接的配置方法
  • 原文地址:https://www.cnblogs.com/luyouqi233/p/8447569.html
Copyright © 2011-2022 走看看