如何成为卷王（上）
云浅知处
·
·
个人记录
卷积&多项式乘法
想要成为卷王，首先需要学会卷积。
两个数列 \bold a 和 \bold b 的卷积 \bold c 定义为：
c_k=\sum_{i+j=k}a_ib_j
类似地定义循环卷积：
c_k=\sum_{(i+j)\bmod n=k}a_ib_j
那么两个多项式 A(x)=\sum\limits_{i=0}^{n-1}a_ix^i 和 B(x)=\sum\limits_{i=0}^{n-1}b_ix^i 的乘积定义为：
C(x)=A(x)\cdot B(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_ib_jx^{i+j}=\sum_{i=0}^{2(n-1)}c_ix^i
其中，\bold c 是 \bold a 和 \bold b 的卷积。因此「多项式乘法」和「计算两序列的卷积/循环卷积」是等价的。
如果直接暴力计算的话，不难发现复杂度是 O(n^2) 的。
这太慢了，卷王需要速度。
怎么优化这个 n^2 的复杂度啊？
用点值表示法试试！
如果用点值表示法，我们发现，如果知道了 C(x)=A(x)B(x) 的足够多的点值，那么就可以计算出 C(x) 的各项系数。
如果朴素地算出这些点值，需要在 O(n^2) 的时间>_<
所以我们优化了个寂寞啊
至少要低于 O(n^2) 吧qwq......
这里我们考虑一组特殊的点值。
单位根
我们定义 \sqrt{-1}=i，即 i^2=-1。
方程 x^n-1=0 的 n 个复根 x 被称为 n 次单位根，记作 \omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}。
其中，\omega_n^k 被定义为 e^{\frac{2k\pi i}{n}}=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}。而 e^{i\theta}=\cos\theta+i\sin\theta 就是欧拉公式，它的证明和 \text{FFT} 没啥关系，就不写啦>_<
例如，$3$ 次单位根有：$\omega_3^0=1,\omega_3^1=\dfrac{\sqrt{3}i-1}{2},\omega_3^2=\dfrac{-\sqrt{3}i-1}{2}$。
我们需要用到单位根的这两个性质：
$$
\omega_{n}^{k}=\omega_{nd}^{kd},d\in \mathbb N\tag{1}
$$
事实上，$\omega_{nd}^{kd}=\cos\dfrac{2kd\pi}{nd}+i\sin\dfrac{2kd\pi}{nd}=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}=\omega_n^k$，故 $(1)$ 成立。
$$
\omega_{2n}^{n+k}=-\omega_{2n}^k\tag{2}
$$
事实上，$\omega_{2n}^{n+k}=\omega_{2n}^n\cdot\omega_{2n}^k=\omega_2^1\cdot\omega_{2n}^k=-\omega_{2n}^k$，故 $(2)$ 成立。
-----------------------
## 离散傅里叶变换
设 $\bold a$ 是长为 $n$ 的序列，对于 $k\in [0,n)$，定义：
$$
b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}
$$
那么我们称 $\bold b$ 为 $\bold a$ 的**「离散傅里叶变换」**，即 $\text{DFT}$ ，记作 $\bold b=\mathcal{F}(\bold a)$ 。
如果我们将序列 $\bold a$ 与 $\bold b$ 的循环卷积记作 $\bold c=\bold a\ast \bold b$，那么不难发现：
$$
\mathcal{F}(\bold c)=\mathcal{F(\bold a\ast\bold b)=}\mathcal{F(\bold a)}\cdot\mathcal{F(\bold b)}
$$
其中 $\cdot$ 表示逐点乘积。
如果我们令多项式 $A(x)=\sum a_ix_i$，不难发现 $b_k$ 就是 $A(x)$ 在 $\omega_n^k$ 处的点值，即：$A(\omega_n^k)=b_k$。
也就是说，想要计算 $\bold a$ 的 $\text{DFT}$，那么我们就需要 $A(x)=\sum a_ix^i$ 的 $n$ 个点值。
同样地，只要我们知道了 $\bold a$ 的 $\text{DFT}$，我们同样可以利用这 $n$ 个点值表示出 $A(x)$，计算出 $A(x)$ 的各项系数。
但是，如果朴素地计算 $\text{DFT}$ 也就是这 $n$ 个点值，那复杂度是 $O(n^2)$ 的>\_<
于是就出现一个问题：如何快速计算一个序列 $\bold a$ 的 $\text{DFT}$ 呢？
------------------------------
## 蝴蝶操作
我们考虑将多项式 $A(x)=\sum\limits_{i=0}^{n-1}$ 的 $n$ 项系数按奇偶性分类：
不妨设 $n=2m$，则
$$
\begin{aligned}
A(x)&=\sum_{i=0}^{n-1}a_ix^i=\sum_{i=0}^{m-1}a_{2i}x^{2i}+\sum_{i=0}^{m-1}a_{2i+1}x^{2i+1}\\
&=\sum_{i=0}^{m-1}a_{2i}x^{2i}+x\sum_{i=0}^{m-1}a_{2i+1}x^{2i}
\end{aligned}
$$
我们记
$$
A_0(x)=\sum_{i=0}^{m-1}a_{2i}x^{i}\\
A_1(x)=\sum_{i=0}^{m-1}a_{2i+1}x^{i}
$$
那么就得到了两个次数不超过 $m$ 的多项式 $A_0(x)$ 和 $A_1(x)$，并且有：
$$
A(x)=A_0(x^2)+xA_1(x^2)\tag{\#}
$$
现在，如果我们有了 $A_0(\omega_m^k)$ 和 $A_1(\omega_m^k)$ 的值，带入 $(\#)$ 式，计算得：
$$
A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_{n}^kA_1(\omega_n^{2k})=A_0(\omega_m^{k})+\omega_{n}^kA_1(\omega_m^k)\\
A(\omega_n^{k+m})=A_0(\omega_n^{2k+n})+\omega_{n}^{k+m}A_1(\omega_n^{2k+n})=A_0(\omega_m^{k})-\omega_{n}^kA_1(\omega_m^{k})
$$
因此，只要我们有了 $A_0(x)$ 和 $A_1(x)$ 在 $\omega_m^0,\omega_m^1,\cdots,\omega_m^{m-1}$ 处的点值，就可以在 $O(n)$ 的时间内计算出 $A(x)$ 在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 处的点值。
那么我们递归下去计算 $A_0(x)$ 和 $A_1(x)$ 的点值，每次将问题规模减半，合并子问题又是 $O(n)$ 的，于是复杂度就是：
$$
T(n)=2T(n/2)+O(n)=O(n\log n)
$$
于是，我们成功在 $O(n\log n)$ 的时间内计算出了 $A(x)$ 在 $n$ 次单位根上的点值，也就是其系数数列 $\bold a$ 的 $\text{DFT}$ 。
这个过程通常被称为**「蝴蝶操作」**。
--------------------
现在我们就可以在 $O(n\log n)$ 的时间内计算出 $\bold a$ 和 $\bold b$ 的 $\text{DFT}$。
再根据 $\mathcal{F}(\bold a\ast \bold b)=\mathcal{F(\bold a)}\cdot\mathcal{F(\bold b)}$，我们可以在 $O(n)$ 的时间内计算出它们的卷积的 $\text{DFT}$。
但是我们要算的是 $\bold a\ast \bold b$ 而不是它的 $\text{DFT}$ 啊=\_=
如果把 $\text{DFT}$ 看作点值，然后用拉格朗日插值，把 $\mathcal{F(\bold a\ast\bold b)}$ 转化为 $\bold a\ast\bold b$，那复杂度还是 $O(n^2)$ >\_<
那现在的问题就是：我们知道了序列 $\bold c$ 的 $\text{DFT}$，如何在低于 $O(n^2)$ 的时间内计算出 $\bold c$？
-----------------
## 离散傅里叶变换的逆变换
对于长为 $n$ 的序列 $\bold a$，若其 $\text{DFT}$ 为 $\bold b$，即：
$$
b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\tag{1}
$$
我们证明 $(1)$ 的逆变换，即其 $\text{IDFT}$ 为：
$$
a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{2}
$$
我们将 $(1)$ 中 $b_k$ 的表达式带入 $(2)$，得到：
$$
\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}=\dfrac{1}{n}\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_i\omega_n^{ij}
$$
交换一下求和顺序（其实就是提取公因式），得到：
$$
\dfrac{1}{n}\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}=\dfrac{1}{n}\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}
$$
我们考虑后面的那个小可爱 $\sum\limits_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}
当 j=k 时，有：
\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\sum_{i=0}^{n-1}\omega_n^{ik-ik}=\sum_{i=0}^{n-1}1=n
当 j\neq k 时，则有：
\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\sum_{i=0}^{n-1}\omega_n^{i(j-k)}
我们发现这东西是个类似于等比数列求和的形式，于是用一下等比数列求和公式：
\sum_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum_{i=0}^{n-1}\left(\omega^{j-k}\right)^i=\dfrac{1-\left(\omega^{j-k}_{n}\right)^n}{1-\omega_n^{j-k}}=\dfrac{1-\left(\omega_n^n\right)^{j-k}}{1-\omega_n^{j-k}}=\dfrac{1-1}{1-\omega_n^{j-k}}=0
于是只要 j\neq k，\dfrac{1}{n}a_j\sum\limits_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_{n}^{-ki} 整个就是个 0，不会对答案产生任何影响。
因此这个式子就被化简成了：
\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}=\dfrac{1}{n}\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\dfrac{1}{n}a_k\cdot n=a_k
即
a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{2}
证毕。
类似地，我们也可以证明 (1) 和 (2) 是等价的。
我们称 (2) 为「离散傅里叶变换的逆变换」，即 \text{IDFT}。如果朴素地按照表达式计算 \text{IDFT}，那么复杂度还是 O(n^2) 的。
但实际上，观察 \text{IDFT} 的表达式，对比一下 \text{DFT} 的表达式：
b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\tag{DFT}
a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{IDFT}
实际上只需要将计算 \text{DFT} 时用到的「蝴蝶操作」中的 \omega_n 全部替换成 \omega_n^{-1} ，类似地递归下去，同样可以在 O(n\log n) 的时间内计算 \text{IDFT}。唯一不同的是算完之后需要乘上 \dfrac{1}{n} 。
所以 \text{IDFT} 同样可以在 O(n\log n) 的时间内被计算。
总结
总结一下我们的 \text{FFT} 流程：
- 利用 \text{DFT}，在 O(n\log n) 的时间内分别计算出 \bold a 和 \bold b 的 \text{DFT}。
- 把算出来的出来的这两个 \text{DFT} 相乘，在 O(n) 的时间内得到 \bold c 的 \text{DFT}。
- 利用 \text{IDFT}，在 O(n\log n) 的时间内，把 \bold c 的 \text{DFT} 转化为 \bold c ，就计算出了 \bold a 和 \bold b 的卷积。
综上，我们在 O(n\log n) 的时间内计算出了两个序列的卷积！
FFT 代码实现&注意事项&技巧（递归版）
先写 \text{DFT}：
typedef complex<double> comp;
#define pi acos(-1)//不要手打 pi QAQ，最好用 acos(-1)
void dft(comp *a,int p){
if(p==1)return;//递归边界
int q=p>>1;
comp a1[q],a2[q];
for(int i=0;i<q;i++){
a1[i]=a[i<<1];
a2[i]=a[(i<<1)|1];
}//按奇偶分类
dft(a1,q);dft(a2,q);//分治
for(int i=0;i<q;i++){
comp omega=comp(cos(2*pi*i/p),sin(2*pi*i/p));
a[i]=a1[i]+omega*a2[i];
a[i+q]=a1[i]-omega*a2[i];//照搬我们推出来的式子qwq
}
}
然后把 \text{DFT} 中的 \omega_n 换成 \omega_n^{-1}，写出 \text{IDFT}：
void idft(comp *a,int p){
if(p==1)return;//递归边界
int q=p>>1;
comp a1[q],a2[q];
for(int i=0;i<q;i++){
a1[i]=a[i<<1];
a2[i]=a[(i<<1)|1];
}//按奇偶分类
idft(a1,q);idft(a2,q);//分治
for(int i=0;i<q;i++){
comp omega=comp(cos(-2*pi*i/p),sin(-2*pi*i/p));//添一个负号表示负一次幂qwq
a[i]=a1[i]+omega*a2[i];
a[i+q]=a1[i]-omega*a2[i];//照搬我们推出来的式子qwq
}
}
我们发现 idft() 和 dft() 两者之间只差了一个负号，因此可以合并成：
void fft(comp *a,int p,int tag){//tag=1 表示 dft，tag=-1 表示 idft 不乘 1/n 之前的操作
if(p==1)return;//递归边界
int q=p>>1;
comp a1[q],a2[q];
for(int i=0;i<q;i++){
a1[i]=a[i<<1];
a2[i]=a[(i<<1)|1];
}//按奇偶分类
fft(a1,q,tag);fft(a2,q,tag);//分治
for(int i=0;i<q;i++){
comp omega=comp(cos(tag*2*pi*i/p),sin(tag*2*pi*i/p));
a[i]=a1[i]+omega*a2[i];
a[i+q]=a1[i]-omega*a2[i];//照搬我们推出来的式子qwq
}
}
由于计算 cos() 和 sin() 都需要相对较为高昂的代价，因此我们提前算出 \omega_n=\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}，然后每次循环乘一遍。
于是：
void fft(comp *a,int p,int tag){//tag=1 表示 dft，tag=-1 表示 idft
if(p==1)return;//递归边界
int q=p>>1;
comp a1[q],a2[q];
for(int i=0;i<q;i++){
a1[i]=a[i<<1];
a2[i]=a[(i<<1)|1];
}//按奇偶分类
fft(a1,q,tag);fft(a2,q,tag);//分治
comp omega=comp(1,0),wn=comp(cos(tag*2*pi/p),sin(tag*2*pi/p));
for(int i=0;i<q;i++){
a[i]=a1[i]+omega*a2[i];
a[i+q]=a1[i]-omega*a2[i];//照搬我们推出来的式子qwq
omega=omega*wn;
}
}
注意到如果我们按照上述的操作进行 \text{FFT}，只能处理长度为 2^k 的数列的卷积，其中 k\in \mathbb N^*。
因此，我们一般会选取一个最小的 2^k\ge n+m 作为序列长度，多出来的地方就看做 0，然后再进行 \text{FFT}。
所以，实际上，在计算卷积的时候，我们通常会：
int k=1;for(;k<=(N+M);k<<=1);
然后再：
fft(x,k,1);
fft(y,k,1);
for(int i=0;i<=k;i++)z[i]=x[i]*y[i];
fft(z,k,-1);
最后才是：
for(int i=0;i<=N+M;i++)cout<<(int)(z[i].real()/k+0.5)<<" ";//不要忘记除以序列长度喵 qwq
此时，由于我们寻找的是最小的 2^k\ge n+m ，因此数组并不能直接开成 10^6 这么大，而要开到最小的 2^k\ge 10^6+10^6 ，大概是 2^{21}=2097152 qwq。
方便起见，开成 2100000 也是没有问题的QAQ。
整个代码十分简洁易懂（所以窝并没有加多少注释），就是把我们之前推导的式子照搬过来了qwq。
迭代实现 FFT
然而递归版常数大到爆炸，尤其是我们递归的时候居然还要递归一下开一遍内存，很容易空间爆炸>_<
（虽然递归版开 O2 就能直接艹过板子题）
而且，\text{FFT} 作为一个经常要使用的算法，如果它的常数过大导致整个程序时间爆炸是很不好的qwq。
能不能把 \text{FFT} 改写成迭代形式啊？
在 \text{FFT} 递归树中，如果我们能够快速地处理出叶子节点所构成的序列 A（本例中为 A=\{a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7\}，那么就可以自底向上地实现整个 \text{FFT} 操作。
具体地，我们现在有了一个序列 A，就是这个递归树的叶子节点所构成的序列。
首先我们从前往后，两个两个地取出 A 中元素，然后对于每两个元素，算出它们的 \text{DFT}，这样就有了叶子节点上一层的 n/2 个 \text{DFT}。
然后用这 n/2 个 \text{DFT} ，可以再算出来再往上一层的 n/4 个 \text{DFT}，直到最后算出来递归树的「根」处的 整个序列的 \text{DFT} 。
假设我们已经有了这个序列 A 的初始值，那么就能写出这样的代码：
typedef complex<double> comp;
comp A[MAXN];
void dft(comp *a,int N){ //N 是 2 的正整数次幂。
Pre(a,N);
//假装我们已经通过一些神奇操作在低于 O(n log n) 的时间内预处理出了 A[] qwq。
int n=log2(N);
for(int i=1;i<=n;i++){
int m=1<<i; //m=2^i，就是递归时这一层中每组的长度
for(int j=0;j<N;j+=m){
for(int k=0;k<m/2;k++){
comp omega=comp(cos(2*k*pi/m),sin(2*k*pi/m));
comp x=omega*A[j+k+m/2],y=A[j+k];
A[j+k]=y+x;A[j+k+m/2]=y-x;
}
}
}
}
所以，到底怎么在一个较好的时间复杂度内确定 A 啊？
位逆序置换
我们仔细观察一下递归的过程，还以 7 次多项式的八个系数为例：
- 一开始，这八个系数好好地放在那里：\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}。
- 我们把这个系数按照奇偶分了个类：\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\}。
- 又把每一个按照奇偶分了个类：\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\}。
- 最后，变成了：\{a_0\},\{a_4\},\{a_2\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\}。
观察 (0,4,2,6,1,5,3,7) 的二进制表示：(000,100,010,110,001,101,011,111)
诶，把这个二进制表示翻转过来：
000 100 010 110 001 101 011 111
翻转后: 000 001 010 011 100 101 110 111
这，这不就是 (0,1,2,3,4,5,6,7) 的二进制表示嘛！
实际上这其中的道理是很好解释的：考虑我们分组时的分组依据。
第一次我们将奇偶分了个组，把偶数放到前面，把奇数放到后面，也就相当于把二进制表示中末尾为 0 的放到前面，末位为 1 的放到后面。
再想想把 (0,1,2,3,4,5,6,7) 的二进制表示：显然前面四个的首位是 0，后面四个的首位是 1，毕竟这个序列是递增的嘛QAQ。
也就是说，我们在末位上把 0 全都提到了前头，把 1 全都扔到了后头，而 (0,1,2,3,4,5,6,7) 则是在首位上把 0 提到前头，把 1 扔到后头，一模一样啊QAQ！
那同理，我们后面的操作也都是把每组中的 0 提到前头，把 1 扔到后头；这恰巧相当于 (0,1,2,3,4,5,6,7) 这一顺序排列在每组中的操作。
如果我们令 \text{rev}(n) 为 「将 n 的二进制表示翻转过来之后得到的值」，根据上面的结论，我们就有：
A_{\text{rev}(i)}=a_i
那我们就可以这么处理序列 A
for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev(i)]);
其中 rev(n) 是一个计算「将 n 的二进制翻转后得到的值」的函数。
这样一来，经过处理后的序列 \bold a 就是我们要求的叶子节点了。
计算 rev(n) 显然可以做到 O(\log n)，毕竟 n 在二进制下的位数也一共只有 O(\log n) ，随便你怎么搞都行嘛qwq。
这样的话就做到了 O(n\log n) 计算序列 A。
实际上，如果 n 在二进制表示下的后 L 位（本例中 L=3），那么将 n 的这 L 位翻转，其实就相当于：把前面的 L-1 位拽出来，翻转一下，再把最后一位搞到最前头去qwq！
也就是说，如果我们知道了 \text{rev}(\texttt{n>>1}) 即 \text{rev}(n/2) 的值，那么就可以直接 O(1) 算出来 \text{rev}(n) 的值。
来构造一下这个递推式吧qwq。
首先要把前 L-1 位单拉出来并翻转：rev[i>>1]。
注意，这里我们相当于把前面的 L-1 全部右移一位，把最后一位「顶掉」，这会导致最前面空出来一个 0。
那翻转之后这个最前面的 0 就被翻到了最后一位，也就是说，我们需要再右移一位，把这个 0 顶掉，即：rev[i>>1]>>1。
现在，这个 rev[i>>1]>>1 就是 rev[i] 的后面 L-1 位翻转得到的结果了，而且位置是没有问题的qwq。
现在只需要把 i 的最后一位搞到前头去就行啦！
i 的最后一位，就是 i&1；搞到前头去就相当于左移 L-1 位（这样它变成了第一位），再和刚才的 rev[i>>1]>>1 「或」一下即可qwq。
那么递推式就是：rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1))
因此，我们就利用递推做到了 O(n) 预处理 \text{rev}(i)，然后 O(n) 计算序列 A qwq！
FFT 代码实现&注意事项&技巧（非递归版）
此时我们需要算出这些数在二进制下的位数，也就是把这一步改一下：
int k=1,L=0;
for(;k<=(n+m);k<<=1,L++)//L 就是我们需要考虑的位数。
其他的都基本没有什么区别qwq。
typedef complex<double> comp;
int rev[MAXN],k,L;
void init(){
rev[0]=0;
for(int i=0;i<k;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
}//预处理 rev[]
void fft(comp *a,int N,int tag){
init();
for(int i=1;i<N;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
int n=log2(N);
for(int i=1;i<=n;i++){
int m=1<<i; //m=2^i，就是递归时这一层中每组的长度
comp wn=comp(cos(tag*2*pi/m),sin(tag*2*pi/m));
for(int j=0;j<N;j+=m){
comp omega=comp(1,0);
for(int l=0;l<m/2;l++){
comp x=omega*a[j+l+m/2],y=a[j+l];
a[j+l]=y+x;a[j+l+m/2]=y-x;
omega=omega*wn;
}
}
}
}