浅谈FFT

Piwry · 2020-01-13 13:46:55 · 个人记录

-- 数论的模板一大堆，但做着做着就会发现所谓"模板"，就是如何套"模板"的经典例子（以及拓展理论和原理论的做法没有半毛钱关系）。但FFT从其出现到现在，其在OI中的重要性是母庸质疑的，一定要深入的学习并熟练的应用，否则就会难以应付多元的以FFT为起点的"模板"与题目。

• p.s.这次我写了整整4天qwq，中途码着码着又不满意还删了重来一次，所以文章多少有一些错位感就放过我吧。FFT真的和其他模板不同，本来以为是挺良心的东西结果这么抽搐，却把一个人的自学能力锻炼到了极致。下面，就请食用把。

• 进阶部分可以参考这篇日报文章。尤其是有关 DFT含义 和 卷积 部分，这篇更加准确。

目录

• 参考和推荐

• 正文：

1. FT, DFT

2. FFT

3. FFT的实现

4. 序列的操作

5. DFT与多项式乘法

6. FFT与NTT

7. *trick-打包FFT

• 后记

参考和推荐

推下花姐姐\color{#00bbbb}\texttt{皎月半洒花}的\color{#ffb464}\texttt{博文}，这是我的FFT启蒙QwQ（同时正因为是启蒙，这篇博文中的一些证明方式其实与DFT并无关系，它们也不会出现在这里；但不可否认的是这些方式是正确的，也是初学时易于理解的）；

还有这篇(《再探快速傅里叶变换》--毛啸)，建议在初步掌握算法后认真研读一遍（果然国家队都是神仙，这篇论文写得比其他教程要严谨清楚多了ORZ）；

最后是这篇，应该是大学教材ppt摘录，其中的分析方式是较现代化的。

正文

1. FT, DFT

傅里叶变换( FT )，是指将函数描述为三角函数的线性组合，这样就可以方便的对多个非线性函数或序列操作，通常如计算卷积；

离散傅里叶变换( DFT )，其实就是在转换前后都成离散形式的傅里叶变换，但其本质是和傅里叶变换一样的，其用也一样是为了方便对多个非线性序列的操作

.

DFT/IDFT公式：

其中的j指j^2=-1（为了和电学的"i"区分）；整个公式是从连续傅里叶变换推得的。

同时，式中的j,-j是可以互换的，只要保证对称，两者是完全等价的（这意味着在代码中也可以这样做）

.

但其实在OI竞赛中往往只用到它的一种变式：

可以发现，这其实就是对"系数"非线性的多项式系数序列做转换，这样我们往往就可以方便地对这些系数序列做操作。

1.5 注：

• 从这里开始不再以"j"表示\sqrt {-1}，而采用"i"表示

• 对于缺省的N,n总是指项数或数据规模

• 对于缺省的k总有0<k<N

• 下标总是从"0"开始

• 如果还有解释不到位的麻烦意会下QwQ

• 通常 DFT/IDFT 指离散傅里叶变换，而 FFT/IFFT 指快速（离散）傅里叶变换；前者强调转换，后者强调算法。不过其本质都是傅里叶变换

2. FFT

各位应该都熟知这样一个公式：

这即是欧拉公式；它将数学中一些常见的数联系在了一起。其实，傅里叶变换也正是借助它来联系三角函数的。

而我们还可以发现这个公式的一个特别之处：它可以将自然底数的指数形式化为三角函数的线性组合；这也正是傅里叶变换的原理。

.

接下来我们研究计算DFT的算法。

.

我们设\omega_n满足\omega^n=1，即\omega是n次单位根，同时n是最小的使w^n=1的数。

于是我们便有：

并设n表示多项式的项数

再将其代入DFT公式：

由于这两个式子是对称的，下文我们将重点研究上面的式子( DFT )，其结论对下面的式子( IDFT )也是适用的

.

我们需要先推出一些引理（引理引用的格式是在用时紧接着有"[数字]"）：

.

现在考虑按项数的奇偶性分组：

其中F_{x0},F_{x1}的规模为\frac n 2；这即是一个经典的递归模型。

.

但这样的递归式并不能减少运算，我们还要推出其他性质：

这即表明的变换序列中有某种映射关系，且这映射关系是一一对应的；同时我们还可以用O(1)计算它。

这样，我们就令运算规模就减小到了O(nlogn)。

.

我们整理下上面得到的式子，并设一些限定条件避免悖论：

其中n的规模还要求我们在算法实际应用时控制。

这递归式因其对称性也常被叫做蝴蝶变换。

.

这样，我们就可以根据这个递归式和变换序列间的映射关系设计一个O(nlogn)的算法计算 DFT/IDFT 了，

这即是 FFT/IFFT 。

3. FFT的实现

我们再讨论下如何实现这个代码，并做一些代码底层的小优化。

.

先观察上面的递归式，我们发现：

• 递归时，按奇偶分配值，下表示值的下标

• 回归时，下标代表变换后值的下标，直接顺序读取即可

这样我们可以初步写出代码：

//递归版 
//T指值类，cp是复数类，pi是圆周率 
void FFT(int n/*即当前规模*/, T *x){
    if(n == 1) return;
    T x0[n/2], x1[n/2];
    for(int k =0; k < n; k +=2) x0[k/2] =x[k], x1[k/2] =x[k+1];//奇偶分配值 
    FFT(n/2, x0), FFT(n/2, x1);

    cp wn(cos(2.0*pi/n), sin(2.0*pi/n))/*这里带入欧拉公式计算*/, w(1, 0);
    for(int k =0; k < n/2; ++k, w *=wn){
        //这里回归时已经是变换后的值了 
        x[k] =x0[k]+x1[k]*w;
        x[k+n/2] =x0[k]-x1[k]*w;
    }
}

但我们发现代码中多次转储数据，且这样的操作是随数据规模增长的，对算法效率很不利。

.

现在我们考虑递归时序列值序的变化规律。

先发现：

• 值唯一分配给递归底层的代码，返回时也不需要知道原来的值

于是考虑预处理数据序；

再打表观察：

• 设数据规模8，序号从0开始

• 递归前值的二进制序号：000 001 010 011 100 101 110 111

• 递归至底层后值的二进制序号：000 100 010 110 001 101 011 111

可以发现它是原来序号的二进制反序！

于是我们可以用dp预处理序号：

//l是最大序号的二进制位长度 
for(int i =0; i < n; ++i) s[i] =(s[i>>1]>>1)|((i&1)<<(l-1));

（仔细观察每个位运算符的作用应该很快就能理解）

这也叫做蝴蝶定理（但我真不会证QAQ，貌似是用归纳法）。

.

最后我们考虑下模拟回归过程时是否需要转储值。

观察上面的这段代码：

x[k] =x0[k]+x1[k]*w;
x[k+n/2] =x0[k]-x1[k]*w;

再结合我们之前做优化时已经按递归最底层的序列排序，这里我们只要在每次操作时转存下x[k]，x[k+n/2]的值就行了。

.

于是我们就可以写出非递归的最终代码：

//非递归版 
//T指值类，cp是复数类，pi是圆周率, op标记是否是逆变换 
inline void FFT(int N, T *x, int op){
    for(int i =0; i < N; ++i) if(i < s[i])//保证只交换一次 
        swap(x[i], x[s[i]]);

    //这里的lim是规模的1/2 
    for(int lim =1; lim < n; lim <<=1){
        cp wn(cos(pi/lim/*2约去了*/), op*sin(pi/lim));
        for(int s =0; s < N; s +=lim<<1){
            cp w(1, 0);
            for(int k =s; k-s < lim; ++k, w *=wn){
                cp res0 =x[k], res1 =x[k+lim]*w;//转存值 
                x[k] =res0+res1;
                x[k+lim] =res0-res1;
            }
        }
    }
}

4. 序列的操作

前文我们曾提到过 DFT/IDFT 是为了方便地对序列做操作，可这些操作究竟是那些？

首先， DFT 实际上就是把 信号 从 时域 转化为 频域 ，这在信号分析时可能有特殊含义。

不过上面的看不懂也没关系w。事实上，在竞赛中我们需要知道的常常是它的另外一点，用来做卷积重要的性质：对两个 DFT 后的序列做 逐项相乘 ，和对它们的原序列做 循环卷积 等价；

这里直接给出 循环卷积 的两种形式：

或C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}A[(n+k-i)\mod n]B[j](0\leq k < n)

或者可以形象地理解为将两个序列分别按 逆时针 和 顺时针 顺序放到两个同心圆上，逐位旋转，每种旋转状态下每对相对应的序列元素的积的和就是该项循环卷积的值；

还有一点，两个 DFT 后的序列做 逐项相加 也是和它们的原序列做 逐项相加 等价的。

.

而我们平常做的多项式积，其系数序列的变化就是在做线性卷积。我们直接按照多项式积理解就行了。

按做多项式积的结果最高次数而扩充系数序列（即补零），再对系数序列做循环卷积，结果和直接做线性卷积是等价的。（这段建议多读几次）

以上结论的证明可以到此文第三节查询，限于深度这里就不再多言

5. DFT与多项式

当然我们可以直接套卷积的结论。

不过为了给出一个更简单的推论方式，这里还是多写几笔：

.

我们设多项式A(x) = \sum\limits_{n=0}^{N-1}{a_nx^n},其系数序列为A[k]，

代入DFT，可得：

或者写成：

于是我们便可以用 FFT/IFFT O(nlogn)地对A[x]（即系数序列）做变换，之后就可以视转换后序列为给多项式带入不同的值（\omega）的结果的序列，O(n)地对系数序列操作。

但在实际操作时，要注意多项式的系数序列相乘后会增多，要拓充省略的0系数到系数序列。

对于序列规模n无法直接做FFT的解决方法也是直接扩充零系数。

就像这样计算：

while(N < la+lb) N <<=1;

（数组的默认值为零，就可以不用做“补零操作”了）

完整程序和上文其实也差不多，这里就不贴了（顺便实现下加强理解呗QwQ）。

6. FFT与NTT

上文中我们曾多次用到单位复根(w_n)：

观察这张图，可以发现 消去引理 的本质就是做复数指数时单位根绕原点不断地旋转，重复取值，从而可能进行简化运算。

.

(科普1：\phi(p)即是欧拉函数。)

(科普2：\delta_p(g)是指g模p的阶 (且g, p互质)；即对于 g^n \equiv 1 \pmod{p} 最小的正整数n。)

现在再引入 原根 的概念：设p是正整数，g是整数，且满足它们互质；

若\delta_p(g)=\phi(p)，则称g为模p意义下的一个原根。

同时对于质数还特别地有每个原根g的阶\delta_p(g)=p-1。（可由定义推得）

.

之所以用 原根 替代 单位根，是因为它的两点性质：

1. 根据上文 阶 的定义，原根也是“循环”的。

2. 原根的取值也会像单位复根一样不断重复。因为对于 g^s \equiv b \pmod{p}，且 0<s<\delta_p(g)，这些式子可以证明是 互不同余 的，即它们构成了模p的 简化剩余系 。

.

于是我们尝试找到一些便于二分的（质）数：998244353=2^{23}\cdot7\cdot17+1,1004535809=2^{21}\cdot479+1,469762049=2^{26}\cdot7+1。

这样我们就可以用g^{\frac {p-1} n}在模p意义下代替单位复根w_n做FFT/IFFT，即 NTT ( 快速数论变换 )。

但实际操作时，一定要注意：NTT 可以计算的答案是不会超过模数大小的。

.

另外998244353,1004535809,469762049这几个数都是质数，它们都有原根3；

这三个数是竞赛常用的，一般的NTT题目都会以它们做模数，是要背住的。

.

NTT的代码也很简单，只需改动下关于w_n的部分，所谓负指数其实就是分数。不过一定要注意转换前后的值不得大于模数。

7. *trick-打包FFT

仍是这张图，我们这次形象地理解一遍DFT/IDFT的过程：

我们可以以较生动的向量形式表示为

.

变换后

其实这个等式是转换前后始终满足的；因为DFT/IDFT就是对这些向量做旋转，而P,Q的变换方向我们保证是一致的。

.

现在我们要简化F_P,F_Q为一次DFT运算。

我们观察F_P[N-1]和F_Q[N-1]的值，发现它们正好是一对共轭复数；而我们又联想到F_P,F_Q的定义，便可以推测它们的值序列一定存在着某种共轭映射。

.

于是我们尝试做F_Q[k]的共轭复数

这不就是将P反向旋转么？！

即F_Q[k]=conj(F_P[N-1-k])（conj(z)为取z的共轭复数）（且对于k=N-1特别的有F_Q[N-1]=F_P[N-1]，就是我们上文提到的）。

.

整理一下

1.设P[k]=A[k]+iB[k],Q[k]=A[k]-iB[k]

2.有F_A[k]=\frac {F_P[k]+F_Q[k]} 2,F_B[k]=i\frac {F_Q[k]-F_P[k]} 2

3.有F_Q[k]=conj(F_P[N-1-k]) (0<k<N-1)（conj(z)为取z的共轭复数）, F_Q[N-1]=F_P[N-1]

这样我们便成功将DFT也两两打包为一次。

.

其实这种trick在 《再探快速傅里叶变换》--毛啸 中也有讲到；但原文是用公式推导的。

个人觉得还是图形的解释更利于理解和认识算法的本质。

后记

从初步接触到略有所成，我总共用了5天的时间（窝太菜了）才结束了FFT的学习。但越学习我才越发现自己的学识疏浅。

这是我写过（目前）最长的一篇博文（毕竟FFT实在是太多了），但毕竟我是个蒟蒻还有许多不足，文中仍多少有不严谨的地方（有些证明我是真不知道标准格式怎么写，各位大佬轻喷啊ORZ）。

另外一篇进阶技巧在这里（其实就是文章开头的QWQ）