从零开始的FFT（ACM/OI）
1. 什么是DFT，FFT
离散傅里叶变换（Discrete Fourier Transform，DFT）傅里叶分析方法是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，通过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。
——百度百科
关于DFT在信号分析上的应用，本文并不涉及（因为我也不会）。
在我个人看来，DFT是一个将多项式从系数表示法转换为点值表示法的算法。而FFT（Fast Fourier Transform）顾名思义，即快速计算DFT的算法。
2. DFT，FFT有什么用
在算法竞赛中，我们常需要解决多项式A乘以多项式B并得到新多项式C的问题。
使用常规的计算方法，将多项式A的每一项乘以多项式B的每一项并相加，所得的新多项式即为多项式C。假设多项式A和B均有n项，则两个多项式相乘的时间复杂度为\(O(n^2)\)
常规方法计算两多项式乘积的速度令人并不满意，有没有更快的计算方法呢，答案是：有。
而这就需要用到接下来要讲的DFT和FFT。
3. 一些前置知识
学习DFT和FFT之前，需先了解以下概念：
1. 点值表示法
2. 复数
3.1 点值表示法
通常，我们表示一个有n项的多项式常使用以下方法：
这种表示方法我们称为系数表示法。
除此之外，还有另一种表达多项式的方法，即上文提到的点值表示法。
由数学知识可知：两点确定一条直线。即2个点可以确定一个1次多项式，相应的，n+1个点可确定一个n次多项式(带上常数项共有n+1个系数)。那么我们可以通过一系列的点表示一个多项式。
表示方式如下：
举个例子：
这里有一个二次多项式，其系数表示法为：
我们任意代入三个数：1、2、3。则
换个写法，(1,0) , (2,7) , (3,18) (没错，就是三个点的坐标)
所以，点值表示法为
3.2 复数
3.2.1 复数的定义
复数分为两部分，分别为实数和虚数，实数不多做赘述，首先定义虚数单位 \(i\)，对于虚数单位 \(i\)，有\(i^2 = -1\) 。
则复数表示为 \(a+bi\) ，其中\(a\)为实数，\(bi\)为虚数。
将\(a\)视为坐标轴上的\(x\)坐标，\(b\)视为\(y\)坐标，复数可表示为平面上的点。
（偷个图）
图中横轴代表实数，纵轴代表虚数，则图中所表示复数为\(2+3i\)
我们称图中红线所代表长度为复数的模，红线与x轴的夹角称为复数的辐角。
3.2.2 复数的运算
加法 ：\(\large (a+bi)+(c+di)=(a+c)+(b+d)i\)
减法 ：\(\large (a+bi)-(c+di)=(a-c)+(b-d)i\)
乘法 ：\(\large (a+bi)\times(c+di)=(ac-bd)+(ad+bc)i\)
除法 ：\(\large \dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)+(bc-ad)i}{c^2+d^2}\)
重点来了
重点来了
重点来了
再看看这张图，不难发现，复数也可由红线长度\(L\)和角度\(\theta\)表示（类似于极坐标）
对于复数的乘法运算，两复数的模相乘，辐角相加 （重点）
使用三角函数，一个复数也可表示为： \(\large Lcos\theta+L(isin\theta)\) （重点）
4. DFT
终于到正片了。。。
前文提到，使用常规的算法计算两多项式相乘的时间复杂度是\(O(n^2)\)，但将多项式转换为点值表示法再计算，如：
和
相乘，则得到的多项式为：
不难发现，使用点值表示法进行计算，时间复杂度为\(O(n)\)。
虽然两多项式相乘的时间复杂度降为\(O(n)\)，但为了将系数表示法转换为点值表示法，一个n次多项式需要代入n+1个数进行计算，其时间复杂度仍为\(O(n^2)\)。
有没有更好的办法呢。
从此处开始，默认n为2的整数次幂（实际运算中若n不是2的整数次幂就补系数为0的项，直到n是2的整数次幂）
从此处开始，默认n为2的整数次幂（实际运算中若n不是2的整数次幂就补系数为0的项，直到n是2的整数次幂）
从此处开始，默认n为2的整数次幂（实际运算中若n不是2的整数次幂就补系数为0的项，直到n是2的整数次幂）
其实可以代入一组神奇的x，代入以后不用做那么多的次方运算
这些x当然不是乱取的，而且取这些x值应该就是 傅里叶 的主意了。考虑一下，如果我们代入一些x，使每个x的若干次方等于1，我们就不用做全部的次方运算了。
代入1和-1可以达到我们的目的，若考虑复数，\(i\) 和\(-i\) 同样可以，但是只有四个值远远不够。
傅里叶提出：在这个圆上都可以
以原点为圆心，画一半径为1的圆，同样，x轴代表实数，y轴代表虚数。
取值在该圆上的复数，其若干次方必可取值为1。不难理解，因取值在圆上，则模为1，其若干次方的模也必定为1，又由 \(k\theta\equiv0(mod\ 360)\) 可得，必定存在一整数\(k\)使得该等式成立，即经过\(k\)次旋转后辐角为0。
接下来的问题，就是在这个圆上都要取哪些点，答案很简单，只需将这个圆等分成n份，每份之间的分界点既是我们要取的点。
举个例子：假设\(n=8\)，则需将圆等分成8份，如图：
同时，我们将图中8个点所代表的复数按照辐角大小的顺序命名为:
我们称第k个n次单位根为 \(\omega_n^k\) ，其中 \(\omega_n^0\)为1。
若仅是将取到的这些特殊值代入，的确可以减少一定的计算量，但在渐进意义上直接代入求值的时间复杂度仍为\(O(n^2)\) ，这时，就该FFT出场了，它可以将两多项式求乘积的时间复杂度降为\(O(nlogn)\) 。
5. FFT
首先是给出两个单位根的性质：
\(\large\omega_{2n}^{2k}=\omega_n^k\)
证明： \(\large \omega_n^k=cos(\frac{k}{n}\cdot2\pi)+isin(\frac{k}{n}\cdot2\pi)=cos(\frac{2k}{2n}\cdot2\pi)+isin(\frac{2k}{2n}\cdot2\pi)=\omega_{2n}^{2k}\)
\(\large \omega_n^{k+\frac{n}{2}}=-\omega_n^k\)
证明： \(\large \omega_n^{k+\frac{n}{2}}=\omega_n^k\cdot\omega_n^{\frac{n}{2}}=\omega_n^k\cdot-1=-\omega_n^k\)
对于多项式
将里面的项换一下位置，将\(x\)的指数为偶数的项放在左边，\(x\)的指数为奇数的项放右边，可以得到
设两多项式\(A_0(x)\)和\(A_1(x)\)，使得
也就是\(A_0(x)\) 为\(A(x)\) 中\(x\)的指数为奇数次的，\(A_1(x)\) 为\(A(x)\) 中\(x\)的指数为偶数次的。
不难发现，\(\large A(x)=A_0(x^2)+xA_1(x^2)\)
也就是说\(A(x)\)变为两个项数为原来一半的多项式的和。
将我们之前取到的特殊值代入：
又因为\(\large \omega_{2n}^{2k}=\omega_n^k\) ，
没错，这意味着我们可以用递归计算一个多项式的值了，递归到\(n=1\)时直接返回就行。
而且又因为\(\large \omega_n^{k+\frac{n}{2}}=-\omega_n^k\) 这个性质，我们计算一个n次多项式时仅递归求出前\(\large \frac{n}{2}\) 项即可，后\(\large \frac{n}{2}\) 项可直接由前\(\large \frac{n}{2}\) 求出。
这样，我们计算代入一个数的时间复杂度为\(O(logn)\)，我们需要代入\(n\)个数计算，则最终时间复杂度为\(O(nlogn)\)。
附上代码：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
double pi=acos(-1);
complex<double> a[5000005],b[5000005];
void fft(complex <double> *arr,int len)//len是数组arr的长度
{
if(len==1) return;//数组长度为1时直接返回
complex<double> a0[(len>>1)+5],a1[(len>>1)+5];
for(int i=1;i<=(len>>1);i++)//将arr数组中的偶次项系数和奇次项系数分别赋给a0数组和a1数组
{
a0[i]=arr[(i<<1)-1];
a1[i]=arr[i<<1];
}
fft(a0,len>>1);//递归处理a0和a1
fft(a1,len>>1);
complex<double> w0(1,0);//初始值
complex<double> w(cos(2*pi/len),sin(2*pi/len));//w0每次需旋转的角度
for(int i=1;i<=(len>>1);i++,w0*=w)//循环每次都使w0旋转一次
{
arr[i]=a0[i]+w0*a1[i];
arr[i+(len>>1)]=a0[i]-w0*a1[i];//使用单位根的性质直接求出数组后半部分
}
return;
}
6. IDFT（傅里叶逆变换）
我们已经将多项式转换为点值表示，并经过了相乘得到新多项式，但它现在仍是点值表示法，我们要将新多项式从点指表示转换回系数表示。
对于前面的DFT的过程，可以用矩阵表示为
最右边和最左边的矩阵已知，我们需要求中间那个矩阵。也就是先求最左边那个矩阵的逆矩阵\(R\)，再将\(R\)乘右边的矩阵可得到中间的矩阵。
对于左边矩阵的逆矩阵\(R\)，有一个特殊的性质，\(R\) 中的每一项恰等于左边矩阵每一项的倒数再除以\(n\) ，即：
以上这些操作，在代码中实现十分简单，只需将前面给出的FFT代码稍作修改即可。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
double pi=acos(-1);
complex<double> a[5000005],b[5000005];
void fft(complex<double> *arr,int len,int o)//传入的参数多了一个o，o的值决定了是进行DFT操作还是IDFT操作
{
if(len==1) return;
complex<double> a0[(len>>1)+5],a1[(len>>1)+5];
for(int i=1;i<=(len>>1);i++)
{
a0[i]=arr[(i<<1)-1];
a1[i]=arr[i<<1];
}
fft(a0,len>>1,o);
fft(a1,len>>1,o);
complex<double> w0(1,0);
complex<double> w(cos(2*pi/len),o*sin(2*pi/len));//当o=1时代表是系数转点值，o=-1时代表点值转系数
for(int i=1;i<=(len>>1);i++,w0*=w)
{
arr[i]=a0[i]+w0*a1[i];
arr[i+(len>>1)]=a0[i]-w0*a1[i];
}
return;
}
计算两个多项式的乘积主函数这么写就行：
（摘自 洛谷P3803 的代码）
int main()
{
ios::sync_with_stdio(0);
cin.tie(0);
int n,m;//两个多项式的次数
while(cin>>n>>m)
{
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
for(int i=1;i<=n+1;i++)
cin>>a[i];
for(int i=1;i<=m+1;i++)
cin>>b[i];
int len=1;
while(len<=n+m) len<<=1;//长度补足至2的整数次幂
fft(a,len,1);//DFT操作o=1
fft(b,len,1);
for(int i=1;i<=len;i++)
a[i]*=b[i];
fft(a,len,-1);//IDFT操作o=-1
for(int i=1;i<=(n+m+1);i++)
cout<<(int)(a[i].real()/len+0.5)<<" ";//防止浮点数误差加一个0.5
cout<<"\n";
}
return 0;
}
参考资料：
FFT的迭代写法以后再更
浙公网安备 33010602011771号