2017-03-12

快速傅立叶变换及逆变换

本文是对FFT及iFFT的简要介绍, 包括原理的证明和代码实现 (不包括应用). 基层的原理我也不清楚啦......比如复数次幂到底是怎么回事, 比如范德蒙德矩阵的行列式怎么算 (不会线代TAT), 但是这里给出的证明在假定这些事实的基础上是正确的.

个人认为本文的优点在于比较简洁, 最后给毛爷论文里提到的合并DFT的技巧给了个简单的证明.

参考文献: <算法导论>, 毛爷爷 <再探快速傅立叶变换>, Picks的博客, vfk <炫酷反演魔术>, uoj 34 上的优秀代码.

前置技能: 复数的加减乘法, 多项式的基本概念, 矩阵的基本概念. # 多项式多项式有两种表示法: 系数表示, 点值表示. 次数为n的多项式需要且仅需要曲线上的(n+1)个点 (坐标的值域为复数) 来确定, 这就是点值表示. 两种表示是等价的, 且可以相互转换: 系数 -求值-> 点值, 点值 -插值-> 系数. 系数表示的加法运算是的, 而朴素乘法运算是的:

即, 的系数等于a的前(i+1)项与b的前(i+1)项的卷积 (如果不足(i+1)项则补0).

点值表示中, 加法和乘法都是的, 因为只需要保持横坐标不变, 纵坐标相加相乘.

于是, 如果要把系数表示下的两个多项式相乘, 得到乘积的系数表示法, 可以这样做: 将系数表示转换为点值表示, 求点积, 再转换回去. 如果两次变换的时间复杂度低于, 那么这个算法优于朴素的乘法. 如果朴素地求值, 再用拉格朗日公式暴力插值, 两步的时间复杂度均为, 不能满足我们的要求. 愉快的是, 有FFT和iFFT, 在的时间内完成求值和插值.

复数

简要地介绍单位复根和共轭复数的一点性质. ## 单位复根 n次单位复根是满足的复数, 根据代数基本定理, 它们恰好有n个. 根据某套我不熟悉的理论, 它们是.

称为主n次单位根, 记作.

周期性

这个名称是我yy的......如有更学术的称呼请告知.

消去引理

由定义可知, 设, 则

推论:

折半引理

对正偶数n, n次单位复根的平方的集合即为n/2次单位复根的集合.

由消去引理, . 根据周期性, 当取遍, 每个n/2次单位复根出现两次.

求和引理

等比数列求和公式同样适用于复数. 设整数不能被整除, 则,

共轭复数

称实部相等, 虚部相反的两个复数互为共轭, 记作. 某数取两次共轭等于它自己, 和的共轭等于共轭之和, 积的共轭等于共轭之积.

后面会用到这样一个式子:

FFT

快速傅立叶变换 (FFT) 是离散傅立叶变换 (DFT) 的一种高效实现. DFT在其他领域的应用不是很了解, 只知道OI方面可以快速求卷积.

设是偶数, , , , 则

要往代入个值. 为了有效地减小问题规模, 同时又具有递归的结构, 和作为次数界为的多项式, 最好只代入个值.

只要是个不同的值即可. 这里选择代入, 因为由折半引理, 次单位复根的平方的集合即为次单位复根的集合. 设,

边界情况是, 此时是常数, 返回即可.

注意, 除了边界情况, 其他的都必须为偶数, 即n为2的自然数次幂.

由于次数界为m, n的两个多项式A, B相乘, 乘积的次数界为(n+m-1), 所以A, B都必须通过补0的方式扩充到长度l, l为不小于(n+m-1)的最小的2的幂.

过程FFT接受一个次数界为n (n是2的幂) 的多项式的系数表示, 返回代入后求得的值. 所有运算在复数域下进行. 其递归实现的伪代码如下:

FFT(a, n)
    if n = 1
	    y[0] = a[0]
	    return y
    for k=0 to n/2-1
        a_0[k] = a[2k]
        a_1[k] = a[2k+1]
    y_0 = FFT(a_0, n/2)
    y_1 = FFT(a_1, n/2)
    w = 1
    for k=0 to n/2-1
        t = w y_1[k]
        y[k] = y_0[k] + t
        y[n/2+k] = y_0[k] - t
        w = w w_n
    return y

其中, w称为旋转因子, 这几行称为蝴蝶操作:

w = 1
for k=0 to n/2-1
    t = w y_1[k]
    y[k] = y_0[k] + t
    y[n/2+k] = y_0[k] - t
    w = w w_n

根据主定理, FFT的时间复杂度为.

高效 (迭代) 实现

观察递归FFT怎样划分系数: fft recursion tree

设初始调用为第0层, 则第i层调用将初始编号倒数第i位为0的系数划分到左子树, 倒数第i位为1的系数划分到右子树. 这是一个递归版的基数排序. 所以, 最终得到位逆序置换: 000 -> 000 001 -> 100 010 -> 010 011 -> 110 100 -> 001 101 -> 101 110 -> 011 111 -> 111

可以地暴力计算, 也可以打个表均摊地计算.

重写FFT过程, 先把各个系数摆到递归树上叶结点的位置, 再向上合并.

FFT(a, n)
    for k=0 to n-1
        A[i] = a[rev[i]]
    for (m=1; m<n; m=2m)
        for (j=0; j<n; j=j+2m)
            w = 1
            for k=0 to n/2-1
                t = w A[j+m+k]
                A[j+m+k] = A[j+k] - t
                A[j+k] = A[j+k] + t
                w = w w_(2m)
    return A

m是待合并的两棵子树的大小, j是左子树在A中的位置.

iFFT

把DFT写成矩阵乘积的形式,

是由的适当次幂填充的范德蒙德矩阵, 即.

则. 可以验证.

当时, 上式=1. 当时, 或, 根据求和引理, 上式=0.

所以, 把看作多项式的系数, 将代入求值, 再把结果除以n即得iDFT.

实现上, 把FFT小小地修改一下即可:

FFT(a, n, d)
    for k=0 to n-1
        A[i] = a[rev[i]]
    for (m=1; m<n; m=2m)
        for (j=0; j<n; j=j+2m)
            w = 1
            for k=0 to n/2-1
                t = w A[j+m+k]
                A[j+m+k] = A[j+k] - t
                A[j+k] = A[j+k] + t
                w = w w_(2m)^d
    if d = -1
        for k=0 to n-1
            A[i] = A[i] / n
    return A

循环卷积

求和引理实际上是说:

把它代入循环卷积的表达式:

令上式即为

这给出了(I)DFT的另一种看待方式.

如果我们用FFT做多项式乘法的时候不把0补齐, 得到的错误答案其实就是循环卷积.

多维 (循环) 卷积

多维卷积可以看作多元多项式的求值和插值. 设每一维的长度为 . 固定其他维, 分别对每一维做长度为的DFT即可 (对DFT做DFT). 逆变换, 循环卷积同理.

合并DFT

把两个系数表示下的多项式相乘需要多少次DFT和iDFT? 通过两次DFT把多项式转换为点值表示, 做乘法, 再iDFT回去.

如果两个多项式的系数均为实数, 可以优化.

方法一

设. 平方, 得将DFT, 平方, 再iDFT, 实部即, 虚部即.

经过实验, 这个方法精度不高.

方法二

看到不少神犇这样写, 但不知道原理. 问了下Claris爷, 原来从毛爷的论文里可以找到.

我觉得使用共轭复数的性质会使推导更简单一些, 原文是用单位根的三角表示, 利用sin, cos的奇偶性来搞.

设.

所以, 只需要对DFT, 就能得到, 的点值表示, 乘起来再iDFT回去即可.

模板

由于Hexo不支持多个同级分类, 为了方便检索, 就分两篇文章了: [uoj 34] 多项式乘法