2025 新看法

当时写这篇文章是为了刚上大一时的自己，现在笔者感觉有点 outdate 了，可以用更加优雅的方式去阐述为什么数和数之间的乘法构成的循环卷积，可以通过映射到复数域/有限域上代入单位根后点积做。

代数角度理解

从代数角度上看，两个数作为按照以 $\{10^x\}, x\in \mathbb{Z}$ 为基底的向量分解，从而有数到多项式的同构 $f: \mathbb{Z}\backslash\{0\} \to \mathbb{Z\backslash\{0\}}[x]$。刨去零在于特定的多项式到零是满射。也就是可以通过不断堆高高次项，同时保持高次项系数置零来找出无穷多的零多项式。

读者可以找模论(module theory)来看看，材料综合加上手实操永远是最好的方法。

但如果暂且忽略这一点以及忽略进位时，不妨先考虑两个数乘起来发生了什么。

$$\begin{matrix} &&&1 &2 &3 &4 &5 &6 \\ &&&\times & &4 &3 &2 &1 \\ \hline &&&1 &2 &3 &4 &5 &6 \\ &&2 &4 &6 &8 &10 &12 \\ &3 &6 &9 &12 &15 &18 \\ 4 &8 &12 &16 &20 &24 \\ \hline 4 &11 &20 &40 &40 &50 &32 &17 &6 \end{matrix}$$

感性认识上是计算得出的结果经过某种映射作用之后向左移位。

$\mathbb{Z}/p\mathbb{Z}$ 作用。

现在把数加载到多项式有
$$\begin{matrix} &&&x^5 &2x^4 &3x^3 &4x^2 &5x &6 \\ &&&\otimes & &4x^3 &3x^2 &2x &1 \\ \hline &&&1x^5 &2x^4 &3x^3 &4x^2 &5x &6 \\ &&2x^6 &4x^5 &6x^4 &8x^3 &10x^2 &12x \\ &3x^7 &6x^6 &9x^5 &12x^4 &15x^3 &18x^2 \\ 4x^8 &8x^7 &12x^6 &16x^5 &20x^4 &24x^3 \\ \hline 4x^8 &11x^7 &20x^6 &40x^5 &40x^4 &50x^3 &32x^2 &17x &6 \end{matrix}$$

观察发现，如果将这些多项式按照 $x^p=1$ 的限制来约束在多项式环上，那么就有 $x^p-1\equiv 0\pmod{x^p-1}$。对于不够的项，补零到最高次项 $x^{n-1}$ 即可。从而 $x^p\equiv 1\pmod{x^p-1}$，进而有

$$\begin{align} x^{p-1}(b_0+b_1x+\cdots+b_{p-1}x^{p-1}) &\equiv b_0x^{p-1}+b_1x^{p}+\cdots+b_{p-1}x^{2p-2} \pmod{x^p-1} \\ &\equiv b_1+\cdots+b_{p-1}x^{p-2}+b_0x^{p-1} \pmod{x^p-1} \\ \\ x^{\lambda}(c_0+c_1x+\cdots+c_{p-1}x^{p-1}) &\equiv c_{p-\lambda}+c_{p-\lambda+1}x\cdots+c_{0}x^{\lambda}+\cdots_+c_{p-\lambda-1}x^{p-1} \pmod{x^p-1} \\ \\ \therefore b_0 + b_1x + \cdots + b_{p-1}x^{p-1} + c_{p}x^p + \cdots + c_{m}x^m &\equiv (b_0 + c_p) + (b_1 + c_{p+1})x + \cdots + (b_{m \bmod p} + c_{m})x^{m\bmod p} + \ldots + b_{p-1}x^{p-1} \pmod{x^p-1} \end{align}$$

从这一无穷多个整数多项式打到有限多个整数多项式形成的满射出发，于是就构建了同态映射 $g: \mathbb{Z}[x] \to \mathbb{Z}[x]/(x^p-1)$，用于生成循环左移后的元素以及在元素环内的累加操作。同时，所有满足 $x^p=1$ 的解 $x$ 都能作为循环左移的算子。

另外，如果从 $\mathbb{Z}[x]/(x^p-1)$ 选取某一特定对象以及其所有循环左移的像，来构成新的向量，那么也不难发现：

1. 整数到这一向量的映射是能按向量的第一项分量对应起来的；
2. 向量的生成相当于受到某个特定矩阵的作用。

从而有下列对应。

$$\begin{matrix} \mathbb{Z}_n &\times &\mathbb{Z}_m &\longrightarrow &\mathbb{Z}_{n+m-1} &n + m - 1 < p, p\in \mathbb{P}\\ a & &b & &\boxed{ab=ba}\\ \Bigg\downarrow & &\Bigg\downarrow &&\Bigg\downarrow \\ \mathbb{Z}[x]/(x^p-1) &\times &\mathbb{Z}[x]/(x^p-1) &\longrightarrow &\mathbb{Z}[x]/(x^p-1)\\ \displaystyle\sum_{i=0}^{p-1}\alpha_ix^i & &\displaystyle\sum_{j=0}^{p-1}\beta_jx^j & &\boxed{\displaystyle\sum_{0\leqslant i,j< p}^{p-1}\alpha_j\beta_{i-j}x^i}\\ f(x)=\alpha_0 + \alpha_1x + \cdots + \alpha_{p-1}x^{p-1} &&g(x) &&\boxed{f(x)g(x) = g(x)f(x)}\\ \Bigg\downarrow & &\Bigg\downarrow &&\Bigg\downarrow \\ \mathbb{Z}[\mathbf{J}_{p\times p}] &\times &\mathbb{Z}^p &\longrightarrow &\mathbb{Z}^p \\ \\ \begin{pmatrix} \alpha_0 &\alpha_{p-1} &\cdots &\alpha_2 &\alpha_1 \\ \alpha_1 &\alpha_0 &\ddots &\alpha_3 &\alpha_2 \\ \vdots &\ddots &\ddots &\ddots &\vdots \\ \alpha_{p-2} &\alpha_{p-3} &\ddots &\alpha_0 &\alpha_{p-1}\\ \alpha_{p-1} &\alpha_{p-2} &\cdots &\alpha_1 &\alpha_0 \end{pmatrix} &&\begin{pmatrix}\beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-2} \\ \beta_{p-1}\end{pmatrix} &&\Bigg\downarrow\\ \\ \mathbb{Z}[\mathbf{J}_{p\times p}] &\times &\mathbb{Z}[\mathbf{J}_{p\times p}] &\longrightarrow &\mathbb{Z}[\mathbf{J}_{p\times p}] \\ f(\mathbf{J})=\alpha_0 + \alpha_1\mathbf{J} + \cdots + \alpha_{p-1}\mathbf{J}^{p-1} &&g(\mathbf{J}) &&\boxed{f(\mathbf{J})g(\mathbf{J}) = g(\mathbf{J})f(\mathbf{J})}\\ \mathbf{A}=\begin{pmatrix} \alpha_0 &\alpha_{p-1} &\cdots &\alpha_2 &\alpha_1 \\ \alpha_1 &\alpha_0 &\ddots &\alpha_3 &\alpha_2 \\ \vdots &\ddots &\ddots &\ddots &\vdots \\ \alpha_{p-2} &\alpha_{p-3} &\ddots &\alpha_0 &\alpha_{p-1}\\ \alpha_{p-1} &\alpha_{p-2} &\cdots &\alpha_1 &\alpha_0 \end{pmatrix} & &\mathbf{B}=\begin{pmatrix} \beta_0 &\beta_{p-1} &\cdots &\beta_2 &\beta_1 \\ \beta_1 &\beta_0 &\ddots &\beta_3 &\beta_2 \\ \vdots &\ddots &\ddots &\ddots &\vdots \\ \beta_{p-2} &\beta_{p-3} &\ddots &\beta_0 &\beta_{p-1}\\ \beta_{p-1} &\beta_{p-2} &\cdots &\beta_1 &\beta_0 \end{pmatrix} &&\boxed{\mathbf{A}\mathbf{B}=\mathbf{B}\mathbf{A}} \end{matrix}$$
其中 $\mathbf{J}$ 为满秩的循环矩阵，满足

$$\mathbf{J} = \begin{pmatrix} 0 &0 &\cdots &0 &1 \\ 1 &0 &\cdots &0 &0 \\ 0 &1 &\cdots &0 &0 \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 &0 &\cdots &1 &0 \end{pmatrix}_{p\times p} = \left(\begin{array}{c|c}\boldsymbol{0} &1 \\ \hline \mathbf{I}_{(p-1)\times (p-1)} &\boldsymbol{0}\end{array}\right), \mathbf{J}^p = \mathbf{E}$$

而构成了乘法过程中在矩阵意义下的卷积作用。

格外注意框起来的部分，进而得到数是向量，是多项式，是模过了 $x^p-1$ 后的多项式，是模过之后构成的向量组/矩阵的理解。此外，将单个向量转换为向量组是从卷积矩阵乘法可交换出发得到的。

另外又能注意到同态核 $\ker(x^p-1)$ 也就是满足 $x^p=1$ 的所有根，如果放到 $\mathbb{C}$ 来看，就是按角度均匀分掉了整个单位圆上的根(分圆)。还能注意到这个分圆多项式本身并不是可约的，而是可以由多个在实数域上不可约多项式（分圆多项式）通过多项式乘法计算得到。即满足：

$$\begin{align} x^p-1 &= \prod_{d \mid p}^p \Phi_d(x) \\ \Phi_d(x) &= \prod_{k\in \mathbb{Z}_d} (x - e^{\frac{2\pi ik}{d}}) % \backslash\{0\} \end{align}$$

但这并不妨碍它作为模来使用。

注意到 $p\in \mathbb{P}$ 时，$\Phi_p(x)=\dfrac{x^p-1}{x-1} = 1 + x +\cdots + x^{n-1}$，此时在整系数域或有理数域上不可约。

往后的探讨就可以从分圆扩张出发，但此处仍借助循环左移的视角——

$$\begin{pmatrix} \alpha_0 &\alpha_{p-1} &\cdots &\alpha_2 &\alpha_1 \\ \alpha_1 &\alpha_0 &\ddots &\alpha_3 &\alpha_2 \\ \vdots &\ddots &\ddots &\ddots &\vdots \\ \alpha_{p-2} &\alpha_{p-3} &\ddots &\alpha_0 &\alpha_{p-1}\\ \alpha_{p-1} &\alpha_{p-2} &\cdots &\alpha_1 &\alpha_0 \end{pmatrix} \begin{pmatrix} 1 \\\omega^{p-1} \\ \cdots \\\omega \end{pmatrix} = \begin{pmatrix} f(\omega) \\ \omega^{p-1} f(\omega) \\ \vdots \\ \omega f(\omega) \end{pmatrix} = f(\omega)\begin{pmatrix} 1 \\ \omega^{p-1} \\ \vdots \\ \omega \end{pmatrix}$$

于是，就证明了 $f(\omega)$ 是卷积矩阵的特征值，$(1,\omega^{p-1},\ldots,\omega)^T$ 为其中的一个特征向量。注意到 $\langle\{1,\omega,\ldots,\omega^{p-1}\}, \times\rangle\cong\mathbb{Z}/p\mathbb{Z}$，再用同样的手法作用所有 $w \in \ker(x^p-1)$，从而可以得到

$$\mathbf{A}w = f(\omega)w \Rightarrow \mathbf{A}\sim\begin{pmatrix} f(1) & & & \\ &f(\omega) & & \\ & &\ddots \\ & & &f(\omega^{p-1}) \end{pmatrix} \Rightarrow \mathbf{A}\mathbf{B}=\mathbf{B}\mathbf{A}\sim\begin{pmatrix} f(1)g(1) & & & \\ &f(\omega)g(\omega) & & \\ & &\ddots \\ & & &f(\omega^{p-1})g(\omega^{p-1}) \end{pmatrix}$$

至此可以有两种视角看待乘法：

• 整数的十进制分量被加载到整多项式环中，其乘法为多项式环上的卷积。
• 整数的十进制分量被加载到整系数多项式循环矩阵环上，此类矩阵的行列受模 $p$ 循环群作用，且生成元可为任意一行列，整数的乘法等价于矩阵环中的乘法，两者均可交换；
  • 傅里叶变换在代数上看就是通过对角化的手段，将卷积变为原有整数的十进制分量构成的多项式对 $x^p-1$ 的根所映射到的值的点乘。（即所谓时域卷积，频域相乘）。
  • 如果继续深入下去，实际上整个过程是借助对十进制分量所在的指标产生某种映射而形成的。感兴趣的可以去学群表示（主要是看特征标）和傅里叶分析课，也可以手动玩一玩二维FFT，此处不过多展开。

进而可以运用到有限域中，构成 NTT 的基础¹′²。
进而可以在某些条件下推广到 $\mathbb{R}$ 以及 $\mathbb{R}^{\infty}$ 上，构成连续傅里叶变换。

结合前面提及的特征向量组，现不难直接列出所有的特征向量。

$$\mathbf{F} = \begin{pmatrix} 1 &1 &\cdots &1 \\ 1 &\omega^{p-1} &\cdots &\omega \\ \vdots &\vdots &\ddots &\vdots \\ 1 &\omega &\cdots &\omega^{p-1} \end{pmatrix}$$

又由于该矩阵现处于复数域上，故为复矩阵。考虑其共轭阵 $\mathbf{F}^H$，则

$$\mathbf{F}^H = \begin{pmatrix} 1 &1 &\cdots &1 \\ 1 &\omega^{-p+1} &\cdots &\omega^{-1} \\ \vdots &\vdots &\ddots &\vdots \\ 1 &\omega^{-1} &\cdots &\omega^{1-p} \end{pmatrix} = \begin{pmatrix} 1 &1 &\cdots &1 \\ 1 &\omega &\cdots &\omega^{-1} \\ \vdots &\vdots &\ddots &\vdots \\ 1 &\omega^{-1} &\cdots &\omega \end{pmatrix}$$

注意到，

$$f_i^Hf_j = \begin{cases} p & i =j \\ \\ 1+ \omega^{i+j} + \omega^{2(i+j)} + \cdots + \omega^{(p-1)(i+j)} = \frac{\omega^{p(i+j)} - 1}{\omega^{i+j} - 1} = \boxed{\frac{1-1}{\omega^{i+j} - 1}} = 0 &i \ne j \end{cases}$$

事实上如果纳入 $\pmod{p}$ 就得不出单位阵。

这说明 $\mathbf{F}\mathbf{F}^H = \mathbf{F}^H\mathbf{F} = p\mathbf{E}$，且单位化后的 $\mathbf{F}$ 为酉矩阵。总之，此时得到正交的、可参与到对角化两卷积矩阵的傅里叶矩阵。

$$\mathbf{F}^\# = \frac{1}{\sqrt{p}}\begin{pmatrix} 1 &1 &\cdots &1 \\ 1 &\omega &\cdots &\omega^{-1} \\ \vdots &\vdots &\ddots &\vdots \\ 1 &\omega^{-1} &\cdots &\omega \end{pmatrix}$$

从而有

$$(\mathbf{F}^\#)^H\mathbf{AB}\mathbf{F}^\# = \dfrac{\mathbf{F}^H\mathbf{AB}\mathbf{F}}{p} = \mathbf{\Lambda}_{\mathbf{A}}\mathbf{\Lambda}_{\mathbf{B}} = \begin{pmatrix} f(1)g(1) & & & \\ &f(\omega)g(\omega) & & \\ & &\ddots & \\ & & &f(\omega^{-1})g(\omega^{-1}) \end{pmatrix}$$

而如何将整数乘法优化为 $O(n\log n)$ 则在下文介绍。

[1] Music Through Fourier Space
[2] 代数同构视角下的离散 Fourier 变换
[2] Discrete Fourier transformation in the perspective of isomorphism
[3] 信息论与编码
[4] 傅里叶变换&环同构
[5] gröbner-基学习
[6] 「笔记⑤」又双叒叕谈傅里叶变换
[7] 傅里叶分析（七）有限傅里叶分析
[8] 有限群上的傅里叶分析

2022 正文

不难注意到，复数单位根均满足下列几个等式。

$$\begin{align}&\omega^p=1\tag{1} \\ &\omega^m=-\omega^{ m+\frac{p}{2} }\tag{2}\\ &\omega_{2p}^{2m}=\omega_p^m. \tag{3}\\ &\omega= \cos\theta+i\sin\theta,\omega^{-1}=\cos\theta-i\sin\theta \tag{4} \end{align}$$

为了方便，我们将 $(3)$ 式折半运算的记法变为 $\omega^{2^p} \rightarrow \omega_1^{2^{p-1}} \rightarrow \ldots \rightarrow \omega_{p-1}^2\rightarrow \omega_{p}$ 的形式。

中间两个式子与增了一半的单位根、折了一半的单位根有关，直观上是能察觉到这两个性质是能让计算机减少重复计算的。而折半的思路又启发我们，序列的长度需要是 $2$ 的非负整数次幂。所以最少就应该需要 $p = 2^{\left \lceil \log (n+m) \right \rceil }$。

具体的计算可以参考下面几行代码。

    int n1, n2;
    cin >> n1 >> n2;
    int bits=0, m=0;
    while ((1<<bits) < n1+n2) {
        bits ++;
    }
    m = 1 << bits;
/*
 * 又或者是下面一步到位的写法。
 * int bits = -1 , m = 1;
 * while (m < n1+n2) {
     m <<= 1, bits ++;
 }
 */

但关键是怎么用上折半这个条件？

人们研究多项式函数图像以及多项式乘法的时候注意到多项式函数是有奇偶性的。也就是只要知道正半轴上分布的一半的点，负半轴只需要根据多项式函数的奇偶性就能确定。

基于这一点，人们发明出了奇偶划分的做法。

以多项式 $p(x)=3x^5+9x^4+6x^3+x^2+4x+1$ 为例，这个函数可以被划分成

$$p(x)=(9x^4+x^2+1)+x(3x^4+6x^2+4)=p_1(x^2)+xp_2(x^2)$$

这样我们就可以利用一半的点值完成完整序列的计算了。

具体来说，已知 $p(x)=p_1(x^2)+xp_2(x^2)$，即可算出 $p(-x)=p_1(x^2)-xp_2(x^2)$。

假如最后算到的长度为 $p=8$ ，那么对于给定的因数 $f$ 运用奇偶划分，结合单位根的性质，就可以有下面的式子。（注意，这里只是我自己的一家之词，个人认为并不够正式）。

$$\begin{aligned} f(\omega^m)&=a_0+a_1\omega^{m}+a_2\omega^{2m}+a_3\omega^{3m}+a_4\omega^{4m}+a_5\omega^{5m}+a_6\omega^{6m}+a_7\omega^{7m} \\&= \left(a_0+a_2\omega^{2m}+a_4\omega^{4m}+a_6\omega^{6m}\right) + \omega^{m}\left(a_1+a_3\omega^{2m}+a_5\omega^{4m}+a_7\omega^{6m}\right) \\&= A_0(\omega^{2m}) +\omega^{m}A_1(\omega^{2m}) \\ \\&= \left(a_0+a_2\omega_1^{m}+a_4\omega_1^{2m}+a_6\omega_1^{3m}\right) + \omega^{m}\left(a_1+a_3\omega_1^{m}+a_5\omega_1^{2m}+a_7\omega_1^{3m}\right) \\&=A_0(\omega_1^{m}) +\omega^{m}A_1(\omega_1^{m}) \\ \\&= \left((a_0+a_4\omega_1^{2m})+\omega_1^{m}(a_2+a_6\omega_1^{2m})\right) + \omega^{m}\left((a_1+a_5\omega_1^{2m})+\omega_1^{m}(a_3+a_7\omega_1^{2m}) \right) \\&=\left(A_{00}(\omega_1^{2m})+\omega_1^{m}A_{01}(\omega_1^{2m})\right) +\omega^{m}\left(A_{10}(\omega_1^{2m})+\omega_1^{m}A_{11}(\omega_1^{2m}) \right) \\&=\left[A_{00}(\omega_2^{m})+\omega_1^{m}A_{01}(\omega_2^{m})\right] +\omega^{m}\left[A_{10}(\omega_2^{m})+\omega_1^{m}A_{11}(\omega_2^{m})\right] \\ \\&=\left((a_0+\omega_2^{m}a_4)+\omega_1^{m}(a_2+\omega_2^{m}a_6)\right) + \omega^{m}\left((a_1+\omega_2^{m}a_5)+\omega_1^{m}(a_3+\omega_2^{m}a_7) \right) \\&=\left[(A_{000}+\omega_2^{m}A_{001})+\omega_1^{m}(A_{010}+\omega_2^{m}A_{011})\right] +\omega^{m}\left[(A_{100}+\omega_2^{m}A_{101})+\omega_1^{m}(A_{110}+\omega_2^{m}A_{111})\right] \\ \\ \Rightarrow f(\omega^{m+\frac{p}{2}})&=a_0-a_1\omega^{m}+a_2\omega^{2m}-a_3\omega^{3m}+a_4\omega^{4m}-a_5\omega^{5m}+a_6\omega^{6m}-a_7\omega^{7m} \\&= A_0(\omega^{2m}) -\omega^{m}A_1(\omega^{2m}) \\&\ldots \\&=\left[(A_{000}+\omega_2^{m}A_{001})+\omega_1^{m}(A_{010}+\omega_2^{m}A_{011})\right] \boldsymbol{\color{red}-}\omega^{m}\left[(A_{100}+\omega_2^{m}A_{101})+\omega_1^{m}(A_{110}+\omega_2^{m}A_{111})\right] \end{aligned}$$

也就是说，我们可以将其写成以下的形式。

$$\begin{cases} f(\omega^{m})=A_0(\omega^{2m})+\omega^{m}A_1(\omega^{2m})\\ \\ f(\omega^{m+\frac{p}{2}})=A_0(\omega^{2m})-\omega^{m}A_1(\omega^{2m}) \end{cases} \iff \begin{bmatrix}f(\omega^{m}) \\f(\omega^{m+\frac{p}{2}})\end{bmatrix} = \begin{bmatrix}1 &\omega^{m} \\1 &-\omega^{m}\end{bmatrix}\begin{bmatrix}A_0(\omega^{2m}) \\A_1(\omega^{2m})\end{bmatrix}$$

其中 $m\in [0,\dfrac{p}{2})$。所以只要正变换能解决，逆变换自然不在话下。

结合刚刚得到的奇偶划分式，只需要扫完前一半的序列并覆盖写入对应的位置，后一半序列做除了加法不同以外的相同操作，就自然能在当前一层中计算得到。而这便是整个算法最抽象的部分。

对于逆变换，这显然也是能做到的。因为对整个序列做两次奇偶划分相当于两次交换。所以只需要分别对两个序列做好正变换、逐项相乘够 $p$ 项后，施行逆变换再逐项除以序列长度即可。

接下来就是写代码的时光了。

实现

using db = long double;
const db pi=(db)acos(-1);
const int N = 4e5 + 10;
struct cp {
    db x,y;
    cp(db xx=0,db yy=0):x(xx),y(yy){}
    cp operator + (const cp& a){ return cp(a.x+x, a.y+y);}
    cp operator - (const cp& a){ return cp(x-a.x, y-a.y);}
    cp operator * (const cp& a){ return cp(x*a.x-y*a.y, x*a.y+y*a.x);}
};
void fft(cp a[], int len, int sig = 1) {
    if(len == 1) { return; }                                                    /* 既然是划分，划分长度到 1 为止。*/
    int mid = len >> 1;
    static cp f[N];
    for (int i = 0; i < mid; ++ i) { /* 将 a 中的数按奇偶分到 b 里。   */
        f[i] = a[i << 1];
        f[i + mid] = a[i << 1 | 1];
    }
    for (int i = 0; i < len; ++ i) { a[i] = f[i]; }                                  /* 然后覆盖                      */
    fft(a, mid, sig), fft(a + mid, mid, sig);                                   /* 能划分就划分左右子区间         */
    for(int i = 0; i < mid; ++ i) {                                              /* 只用一半循环                  */
        cp w(cos(2 * pi * i / len), sig * sin(2 * pi * i / len));               /* sig 标志正或逆变换            */
        f[i] = a[i] + w * a[i + mid], f[i + mid] = a[i] - w * a[i + mid];       /* 完成当前整个区间的计算         */
    }
    for(int i = 0;i < len; ++ i) { a[i] = f[i]; }                               /* 算完以后要记得覆盖             */
}

[HTML_REMOVED]

然后考虑优化。static cp b[N] 和递归有时候会搞人。所以可以丧心病狂的省略掉它们。具体的做法如下。

注意到算的时候是把原始序列，如 $\{0,1,2,3,4,5,6,7\}$ 转变成 $\{0,4,2,6,1,5,3,7\}$，然后再合并回来。所以可以在处理之前先把这工作做了。然后再整体算回去。

对于变换的规律，可以展开他们的二进制数来看。对于存储数的变换，有：
$$\begin{aligned} \text{[000]} \gets 000 \\ [001] \gets 100 \\ [010] \gets 010 \\ [011] \gets 110 \\ [100] \gets 001 \\ [101] \gets 101 \\ [110] \gets 011 \\ [111] \gets 111 \end{aligned}$$
显然只是二进制数上的反转。那么，可以给出最初的反转方案。

int I=i; // rev[N] 初始化为0
while(i) { rev[I]<<=1, rev[I]|=i&1, i>>=1; }

但，显然不可能对每个数都这么操作一次。我们需要更快的办法。我们可以利用先前求得的 rev[i>>1] 来 dp，毕竟我们计算机在算好了先前反转的块后就能偷懒了。这个时候，就可以有下面的式子。

for(int i = 0;i < m; ++ i) { rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bits - 1)); }

先前算好的 rev[i>>1] 右移一位后，或上当前左移了 bits-1 位的 i&1，即为我们要求的 rev[i] 了。相应地，快速傅里叶变换的函数就能写成下面的样子。

const int N = 4e5 + 10;
cp a[N], b[N];
void fft(cp a[], int len, int sig = 1) {
    for (int i = 0; i < len; ++i) {
        if (i < rev[i]) {
            swap(a[i], a[rev[i]]); 
        }
    }
    for (int mid = 1; mid < len; mid <<= 1) {     
        // 从分到最后的最小长度枚举回整个序列的一半
        cp w(cos(pi / mid, sig * sin(pi / mid)));
        for (int i = 0; i < len; i += mid * 2) {
            // 枚举合成一块的长度，i是处理到当前块的下标
            cp w0(1, 0);
            for (int j = 0; j < mid; j++, w0 *= w) {
                // 只需要枚举前半部分，后一半序列自然会算出来。
                auto x = a[i + j], y = w0 * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}
int main() {
    int n1, n2;
    // 若给出的格式是先输入两数字长度然后从低到高位给出序列，就可以写循环了。
    cin >> n1 >> n2;
    for (int i = 0; i < n1; ++i) { cin >> a[i].x; }
    for (int i = 0; i < n2; ++i) { cin >> b[i].x; }
    int bits = 0, m = 0;
    while ((1 << bits) < n1 + n2) { bits++; }
    m = 1 << bits;
    for (int i = 0; i < m; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bits - 1));
    fft(a, m), fft(b, m);
    for (int i = 0; i < m; ++i) { a[i] = a[i] * b[i]; }
    // 请读者自行写好逐项除以序列长度、完成进位与前导零去除的代码。
    return 0;
}

到这里，FFT应该不难了。（大雾）

2023 拓展

对于连续型的傅里叶变换，即对 $\mathbb{R}$ 上的卷积有

$$C(x) =(A*B)(x)= \int_{-\infty}^{+\infty} A(t)B(x-t) \mathrm{d}t$$

对函数 $f(x)$ 施行傅里叶变换，定义为以下的形式。

$$\mathscr{F}[f(x)]=\int_{-\infty}^{+\infty}f(x)\exp(-j\omega x)\mathrm{d}x = F(\omega)$$

逆变换回来有

$$\begin{aligned} \mathscr{F}^{-1}[F(\omega)] &=\dfrac{1}{2\pi }\int_{-\infty}^{+\infty}F(\omega)\exp(j\omega x)\mathrm{d}\omega \\ &= \dfrac{1}{2\pi }\int_{-\infty}^{+\infty}(\int_{-\infty}^{+\infty}f(t)\exp(-j\omega t)\mathrm{d}t)\exp(j\omega x)\mathrm{d}\omega \\ &=\dfrac{1}{2\pi }\int_{-\infty}^{+\infty}(\int_{-\infty}^{+\infty}\exp[-j\omega (t-x)]\mathrm{d}\omega)f(t)\mathrm{d}t \\ &=\dfrac{1}{2\pi }\int_{-\infty}^{+\infty}2\pi \delta(t-x)f(t)\mathrm{d}t = f(x). \end{aligned}$$

其中 $\delta(t-x)$ 为狄拉克 $\delta$ 函数，相关含义不在此叙述。因而，对卷积后的式子有

$$\begin{aligned}\mathscr{F}[(f*g)(x)] &= \int_{-\infty}^{+\infty}\left(\int_{-\infty}^{+\infty}f(t)g(x-t)\mathrm{d}t\right)\exp(-j\omega x)\mathrm{d}x \\&=\int_{-\infty}^{+\infty}f(t)\left(\int_{-\infty}^{+\infty}g(x-t)\exp(-j\omega x)\mathrm{d}x\right)\mathrm{d}t \\&=\int_{-\infty}^{+\infty}f(t)\exp(-j\omega t)\mathrm{d}t\left(\int_{-\infty}^{+\infty}g(x-t)\exp[-j\omega (x-t)]\mathrm{d}(x-t)\right) \\ \\&=\mathscr{F}[f(t)]\cdot\mathscr{F}[g(x-t)] = F_1(\omega)F_2(\omega).\end{aligned}$$

然后就发现，加法卷积（乘法）要算的东西，可以是经变换后点值相乘再求逆的结果。只需要选取合适数量的单位根，能算完结果即可。

同时也能注意到，积分和求导操作实际上也是一个矩阵。

如 $f(x+T) = \exp(T\dfrac{\mathrm{d}}{\mathrm{d} x}) f(x)$。感兴趣的读者可以去找算子代数来读，此处不多介绍。