快速傅里叶变换（FFT）（下）

关于学习FFT算法的资料个人最推荐的还是算法导论上的第30章(第三版), 多项式与快速傅里叶变换, 基础知识都讲得很全面。

FFT算法基本概念：

FFT(Fast Fourier Transformation)即快速傅里叶变换, 是离散傅里叶变换的加速算法, 可以在O(nlogn)的时间里完成DFT, 利用相似性也可以在同样复杂度的时间里完成逆DFT。
DFT(Discrete Fourier Transform) 即离散傅里叶变换, 这里主要就是多项式的系数向量转换成点值表示的过程。
在ACM-ICPC竞赛中, FFT算法常被用来为多项式乘法加速, 即在O(nlogn)复杂度内完成多项式乘法, 当然实际应用不仅仅限于这些, 时常出现需要构造多项式相乘来进行计数的问题, 也需要用FFT算法来解决, 相关的几个问题在本文中也会提及。

FFT算法需要的基础数学知识：
多项式相关：
多项式相关的定义：一个以x为变量的多项式定义在代数域F上, 将函数A(x)表示为形式和：

称

为系数, 所有系数属于代数域F, 如果一个多项式的最高非零系数是

, 那么成这个多项式的次数是k, 记做degree(A)=k, 任何一个严格大于一个多项式次数的整数都是该多项式的次数界.

关于多项式的加法和乘法相信看到这篇博客的读者都会最基本的中学的算法, 在计算的时候, 如果采用传统的中学的计算方法, 多项式加法的时间复杂度是O(n), 乘法的时间复杂度是O(n^2) (n是两个多项式A和B的次数)。

多项式的表示:

在平常的学习中, 最常见的是多项式的系数表达方式, 次数界为n的

多项式的系数表达为一个由系数组成的向量

。 但是多项式还有一个比较常用的表示方法, 即多项式的点值表达方式。

一个次数界为n的多项式

的点值表达就是一个由n个点对组成的集合

使得对任意的整数k=0,1,..n−1，

各不相同, 且

=A(

)。

关于点值表达的正确性证明：对于任意n个点值对组成的集合

, 如果存在一个次数界为n的多项式A(x)过这n个点, 那么：

最左边这个n*n的矩阵称为范德蒙德矩阵, 可以用数学归纳法证明它的行列式值为：

。当

两两不相同时明显这个行列式的值不为0, 该矩阵可逆, 于是存在唯一解, 所以多项式的点值表达是合理的。

相应的通过n个点的坐标直接确定多项式各个系数的值的方法是存在的, 感兴趣的读者可以查询拉格朗日公式的相关资料, 利用拉格朗日插值公式可以在O(n^2)的时间复杂度内得到多项式的系数表达, 这个也是算法导论中的一个习题

点值表达方式下多项式的乘法, 不难发现, 如果多项式A(x),B(x)的点值表示分别是：

和

那么如果多项式C(x)=A(x)B(x), 那么C(x)的点值表达是：

卷积：

对于两个多项式的系数向量

和

, 两个多项式相乘得到的多项式的系数向量

满足

，称系数向量c是输入向量a和b的卷积, 记作c=a⊗b。

简单的多项式乘法的计算方法中, 每一个多项式的系数都通过系数表示方式下卷积的方式来进行计算, 时间复杂度是O(n^2), 但是FFT是先将多项式的从系数表示法转换成点值表示法(可以在O(nlogn)的时间复杂度下完成, 也就是加速的DFT变换, 然后在点值表示法下进行乘积计算, 在O(n)的时间复杂度内得到结果的点值表示法, 然后进行逆DFT变换, 在O(n*logn)的时间复杂度下完成逆DFT变换得到系数表示法。

而要理解DFT, 则需要一定复数上的数学知识：

复数相关的基础知识：

单位复数根：n次单位复数根指的是满足

的所有复数ω, n次单位复数根刚好有n个, 他们是

, 其中i是复数单位, k=0,1,2...n−1, 在复平面上这n个根均匀的分布在半径为1的圆上, 关于复数指数的定义如下：

其中

倍称为主n次单位根(这个定义好像接下来没用到)

关于复数根的几个定理和引理：消去引理： 对任何整数

有

证明： 

一个推论： 对任意偶数 n > 0 有

证明：设n = 2*k那么 

折半引理：如果n > 0是偶数, 那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合

证明：实际上这个引理就是证明了

折半引理对于采用分治对多项式系数向点值表达的转换有很大作用, 保证了递归的子问题是原问题规模的一半。

求和引理：对任意整数n≥1和不能被n整除的非负整数k, 有

这个问题通过等比数列求和公式就可以得到：

DFT与FFT, 以及逆DFT：

DFT:

在DFT变换中, 希望计算多项式A(x)在复数根

处的值, 也就是求

称向量y=(y0,y1,...,yn−1)是系数向量a=(a0,a1,...,an−1)的离散傅里叶变换, 记为

FFT:

直接计算DFT的复杂度是O(n^2)， 而利用复数根的特殊性质的话, 可以在O(n*logn)的时间内完成, 这个方法就是FFT方法, 在FFT方法中采用分治策略来进行操作, 主要利用了消去引理之后的那个推论。

在FFT的策略中, 多项式的次数是2的整数次幂, 不足的话再前面补0, 每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分：

于是就有了

那么我们如果能求出次数界是n/2的多项式

和

在n个n次单位复数根的平方处的取值就可以了, 即在

处的值, 那么根据折半引理, 这n个数其实只有n/2个不同的值, 也就是说, 对于每次分出的两个次数界n/2的多项式, 只需要求出其n/2个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和的复杂度内解决, 那么这样递归问题的复杂度将会是O(nlogn)的, 用a=(a0,a1,...,an−1)表示系数向量, y=(y0,y1,...,yn−1)表示离散变换之后的向量, 这里给出将算导上的代码翻译出来的C++代码实现(以解决算法第三版大论第三十章的一个习题, 求(0, 1, 2, 3)的DFT为例)：

#include<bits/stdc++.h>
using namespace std;
const double eps(1e-8);
typedef long long lint;

const double PI = acos(-1.0);
/*
 * 这是一个递归实现FFT的测试, 测试习题中求DFT(0, 1, 2, 3)
 */

struct Complex
{
    double real, image;
    Complex(double _real, double _image)
    {
        real = _real;
        image = _image;
    }
    Complex(){}
};

Complex operator + (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real + c2.real, c1.image + c2.image);
}

Complex operator - (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real - c2.real, c1.image - c2.image);
}

Complex operator * (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);
}

Complex* RecursiveFFT(Complex a[], int n)//n表示向量a的维数
{
    if(n == 1)
        return a;
    Complex wn = Complex(cos(2*PI/n), sin(2*PI/n));
    Complex w = Complex(1, 0);
    Complex* a0 = new Complex[n >> 1];
    Complex* a1 = new Complex[n >> 1];
    for(int i = 0; i < n; i++)
        if(i & 1) a1[(i - 1) >> 1] = a[i];
        else a0[i >> 1] = a[i];
    Complex *y0, *y1;
    y0 = RecursiveFFT(a0, n >> 1);
    y1 = RecursiveFFT(a1, n >> 1);
    Complex* y = new Complex[n];
    for(int k = 0; k < (n >> 1); k++)
    {
        y[k] = y0[k] + w*y1[k];
        y[k + (n >> 1)] = y0[k] - w*y1[k];
        w = w*wn;
    }
    return y;
}

int main()
{
    Complex* a = new Complex[10];
    a[0] = Complex(0, 0);
    a[1] = Complex(1, 0);
    a[2] = Complex(2, 0);
    a[3] = Complex(3, 0);
    Complex* ans = new Complex[10];
    ans = RecursiveFFT(a, 4);
    for(int i = 0; i < 4; i++)
        cout<<"("<<ans[i].real<<")"<<"+"<<"("<<ans[i].image<<")"<<"i"<<endl;

    return 0;
}

卷积定理：

对任意两个长度为n的向量a和b, 其中n是2的幂, 有

其中向量a和b用0填充, 使其长度达到2n, 并用⋅表示两个2n个元素组成向量的点乘(也就是每一维上的数相乘)
这个式子实际上就是多项式的系数表达在乘法时进行的卷积运算得到的结果, 等同于通过将其系数进行DFT变换变成点值表达之后相乘再换回来的过程

关于FFT算法的迭代实现：

在递归实现DFT过程的FFT算法中, 我们每次将系数向量a分成两个部分利用折半引理来降低计算的规模, 可以发现在每次分组当中他们满足这样一个完全二叉树的分组(n是2的幂)：

FFT递归系数分组

通过上图的流程可以看出, 最后一层的子节点下标的顺序实际上就是其下标转换成二进制串的倒序的字符串按照字典序排列的顺序。

比如a0~a7得到的序列a0, a4, a2, a6, a1, a5, a3, a7下标的二进制是000, 100, 010, 110, 001, 101, 011, 111对应的串的倒序是000, 001, 010, 011, 100, 101, 110, 111这个倒序的二进制刚好是0, 1, 2, 3, 4, 5, 6, 7的二进制表示, 于是我们可以在O(nlogn)的复杂度内得到做下面一层的下标顺序, 然后可以根据子节点的结果向上迭代得到父亲结点的值, 这样计算的话直接避免了递归, 如果直接递归的话在一些OJ上可能会造成爆栈的错误, 所以还是采用迭代的方式进行比较好。

关于迭代形式的FFT算法, C++代码实现如下(同样以求(0, 1, 2, 3)的DFT变换为例)：

这段代码的话同时也进行了逆DFT, DFT和逆DFT的过程相似, 加上一个标记判断当前执行的是哪一种就行了。

#include<bits/stdc++.h>
using namespace std;
const double eps(1e-8);
typedef long long lint;

const double PI = acos(-1.0);

struct Complex
{
    double real, image;
    Complex(double _real, double _image)
    {
        real = _real;
        image = _image;
    }
    Complex(){}
};

Complex operator + (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real + c2.real, c1.image + c2.image);
}

Complex operator - (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real - c2.real, c1.image - c2.image);
}

Complex operator * (const Complex &c1, const Complex &c2)
{
    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);
}

int rev(int id, int len)
{
    int ret = 0;
    for(int i = 0; (1 << i) < len; i++)
    {
        ret <<= 1;
        if(id & (1 << i)) ret |= 1;
    }
    return ret;
}

//当DFT= 1时是DFT, DFT = -1则是逆DFT
Complex* IterativeFFT(Complex* a, int len, int DFT)//对长度为len(2的幂)的数组进行DFT变换
{
    Complex* A = new Complex[len];//用A数组存储数组a分组之后新的顺序
    for(int i = 0; i < len; i++)
        A[rev(i, len)] = a[i];
    for(int s = 1; (1 << s) <= len; s++)
    {
        int m = (1 << s);
        Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m));
        for(int k = 0; k < len; k += m)//这一层结点的包含数组元素个数都是(1 << s)
        {
            Complex w = Complex(1, 0);
            for(int j = 0; j < (m >> 1); j++)//折半引理, 根据两个子节点计算父亲节点
            {
                Complex t = w*A[k + j + (m >> 1)];
                Complex u = A[k + j];
                A[k + j] = u + t;
                A[k + j + (m >> 1)] = u - t;
                w = w*wm;
            }
        }
    }
    if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len;
    return A;
}

int main()
{
    Complex* a = new Complex[4];
    a[0] = Complex(0, 0);
    a[1] = Complex(1, 0);
    a[2] = Complex(2, 0);
    a[3] = Complex(3, 0);
    a = IterativeFFT(a, 4, 1);
    cout<<"----------After DFT----------"<<endl;
    for(int i = 0; i < 4; i++)
        printf("%.9f + (%.9f) i\n", a[i].real, a[i].image);
    cout<<"----------After DFT-1----------"<<endl;
    a = IterativeFFT(a, 4, -1);
    for(int i = 0; i < 4; i++)
        printf("%.9f + (%.9f) i\n",a[i].real, a[i].image);
    return 0;
}

通过上面这个代码的示例, FFT算法的实现基本没有什么问题了, 另外算法导论中的习题有一些很不错, 便于熟悉这一算法的很多细节, 这里就不一一提及了。

经过这些学习之后, 进行一些实战演练是很有必要的, 接下来是相关习题的练习部分。

HDU 1402 A*B Problem Plus 大整数乘法
HDU 4609 3-idiots FFT计数
UVALive 4671 K-neighbor Substrings FFT算字符串Hamming距离 
UVA 12298 Super Poker II FFT计数, long double
URAL 1996 Cipher Massage 3 FFT + KMP 
CodeChef COUNTARI Arithmetic Progressions FFT + 分块
ZOJ 3856 Goldbach FFT计数
UVALive 6886 Golf Bot FFT模板题
HDU 4093 Xavier is Learning to Count 容斥原理 + FFT,(2011年上海现场赛C题)

HDU 5751 BestCoder Round #84 Eades（线段树+FFT）

HDU 5730 2016多校1 Shell Necklace （CDQ分治+FFT）

2016 acm香港网络赛 A题 A+B Problem （FFT）

以后还有题目可以看我博客，可能有。最后感谢Ichimei对FFT的讲解。Orz....对理解FFT真的有很大帮助。