FFT学习笔记
0、前言
首先，我们要知道\(FFT\)是用来解决多项式乘法问题的
比如对于两个多项式\(f(x)\),\(g(x)\)
\(f(x)=\sum_{i=0}^{n-1}a_ix^i\)
\(g(x)=\sum_{i=0}^{m-1}b_ix^i\)
求\(f\times g\)
1、系数表示法
一个\(n-1\)次\(n\)项多项式\(f(x)\)可以表示为\(f(x)=\sum_{i=1}^{n-1}a_ix^i\)
所以我们可以用每一项的系数来表示\(f(x)\)，即
2、点值表示法
我们对于每一个\(x\)都可以求出一个对应\(x\)的\(f(x)\)值，我们把这样的\(x\)与\(f(x)\)当成\((x,y)\)，就会有很多个点
显然，这\(n\)个点是可以确定一个这个多项式的(\(n\)元\(n\)次方程的解)
因此，这个多项式还可以用\(n\)个点来表示
\(f(x)= \left\{ (x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1})) \right\}\)
点值表示法相比系数表示法有什么优势呢？
我们考虑对于系数表示法我们怎么求两个多项式乘积
必然是枚举\(f\)和\(g\)的每一项分别相乘，时间复杂度\(O(n^2)\)
但对于点值表示法来说
我们记\(h(x)\)为\(f\)和\(g\)的乘积，则有
显然，这样算我们就只需要\(O(n)\)的时间复杂度了
这时候，问题就出在了如何将系数表示法，转化成点值表示法
前置知识
- 复数相乘，模长相乘，极角相加
2、单位圆上的点满足\(z=cos\theta + isin\theta\)(\(\theta\)为直线与\(x\)轴的交点)
Tip:从这里开始之后所有n都默认是2的整数次幂
DFT（离散傅里叶变换）
对于任意系数多项式转点值，当然可以随便取任意\(n\)个\(x\)值代入计算
显然暴力计算的复杂度是\(O(n^2)\)的
那么应该选取那些\(x\)的值可以降低时间复杂度呢
考虑一下，入果我们带入的\(x\)可以是每个\(x\)的若干次方等于1，我们就可以减少很多次运算，那么那些数的若干次方会是1呢
傅里叶给出了答案，他说在单位圆上的点都是满足的
我们把这个圆\(n\)等分，就得到了\(n\)个满足条件的点
将这\(n\)个点从\((1,0)\)开始编号
记\(\omega_n^1\)为\(n\)次单位根，编号为\(k\)的点的复数值为\(\omega_n^k\)
由复数乘法模长相乘，极角相加可得\(\omega_n^k=(\omega_n^1)^k\)
且有\(\omega_n^k=cos\frac{k}{n}2\pi+isin\frac{k}{n}2\pi\)
那么\(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\)即为我们要带入的\(x\)
性质
- \(\omega_n^{x+y}=\omega_n^x\times \omega_n^y\)
证：
\(
\begin{aligned}
&\omega_n^x\times \omega_n^y=(cos\frac{x}{n}2\pi+isin\frac{x}{n}2\pi)\times (cos\frac{y}{n}2\pi+isin\frac{y}{n}2\pi) \\
&~~~~~~~~~~~~~~=cos\frac{x}{n}2\pi cos\frac{y}{n}2\pi-sin\frac{x}{n}2\pi sin\frac{y}{n}2\pi+icos\frac{x}{n}2\pi sin\frac{y}{n}2\pi+isin\frac{x}{n}2\pi cos\frac{y}{n}2\pi\\
&~~~~~~~~~~~~~~=cos\frac{x+y}{n}2\pi+isin\frac{x+y}{n}2\pi\\
&~~~~~~~~~~~~~~=\omega_n^{x+y}
\end{aligned}
\)
- \(\omega_n^k=\omega_{2n}^{2k}\)
证：
\(\omega_n^k=cos\frac{k}{n}2\pi+isin\frac{k}{n}2\pi=cos\frac{2k}{2n}2\pi+isin\frac{2k}{2n}2\pi=\omega_{2n}^{2k}\)
- \(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)
证：
\(
\begin{aligned}
&\omega_n^{\frac{n}{2}}=cos\frac{\frac{n}{2}}{n}2\pi+isin\frac{\frac{n}{2}}{n}2\pi\\
&~~~~~~=cos\pi+isin\pi\\
&~~~~~~=-1
\end{aligned}
\)
\(
\begin{aligned}
&\omega_n^{k+\frac{n}{2}}=\omega_n^k \times \omega_n^{\frac{n}{2}}\\
&~~~~~~~~~=-\omega_n^k
\end{aligned}
\)
- \(\omega_n^0=\omega_n^n=1\)
略
- \((\omega_n^x)^y=(\omega_n^y)^x\)
证：
左右式结果都等于\(\omega_n^{xy}\)
FFT（快速傅里叶变换）
尽管\(DFT\)里\(x\)的取值看起来很牛逼，但实际上还是暴力\(O(n^2)\)
我们考虑进一步优化
我们把多项式写出来
把它按照奇偶性分组
发现两边长得很像，于是再设两个多项式
显然
不妨设\(k<\frac{n}{2}\)，将\(x=\omega_n^k\)带入\(A(x)\)
再将\(x=\omega_n^{k+\frac{n}{2}}\)带入\(A(x)\)
我们发现\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)只有后面那个东西不一样
也就是说，只要我们知道\(A_1(\omega_{\frac{n}{2}}^k)\)和\(A_2(\omega_{\frac{n}{2}}^k)\)，就可以计算\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)的值
显然，这东西是可以递归分治来把系数表示法转换成点值表示法
时间复杂度\(O(nlogn)\)
IFFT（快速傅里叶逆变换）
当我们用点值表示法求出两个多项式的乘积后，我们发现一个问题
这玩意儿怎么转回系数表示法呢？
显然，傅里叶选取\(x=\omega_n^k\)带入多项式肯定不仅仅是因为那些性质
考虑原来的多项式是\(A(x)=\sum_{i=0}^{n-1}a_ix^i\)
记\(b_i=A(\omega_n^i),i\in\left\{0,1,...,n=1\right\}\)
构造一个新的多项式
相当于把\(\left\{b_0,b_1,...,b_{n-1}\right\}\)当作\(f\)的系数表示法
设\(c_i=\omega_n^{-i}\)，则\(A\)的点值表示法可以写成\(\left\{(c_0,A(c_0)),(c_1,A(c_1)),...,(c_{n-1},A(c_{n-1}))\right\}\),则有
记\(S(\omega_n^k)=\sum_{i=0}^{n-1}(\omega_n^k)^i\)
当\(k=0\)时，显然\(S(\omega_n^k)=n\)
当\(k\not =0\)时，有
\((2)-(1)\)得
即
所以有
也就是说对于给定点\(c_i=\omega_n^{-i}\),\(f\)的点值表示法为
因此，我们只需要对于\(\left\{b_0,b_1,b_2,...,b_{n-1}\right\}\)跑一遍\(FFT\)，然后每一项都除以\(n\)即可还原系数表示法
迭代版FFT
朴素的\(FFT\)还可以从分治的角度再优化一步
由于朴素的\(FFT\)要递归地将整个多项式的奇数次项和偶数次项分开知道只剩一个系数，这样既浪费空间又浪费时间
我们可以提前处理出每个系数到最后会拆分到哪个位置，然后把递归的过程变成合并的过程
以8项多项式为例，模拟拆分的过程：
- 初始序列为\(\left\{0(000),1(001),2(010),3(011),4(100),5(101),6(110),7(111)\right\}\)
- 一次递归后\(\left\{0(000),2(010),4(100),6(110)\right\},\left\{1(001),3(011),5(101),7(111)\right\}\)
- 两次递归后\(\left\{0(000),4(100)\right\},\left\{2(010),6(110)\right\},\left\{1(001),5(101)\right\},\left\{3(011),7(111)\right\}\)
- 三次递归后为\(\left\{0(000)\right\},\left\{4(100)\right\},\left\{2(010)\right\},\left\{6(110)\right\},\left\{1(001)\right\},\left\{5(101)\right\},\left\{3(011)\right\},\left\{7(111)\right\}\)
我们惊喜的发现，每个数最后的位置其实就是把二进制对称翻转一下，即\(001\rightarrow 100\)
这样就是FFT的究极体了
void FFT(complex<double> *a,int n,int inv) {
int bit=0;
while((1<<bit)<n) bit++;
for (int i=0;i<n;i++) {
pos[i]=(pos[i>>])>>1|((i&1)<<(bit-1));
if (i<pos[i]) swap(a[i],a[rev[i]]);
}
for (int Mid=1;Mid<n;Mid*=2) {
complex<double> tmp(cos(pi/Mid),inv*sin(pi/Mid));
for (int i=0;i<n;i+=Mid/2) {
complex(double) omega(1,0);
for (int j=0;j<Mid;j++,omega*=tmp) {
complex<double> x=a[i+j],y=omega*a[i+j+Mid];
a[i+j]=x+y;a[i+j+Mid]=x-y;
}
}
}
}