FFT 初步
Hexarhy
2021-03-29 20:25:39
# Introduction
快速傅里叶变换 (Fast Fourier Transform，FFT)，一种支持在 $O(n\log n)$ 时间内计算两个 $n$ 次多项式乘法的算法。
**前置知识：** 复数与复平面，矩阵初步。
# Analysis
我们需要借助一些工具来帮助我们计算多项式。
**符号约定：**
- $a/b$，表示 $a$ 除以 $b$ 向下取整。
- 为了区分，斜体的 $i$ 是一个数学变量，而正体的 $\rm i$ 表示虚数单位。
## 多项式的表示
### 系数表示法
即用系数表示一个多项式。
$$
F(x)=\sum_{i=0}^n a_ix^i \Leftrightarrow F(x)=\{a_0,a_1,\cdots,a_n\}
$$
显然，一个确定的多项式与系数构成双射。
### 点值表示法
在函数 $F(x)$ 上取 $n+1$ 个横坐标不同的点来表示，也就是：
$$
F(x)=\sum_{i=0}^n a_ix^i \Leftrightarrow F(x)=\{(x_0,F(x_0)),(x_1,F(x_1)),\cdots,(x_n,F(x_n))\}
$$
一个确定的多项式也与这 $n+1$ 个点构成双射。
证明：考虑待定系数法。用这 $n+1$ 个点可以组成一个 $n+1$ 元一次方程组。由于所有点的横坐标都互不相同，不存在等价或矛盾的方程，因此能解出唯一解，也就对应了唯一一个多项式。
### 基本原理
有了先进便捷的工具，我们可以得到 FFT 的基本原理：
- 系数表示 $\xrightarrow{\rm DFT} $ 点值表示
- 在点值上运算
- 点值表示 $\xrightarrow{\rm IDFT}$ 系数表示
当然，为了加快速度，我们要取特殊的横坐标，并用**分治**计算。
## 复平面与单位根
### 单位根的表示
对于同样的横坐标选取，用**点值表示**两个多项式 $F(x),G(x)$ 。那么 $H(x)=F(x)\cdot G(x)$ 也容易用点值表示：
$$
H(x)=\{(x_0,F(x_0)G(x_0)),(x_1,F(x_1)G(x_1)),\cdots,(x_n,F(x_n)G(x_n))\}
$$
横坐标要怎么取才能快？$1,-1$ 的幂都很好算。但限制在实数里实在找不到 $n+1$ 个数这么多。联想到方程 $x^n=1$，我们取复数进行运算。
复平面上复数乘法的几何意义告诉我们，$1$ 的 $n$ 次方根即把单位圆等分成 $n$ 份，也就是说，对于第 $k(k=0,1,\cdots,n-1)$ 个根 $\omega_k$，根据复数的三角表示和欧拉公式，有：
$$
\omega_k=\cos\dfrac{2\pi}{k}+\mathrm{i}\sin\dfrac{2\pi}{k}=e^\frac{2\pi\mathrm{i}}{k}
$$
### 单位根的性质
单位根有非常多性质，对于 $\forall n\in \mathbb{N}^+,k\in \mathbb Z$，
- $\omega_{n}^n=1$
证明：用三角表示法和复数乘法的几何意义即得证。
$$
\omega_{n}^n=\cos^n\dfrac{2\pi}{n}+\mathrm{i}\sin^n\dfrac{2\pi}{n}=\cos{2\pi}+\mathrm{i}\sin2\pi=1
$$
- $\omega_{2n}^{2k}=\omega_n^k$
证明：同理。
$$
\begin{aligned}
\omega_{2n}^{2k}
&=\cos^{2k}\dfrac{2\pi}{2n}+\mathrm{i}\sin^{2k}\dfrac{2\pi}{2n}\\
&=\cos^k \dfrac{2\cdot 2\pi}{2n}+\mathrm{i}\sin^k\dfrac{2\cdot 2\pi}{2n}\\
&=\cos^k\dfrac{2\pi}{n}+\mathrm{i}\sin^k\dfrac{2\pi}{n}\\
&=\omega_n^k
\end{aligned}
$$
- $\omega_{2n}^{k+n}=-\omega_{2n}^k$
证明：同理，中间用到了诱导公式。
$$
\begin{aligned}
\omega_{2n}^{k+n}&=\cos^{k+n} \dfrac{2\pi}{2n}+\mathrm{i}\sin^{n+k}\dfrac{2\pi}{2n}\\
&=\cos \dfrac{(k+n)2\pi}{2n}+\mathrm{i}\sin\dfrac{(k+n)2\pi}{2n}\\
&=\cos\left(\dfrac{2k\pi}{2n}+\pi\right)+\mathrm{i}\sin\left(\dfrac{2k\pi}{2n}+\pi\right)\\
&=-\cos\dfrac{2k\pi}{2n}-\mathrm{i}\sin\dfrac{2k\pi}{2n}\\
&=-\omega_{2n}^k
\end{aligned}
$$
- $\omega_{n}^{k+n}=\omega_n^k$
证明：$\omega_{n}^{k+n}=\omega_{n}^{k+n/2+n/2}=-\omega_n^{k+n/2}=\omega_n^k$。当然用三角表示也容易证明。
## DFT
我们取横坐标 $x=\omega_n^k$ 来转化成点值表示。分治思想体现在将多项式按次数奇偶处理。
举例来说，对于 $F(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5$，转化为：
$$F(x)=(a_0+a_2x^2+a_4x^4)+(a_1x+a_3x^3+a_5x^5)$$
$$=(a_0+a_2x^2+a_4x^4)+x(a_1+a_3x^2+a_5x^4)$$
$$=G(x^2)+xH(x^2)$$
$$G(x)=a_0+a_2x+a_4x^2$$
$$H(x)=a_{1}+a_3x+a_5x^2$$
利用单位根的性质得到：
$$
\begin{aligned}
F(\omega_{n}^{k})&=G(\omega_n^{2k})+\omega_{n}^{k}H(\omega_n^{2k})\\
&=G(\omega_{n/2}^k)+\omega_{n}^kH(\omega_{n/2}^k)\\
F(\omega_{n}^{k+n/2})&=G(\omega_n^{2k+n})+\omega_{n}^{k+n/2}H(\omega_n^{2k+n})\\
&=G(\omega_{n}^{2k})-\omega_{n}^kH(\omega_{n}^{2k})\\
&=G(\omega_{n/2}^k)-\omega_{n}^kH(\omega_{n/2}^{k})
\end{aligned}
$$
递归处理即可。每次递归时间复杂度为 $O(\log n)$。
但是递归处理要求多项式次数必须为 $2^m-1$。所以就在多项式后面用系数 $0$ 补齐即可。
代入 $2^m$ 个单位根 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}(n=2^m)$ 即可得到多项式的点值表示。
显然，DFT 的时间复杂度为 $O(n\log n)$。
### 位逆序置换
依然按照奇偶次数分治。我们并不需要像 DFT 那样实际操作，只需要知道递归后的系数对应的位置即可。
以 $7$ 次多项式为例：
- 初始：$\{x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8\}$
- 第一次：$\{x_0,x_2,x_2,x_4\},\{x_1,x_3,x_5,x_7\}$
- 第二次：$\{x_0,x_4\},\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\}$
- 第三次：$\{x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\}$
找规律过程略去，毕竟我们只关心结论。
每个数初始下标用二进制表示，翻转 (reverse)，得到的数字即为最终的下标。这个称为**位逆序置换**。
直接根据这个规律模拟，时间复杂度 $O(n\log n)$。
事实上，还有一种更优的实现方式，采用了递推。设 $R(x)$ 为二进制下 $x$ 翻转后的数。记 $k$ 为二进制下 $x$ 的长度。
已知 $R(0)=0$。从小到大递推，已知 $R(x/2)$ 的值，将 $x$ 右移一位，取反，右移一位，就得到了二进制下 $x$ **个位之外**的翻转结果。例如：$x=4=(100)_2$ 时，翻转后的结果是 $(001)_2=1$。根据刚才的操作演示一遍：
$$
(100)_2\xrightarrow{\rm Move\ Right}(010)_2\xrightarrow{\rm Reverse}(010)_2\xrightarrow{\rm Move\ Right}(001)_2\xrightarrow{\rm Calculate}(001)_2
$$
对于个位，如果是 $0$ ，翻转之后最高位为 $0$；如果为 $1$，则最高位为 $1$，等价于加上了 $2^{k-1}$。综上，得到递推式：
$$
R(x)=\begin{cases}R(x/2)/2&,2\mid x\\
R(x/2)/2+2^{k-1}&,2\nmid x
\end{cases}
$$
时间复杂度 $O(n)$。
## IDFT
将点值转化为系数，有多种理解方式。个人觉得最好理解的是用矩阵运算。
先下结论：**DFT 中 $\omega_{n}^1$ 换成 $\omega_n^{-1}$，再除以 $n$**。
对于所取的 $n$ 个单位根，可以得到 $n$ 个方程，其中第 $k(k=0,1,2,\cdots,n-1)$ 条为：
$$
F(\omega_n^k)=\sum_{i=0}^{n-1} a_i(\omega_{n}^{k})^i
$$
咦，这个形式好熟悉啊。就是矩阵运算！实际上 DFT 本身是线性变换，所以就有了以下理解方式：
$$
\left[\begin{array}{c}
F(\omega_n^0) \\
F(\omega_n^1) \\
F(\omega_n^2) \\
F(\omega_n^3) \\
\vdots \\
F(\omega_n^{n-1})
\end{array}\right]=\left[\begin{array}{cccccc}
(\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \omega_n^0 & \cdots & (\omega_n^0)^{n-1} \\
(\omega_n^1)^0 & (\omega_{n}^{1})^1 & (\omega_{n}^1)^{2} & (\omega_{n}^1)^{3} & \cdots & (\omega_{n}^{1})^{n-1} \\
(\omega_n^2)^0 & (\omega_{n}^{2})^1 & (\omega_{n}^{2})^2 & (\omega_{n}^{2})^3 & \cdots & (\omega_{n}^{2})^{n-1} \\
(\omega_n^3)^0 & (\omega_{n}^{3})^1 & (\omega_{n}^{3})^2 & (\omega_{n}^{3})^3 & \cdots & (\omega_{n}^3)^{n-1} \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
(\omega_n^{n-1})^0 & (\omega_{n}^{n-1})^1 & (\omega_{n}^{n-1})^2 & (\omega_{n}^{n-1})^3 & \cdots & (\omega_{n}^{n-1})^{n-1}
\end{array}\right]\left[\begin{array}{c}
a_{0} \\
a_{1} \\
a_{2} \\
a_{3} \\
\vdots \\
a_{n-1}
\end{array}\right]
$$
现在已经知道了左边的矩阵 $A$ 和中间的左乘的矩阵 $B$，要求系数也就不难。也就是将 $A$ 乘上 $B$ 的逆矩阵。
至于逆矩阵怎么求，一种可行的方法采用了拉格朗日插值，参考[这里](https://www.zhihu.com/question/309243881/answer/575306025)。
我们要求的逆矩阵就是 $B$ 中的每一项取倒数再乘上 $\dfrac{1}{n}$。至于为什么是正确的，代入计算验证即可。
取倒数相当于指数乘上 $-1$，因此我们可以在 DFT 中传参数，$1$ 为 DFT，$-1$ 为 IDFT 来简洁地实现。
# Code
题目链接：[P3803 【模板】多项式乘法（FFT）](https://www.luogu.com.cn/problem/P3803)。
[评测记录](https://www.luogu.com.cn/record/51595951)，看起来常数不算太大（雾）。
```cpp
typedef const int cint;
typedef const double cdb;
cint MAXN=2e6+5;
int n,m;
cdb PI=acos(-1);
cdb PI2=2*PI;
int rev[MAXN<<1];
struct complex //STL 其实有这个库的，但众所周知常数大（吧
{
double a,b;
complex(){}
complex(cdb& a,cdb& b):a(a),b(b){}
complex operator+(const complex& x)const {return complex(a+x.a,b+x.b);}
complex operator-(const complex& x)const {return complex(a-x.a,b-x.b);}
complex operator*(const complex& x)const {return complex(a*x.a-b*x.b,a*x.b+b*x.a);}
}f[MAXN<<1],g[MAXN<<1];
void FFT(complex *f,cint& Len,cint& inv)
{
if(Len==1) return;
cint len=Len/2;
FFT(f,len,inv);
FFT(f+len,len,inv);
const complex tg(cos(PI2/Len),sin(PI2/Len)*inv);
complex b(1,0);
for(int i=0;i<len;i++)
{
const complex& t=b*f[i+len];
f[i+len]=f[i]-t;
f[i]=f[i]+t;
b=b*tg;
}
}
int main()
{
read(n,m);
for(int i=0;i<=n;i++) read(f[i].a);
for(int i=0;i<=m;i++) read(g[i].a);
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
for(int i=0;i<n;i++)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int i=0;i<n;i++)
if(i<rev[i]) swap(g[i],g[rev[i]]);
FFT(f,n,1);
FFT(g,n,1);
for(int i=0;i<n;i++) f[i]=f[i]*g[i];
for(int i=0;i<n;i++)
if(i<rev[i]) swap(f[i],f[rev[i]]);
FFT(f,n,-1);
for(int i=0;i<=m;i++) printf("%d ",int(f[i].a/n+0.5));
printf("\n");
return 0;
}
```
# Optimization
实际上 FFT 有很多种写法和优化方法，有些常数差异还比较大。
# Reference
- [OI-Wiki](https://oi-wiki.org/math/poly/fft/)。
- [command_block 的博客。](https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji)