OI笔记 | 快速傅里叶变换(FFT)
模板题：
本文大部分思路来源于 Reducible 的视频。
推荐观看：形象展示傅里叶变换 - 3Blue1Brown 的视频。
FFT 的概念和作用
快速傅立叶变换（Fast Fourier Transform，FFT）可以在 $O(n\log n)$ 的时间内解决多项式乘法问题，在 OI 中，是用来操纵生成函数，优化效率的重要算法。
当然，FFT 在其他领域也有广泛应用，甚至与核试验的检测有关，不过那是另外的事了。
如果要完全掌握 FFT 的思想与正确性证明，你可能需要熟练掌握线性代数、复数、积分等前置知识；又因为实际上傅里叶变换（Fourier Transform）是一种分析信号的方法，你可能需要了解时域与频域的转换等知识。
其实我也没有熟练掌握上述内容。然而如果只考虑它在 OI 中的应用，我们大可从多项式乘法出发来了解 FFT 的思想。
从多项式乘法谈起
我们最熟悉的表示多项式的方法是 系数表示：
\[A(x)=p_0 + p_1 x +p_2 x^2+ \dots + p_{n-1} x^{n-1} \to A(x)= \left[ p_0,p_1,\dots,p_{n-1} \right]\]用系数表示进行的多项式乘法为 $O(n^2)$ 的，需要把两个多项式的每个系数都相乘。
显然一个 $n$ 次多项式可以被 $t$ 个点确定，当且仅当 $t\ge n + 1$。这里不证明这个结论。我们把用点来表示多项式的方法称作 点值表示：
\[A(x)=[(x_0,A(x_0)),(x_1,A(x_1)),\dots(x_n,A(x_n))]\]用点值表示进行的多项式乘法为 $O(n)$ 的，只需要对相应的点值做乘法。
事实上，FFT 能做的就是把 系数表示 和 点值表示 的相互转换在 $O(n\log n)$ 内完成。
分割成偶函数进行分治
考虑什么情况下，我们已经知道一个点的点值 $(x_0,A(x_0))$，就能在 $O(1)$ 内得到另一个点值。可以想到如果这个多项式是偶函数，那么可以得出 $(-x_0,A(x_0))$ 是另一个点值。于是考虑把多项式分成两个有奇偶性的部分：
\[A(x)=(p_0 + p_2x^2 + p_4x^4 + \dots + p_{n-2}x^{n-2}) +\\ x(p_1 + p_3x ^2 + p_5x^4+\dots +p_{n-1}x^{n-2})\]这里假定 $n$ 是 $2$ 的幂，这么取值方便进行分治。
不妨设 $A_e(x)=(p_0+p_2x+p_4x^2+\dots + p_{n-2}x^{n/2-1})$，$A_o(x)=(p_1 + p_3x+p_5x^2+\dots + p_{n-1}x^{n/2-1})$，则有：
\[A(x)=A_e(x^2)+xA_o(x^2)\]并且这里的 $A_e(x^2), A_o(x^2)$ 都是偶函数，这个性质很妙。
而且上面这个式子说明我们如果可以递归地用 FFT 处理完 $A_e$ 和 $A_o$ 的点值表示，就可以在 $O(n)$ 的时间内合并答案，这样不就把时间复杂度降到了 $O(n\log n)$ 了吗。
用单位根取点
但是还有一个问题，就是你会发现我们括号里的是 $x^2$，如果 $x$ 是实数，则有 $x^2\ge 0$，那么我们递归完一层就递归不下去了，因为我们的设想是通过偶函数的性质取 $x_0$ 和 $-x_0$ 这一个点对来分治，但是 $-x_0\le 0$，不能成为另一个数的平方。
所以我们考虑使用复数，而且取点方式很巧妙，我们直接取 $x^n=1$ 的 $n$ 个单位根作为我们的点值。为什么这样可行呢？
取 $\omega=e^{\frac{2\pi i}{n}}=\cos \frac{2\pi}{n}+i \sin{\frac{2\pi}{n}}$。据复平面的知识，$n$ 次单位根把单位圆 $n$ 等分，容易知道 $\omega ^k(k\in {0,1,\dots, n-1})$ 就是逆时针旋转得到的第 $k$ 个根。
然后单位根有 $3$ 个巧妙的性质，根据复平面的几何直观很好理解：
\[\begin{aligned} \omega_n^n&=1\\ \omega_n^k&=\omega_{2n}^{2k}\\ \omega_{2n}^{k+n}&=-\omega_{2n}^k\\ \end{aligned}\]那么当你递归地处理完 $P_e = [A_e(\omega^0),A_e(\omega^2),A_e(\omega^4),\dots,A_e(\omega^{n-2})]$，$P_o = [A_o(\omega^0),A_o(\omega^2),A_o(\omega^4),\dots,A_o(\omega^{n-2})]$ 后，可以顺利的用上面的公式合并出 $A(x)=P_e(x)+\omega^x P_o(x)$，$A(x + n/2)=P_e(x)-\omega^x P_o(x)$（这是因为单位根的性质， $\omega^{x+n/2}=-\omega^{x}$）。
算法的时间复杂度 $O(n\log n)$。
快速傅里叶逆变换
傅里叶逆变换可以用傅里叶变换表示。对此本文不想提供理解，因为我也不是很理解。
有没有大佬来告诉我一种简单的理解方式！
总之，代码实现上，你把 FFT 中的 $\omega$ 从 $e^{\frac{2\pi i}{n}}$ 改成 $e^{\frac{-2\pi i}{n}}$ ，就变成了 IFFT。很神奇吧。
代码实现
先手写一个复数。
struct cp {
double x, y;
cp(double _x = 0, double _y = 0) : x(_x), y(_y) {}
cp operator + (const cp &t) const {return cp(x + t.x, y + t.y);}
cp operator - (const cp &t) const {return cp(x - t.x, y - t.y);}
cp operator * (const cp &t) const {return cp(x * t.x - y * t.y, x * t.y + y * t.x);}
};
按照上面的思路，容易写出递归版本：
void FFT(cp *P, int x, int type) {
if(x == 1) return;
cp Pe[x], Po[x];
for(int i = 0; i <= x; i += 2)
Pe[i >> 1] = P[i], Po[i >> 1] = P[i + 1];
FFT(Pe, x >> 1, type);
FFT(Po, x >> 1, type);
cp omega = cp(cos(PI * 2.0 / x), type * sin(PI * 2.0 / x)), now = cp(1, 0);
for(int i = 0; i < (x >> 1); i++) {
cp u = Pe[i], v = now * Po[i];
P[i] = u + v;
P[i + (x >> 1)] = u - v;
now = now * omega;
}
}
当然也有基于位逆序置换(蝴蝶变换)的迭代版本，不太好理解，不过效率会高一点。这里直接给出代码。
void pre(){
while(len <= n + m)
len <<= 1, k++;
for(int i = 0; i <= len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
void FFT(cp *P, int type) {
pre();
for(int i = 0; i < len; i++)
if(i < rev[i]) swap(P[i], P[rev[i]]);
for(int i = 1; i < len; i <<= 1) {
cp omega(cos(PI / i), type * sin(PI / i));
for(int k = 0; k < len; k += (i << 1)) {
cp now(1, 0);
for(int l = 0; l < i; l++) {
cp u = P[k + l], v = now * P[k + i + l];
P[k + l] = u + v;
P[k + i + l] = u - v;
now = now * omega;
}
}
}
}
♥