问题：给定两个大整数 $A$ 和 $B$， $A$ 和 $B$ 的长度为 $n$ 和 $m$，求 $A$ 和 $B$ 的乘积

1. 朴素做法

思考小学数学中两个数的乘法列竖式的方式，对于 $B$ 的每一位，分别乘上 $A$ 的每一位，得到的结果相加。

以 $23\times 456$ 为例：

$\begin{array}{rl}23\times 456& =6\times 3+6\times 20+50\times 3+50\times 20+400\times 3+400\times 20\\ & =18+120+150+1000+1200+8000\\ & =10488\end{array}$

如此需要进行 $2\times 3$ 次乘法以及 $2\times 3-1$ 次加法。

对于到长度为 $n$ 的数 $A$ 和长度为 $m$ 的数 $B$ ，则需要进行 $n\times m$ 次乘法和 $n\times m-1$ 次加法

如此时间复杂度为：$O(nm)$，不妨设 $m\le n$，则时间复杂度为 $O(n^2)$

2. 快速傅里叶变换

上述的朴素做法时间复杂度过高，当两者长度都达到 $10^5$ 级别时，就很难快速得到两个大整数的乘积了。

2.1 求解思想

2.1.1 大整数转换为多项式

考虑将两个大整数转换成多项式：

令 $A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$，$B(x)=b_0+b_{1}x+b_2{x}^2+\cdots +b_{m-1}{x}^{m-1}$

以上多项式的表示称为多项式的系数表示法

以 $23\times 456$ 为例：

当 $a_0=3,a_1=2$，$A(x)=3+2x$

当 $b_0=6,b_1=5,b_2=4,B(x)=6+5x+4x^2$

可以看出：

$23$ 可以转换为 $a_0=3,a_1=2$ 的多项式，$A(10)=23$

$456$ 可以转换为 $b_0=6,b_1=5,b_2=4$ 的多项式，当 $B(10)=456$

那么 $A(10)\times B(10)$ 即为 $23\times 456$ 的结果

问题转换成了求解两个多项式 $A(x)$ 与 $B(x)$ 的乘积

令 $C(x)=A(x)\times B(x)=c_0+c_{1}x+c_2{x}^2+\cdots +c_{n+m-2}{x}^{n+m-2}$

则 $c_k=\sum _{i+j=k}{a}_i\times b_j$

求 $c_k$ 的过程也称为卷积，整体的时间复杂度也为 $O(n^2)$

2.1.2 多项式的点值表示法

上述将大整数转换为多项式的系数表示法进行多项式的计算，时间复杂度仍未改进。

事实上，多项式还有另外一种表示法，即点值表示法。

我们思考一条直线 $y=kx+b$，如果转换为上述表达，即为：$A(x)=a_0+a_{1}x$，其中$a_0=b,a_1=k$。而两点确定一条直线，则我们如果找出两个不同点，计算出斜率 $k$ 和截距 $b$ ，则唯一确定一条直线，即唯一确定一个一阶多项式。

推广到更高阶的多项式，对于 $n$ 阶多项式 $A(x)=\sum _{i=0}^n{a}_i{x}^i$，如果找到 $n+1$ 个不同点，则可以唯一确定一个 $n$ 阶多项式。二维平面上 $n+1$ 个不同点可以唯一确定一个 $n$ 阶多项式)。

对于多项式 $A(x)$，其点值表示法为：$\{(x_0,A(x_0)),(x_1,(A(x_1)),\cdots ,(x_n,A(x_n))\}$

要证明 $n$ 阶多项式可以由 $n+1$ 个点唯一确定，只需要证明这 $n+1$ 个系数是唯一的

对于 $A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$

将 $n+1$ 个点带入多项式，得到 $n+1$ 个方程

$A(x_0)=a_0+a_1{x}_0+a_2{x}_0^2+\cdots +a_n{x}_0^n$

$A(x_1)=a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n$

$\cdots$

$A(x_n)=a_0+a_1{x}_n+a_2{x}_n^2+\cdots +a_n{x}_n^n$

写成矩阵形式：

$\left[\begin{array}{c}A\left(x_0\right)\\ A\left(x_1\right)\\ ⋮\\ A\left(x_d\right)\end{array}\right]=\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^n\\ 1& x_1& x_1^2& \cdots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \cdots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]$

令方阵为矩阵 $M$ ，$M$ 也称为范德蒙德矩阵， $x_0,x_1,\cdots ,x_n$ 两两不同。证明 $n+1$ 个系数唯一只需要证明 $M$ 的行列式不为 $0$ 即可，而范德蒙德行列式的值一定不为 $0$，故得证二维平面上 $n+1$ 个不同点可以唯一确定一个 $n$ 阶多项式

至此，对于对于 $n$ 阶多项式 $A(x)=\sum _{i=0}^n{a}_i{x}^i$，我们有两种多项式表示法：

• 系数表示法，给出 $n+1$ 个系数

​ $[a_0,a_1,\cdots ,a_n]$

• 点值表示法，给出 $n+1$ 个点

  $\{(x_0,A(x_0)),(x_1,(A(x_1)),\cdots ,(x_n,A(x_n))\}$

2.1.3 多项式乘积的计算流程

对于两个多项式，$A(x)=\sum _{i=0}^n{a}_i{x}^i$ 和 $B(x)=\sum _{i=0}^m{b}_i{x}^i$
两个多项式的乘积为 $C(x)=A(x)\cdot B(x)=\sum _{i=0}^{m+n}{c}_i{x}^i$ 共有$n+m+1$ 项

因此对于多项式 $A(x)$ 和 $B(x)$ 各取 $n+m+1$ 个点

对于$A(x)$，取点为：$\{(x_0,A(x_0)),(x_1,(A(x_1)),\cdots ,(x_{n+m},A(x_{n+m}))\}$

对于$B(x)$，取点为：$\{(x_0,B(x_0)),(x_1,(B(x_1)),\cdots ,(x_{n+m},B(x_{n+m}))\}$

得到的 $C(x)$ 的 $n+m+1$ 个点为：

$\{(x_0,A(x_0)\times B(x_0)),(x_1,(A(x_1)\times B(x_1)),\cdots ,(x_{n+m},A(x_{n+m}\times B(x_{n+m})))\}$

该乘积的时间复杂度为：$O(n)$

如何得到 $C(x)$ 的系数表示法呢？

• 将 $A(x)$ 和 $B(x)$ 从系数表示法转换成点值表示法，在 $A(x)$ 和 $B(x)$ 中各取 $n+m+1$ 个点， 这 $n+m+1$ 个点的横坐标对应相同
• $O(n)$ 时间计算出每个点的纵坐标乘积，即得到表示 $C(x)$ 的$n+m+1$ 个点
• 将 $C(x)$ 从点值表示法转换成系数表示法

2.2 求解流程

2.2.1 从系数表示到点值表示（FFT）

朴素做法

在 $A(x)$ 和 $B(x)$ 中各取 $n+m+1$ 个点，对于每个点，选择其横坐标后，通过系数表示法计算出其对应的纵坐标，需要的时间复杂度为 $O(n)$，总时间复杂度为 $O(n^2)$，因此时间复杂度不满足要求。

快速傅立叶变换（Fast Fourier Transform）是一种快速进行离散傅立叶变换的方法, 通过单位根的性质, 将 $n$ 次求值的时间复杂度降为 $O(n\log n)$

$\mathrm{FFT}$ 的 $n$ 保证为 $2$ 的幂

$A(x)\to a={\left[a_0,a_1,a_2,\cdots ,a_{n-1}\right]}^{\top }$ , 将其分为偶数项和奇数项两个向量

$$\begin{gathered} 偶a^{[0]}=\left[a_0,a_2,\cdots ,a_{n-2}\right] \\ 奇a^{[1]}=\left[a_1,a_3,\cdots ,a_{n-1}\right] \end{gathered}$$

$\begin{array}{l}A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}^{x}n-1\\ A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a_{n-2}{x}^{\frac{n}{2}-1}\\ A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a_{n-1}{x}^{\frac{n}{2}-1}\end{array}$

$\begin{array}{l}{A}^{[0]}\left(x^2\right)=a_0+a_2{x}^2+a_4{x}^4+\cdots +a_{n-2}{x}^{n-2}\\ A^{[1]}\left(x^2\right)=a_1+a_3{x}^2+a_5{x}^4+\cdots +a_{n-1}{x}^{n-2}\\ x{A}^{[1]}\left(x^2\right)=a_{1}x+a_3{x}^3+a_5{x}^5+\cdots +a_{n-1}{x}^{n-1}\end{array}$

那么有: $A(x)=A^{[0]}\left(x{2}\right)+x A^{[1]}\left(x{2}\right) $

原问题: 求 $A(x)$ 在每个单位根上的值
转化为: 求次数界为 $\frac{n}{2}$ 的两个多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 在每个单位根平方上的值, 再合并。

$\begin{array}{l}A\left(w_n^k\right)=A^{[0]}\left({\left(w_n^k\right)}^2\right)+w_n^k{A}^{[1]}\left({\left(w_n^k\right)}^2\right)\\ A\left({\omega }_n^{k+\frac{n}{2}}\right)=A^{[0]}\left({\left({\omega }_n^{k+\frac{n}{2}}\right)}^2\right)+{\omega }_n^{k+\frac{n}{2}}{A}^{[1]}\left({\left({\omega }_n^{k+\frac{n}{2}}\right)}^2\right)\end{array}$

通过消去引理和折半引理, 将上式化简为:

$\begin{array}{ll}& A\left({\omega }_n^k\right)=A^{[0]}\left({\omega }_{\frac{n}{2}}^k\right)+{\omega }_n^k{A}^{[1]}\left({\omega }_{\frac{n}{2}}^k\right)\\ & A\left({\omega }_n^{k+\frac{n}{2}}\right)=A^{[0]}\left({\omega }_{\frac{n}{2}}^k\right)-{\omega }_n^k{A}^{[1]}\left({\omega }_{\frac{n}{2}}^k\right)\end{array}$

$A^{[0]}(x)$ 和 $A^{[1]}(x)$ 恰好均为次数界为 $\frac{n}{2}$ 的多项式

${\mathrm{DFT}}_n(a)=\{\begin{array}{l}{\mathrm{DFT}}_{\frac{n}{2}}\left(a^{[0]}\right)\\ {\mathrm{DFT}}_{\frac{n}{2}}\left(a^{[1]}\right)\end{array}$

所以问题又转化为了求 次数界为 $\frac{n}{2}$ 的多项式 $A^{[0]}(x)$和 $A^{[1]}(x)$
这两个问题和原问题的描述类似, 原问题的次数界为 $n$ , 子问题的次数界为 $\frac{n}{2}$ , 所以这两个问题是范围缩小一 半的两个子问题, 通过同样的方式求解这两个子问题, 那么就是熟悉的递归求解了。再根据折半引理, 就可以快速 地合并结果。

$T(n)=2T\left(\frac{n}{2}\right)+O(n)$

根据主定理, 时间复杂度为: $O(n\log n)$

2.2.2 点值表示的乘积

对于$A(x)$，取点为：$\{(x_0,A(x_0)),(x_1,(A(x_1)),\cdots ,(x_{n+m},A(x_{n+m}))\}$

对于$B(x)$，取点为：$\{(x_0,B(x_0)),(x_1,(B(x_1)),\cdots ,(x_{n+m},B(x_{n+m}))\}$

得到的 $C(x)$ 的 $n+m+1$ 个点为：

$\{(x_0,A(x_0)\times B(x_0)),(x_1,(A(x_1)\times B(x_1)),\cdots ,(x_{n+m},A(x_{n+m}\times B(x_{n+m})))\}$

该乘积的时间复杂度为：$O(n)$

2.2.3 从点值表示到系数表示（IFFT）

朴素做法

已知 $n+1$ 个点组成的多项式点值表示 { x 0 , A ( x 0 ) , ( x 1 , A ( x 1 ) ) , … , ( x n , A ( x n ) ) } \{x_0,A(x_0),(x_1,A(x_1)),\ldots,(x_n,A(x_n))\} {x0​,A(x0​),(x1​,A(x1​)),…,(xn​,A(xn​))}，求 $A(x)$ 的过程叫作插值

拉格朗日插值公式计算：$A(x)=\sum _{i=0}^{n}A(x_i)\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}$

时间复杂度：$O(n^2)$

顾名思义，快速傅里叶逆变换（Inverse Fast Fourier Transform）是快速傅里叶变换的逆。

$\mathrm{DFT}\to y_i=\sum _{j=0}^{n-1}{w}_n^{ij}{a}_j\to y=V_{n}a$

$V_n=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]$

这里的 $V_n$ 是一个范德蒙德矩阵, 它的第 $i$ 行第 $j$ 列尾主 $n$ 次单位根的 $ij$ 次方

${\left(V_n\right)}_{ij}={\omega }_n^{ij}$

离散傅立叶变换给出 $a$ 求 $y$ , 而离散傅立叶逆变换给出 $y$ 求 $a$ , 所以要用 $V_n$ 的逆矩阵乘向量 $y$ , 即 $a=V_n^{-1}y$ , 其中

${\left(V_n^{-1}\right)}_{ij}=\frac{{\omega }_n^{-ij}}{n}$

${\left(V_n^{-1}{V}_n\right)}_{ij}={\sum }_{k=0}^{n-1}\frac{{\omega }_n^{-ki}}{n}\times {\omega }_n^{kj}={\sum }_{k=0}^{n-1}\frac{{\omega }_n^{k(j-i)}}{n}$

当且仅当 $i=j$ 时, ${\left(V_n^{-1}{V}_n\right)}_{ij}=1$ , 其余情况 ${\left(V_n^{-1}{V}_n\right)}_{ij}=0$ , 所以结果是一个单位矩阵, 即 $V_n^{-1}{V}_n=I_n$ , 所以 ${\left(V_n^{-1}\right)}_{ij}=\frac{{\omega }_n^{-ij}}{n}$ 得证

可以发现

$\begin{array}{l}\mathrm{IDFT}\to a_i=V_n^{-1}y={\sum }_{j=0}^{n-1}\frac{{\omega }_n^{-ij}}{n}=\frac{1}{n}{\sum }_{j=0}^{n-1}{\omega }_n^{-ij}{y}_j\\ \mathrm{DFT}\to y_i={\sum }_{j=0}^{n-1}{\omega }_n^{ij}{a}_j\end{array}$

将 $\mathrm{FFT}$ 中的 ${\omega }_n=e^{\frac{2\pi }{n}i}=\cos \left(\frac{2\pi }{n}\right)+i\sin \left(\frac{2\pi }{n}\right)$
改成 ${\omega }_n=e^{-\frac{2\pi }{n}i}=\cos \left(\frac{2\pi }{n}\right)-i\sin \left(\frac{2\pi }{n}\right)$
如此, 对 $\mathrm{FFT}$ 做部分修改即可由点值表示转换为系数表示，由于仅参数不同，故时间复杂度仍为 $O(n\log n)$。

3. 快速数论变换

之前介绍了FFT，一种利用复数性质加速多项式乘积运算的算法

可以发现，在计算时，实部和虚部因为需要使用 $cos$ 和 $sin$ 以及 $\pi$，必须存储为浮点形式，复数的乘法运算也较为复杂。

计算机中的浮点数不免存在精度误差问题，所以我们考虑是否存在一种整数域内的形式，将复数转换为实数？

3.1 部分概念

3.1.1 欧拉函数

欧拉函数 $\phi (n)$ 表示小于等于 $n$ 的数中与 $n$ 互质的数的个数

3.1.2 欧拉定理

当 $\gcd (a,n)=1$ ，存在 $a^{\phi (n)}\equiv 1(\mathrm{mod}n)$

3.1.3 费马小定理

当 $p$ 是质数，欧拉定理为 $a^{p-1}\equiv 1(\mathrm{mod}p)$

费马小定理就是欧拉定理的特殊形式

3.1.4 阶

由欧拉定理可知, 对 $a\in \mathbb{Z},m\in {\mathbb{N}}^{\ast }$ , 若 $\gcd (a,m)=1$ , 则 $a^{\phi (m)}\equiv 1(\mathrm{mod}m)$ 。因此满足同余式 $a^n\equiv 1(\mathrm{mod}m)$ 的最小正整数 $n$ 存在, 这个 $n$ 称作 $a$ 模 $m$ 的阶, 记作 ${\delta }_m(a)$

3.1.5 原根

当 $a$ 模 $m$ 的阶为 $\phi (m)$ ，即 ${\delta }_m(a)=\phi (m)$ ，则称 $a$ 为模 $m$ 的原根

3.2 用原根代替复数根

复数根：${\omega }_n^k=(e^{\frac{2\pi }{n}i}{)}^k$

用原根来替代复数根的每个实数根：$g_n^k=(G^{\frac{p-1}{n}}{)}^k\mathrm{mod}p$

对原根而言：

$\begin{array}{l}{g}_n^n=g^{p-1}=g^{\phi (p)}=1\equiv (\mathrm{mod}p)\\ g_n^{n/2}=g^{(p-1)/2}\end{array}$

而:

${\left(g^{(p-1)/2}\right)}^2\equiv g^{p-1}\equiv 1(\mathrm{mod}p)$

即求方程 $x^2\equiv 1(\mathrm{mod}p)$ , 在 $p$ 为素数时只有 $\pm 1$ 两种取值。由于 $g$ 是原根, 只有 $g_n^n\equiv 1(\mathrm{mod}p)$ , 故 $g^{(p-1)/2}\equiv -1(\mathrm{mod}p)$

原根满足:

• 具有周期性
• 一个周期内的数均不相同
• $g_n^k=-g_n^{k+n/2}$

如此我们就找到了一个可以替换复数根的实数根，常用原根模数 $9982442353$，原根为 $3$，因为 $998244352=2^{23}\cdot 7\cdot 17$

因为 $2^{23}=8388608$，对于长度在 $8\times 10^6$ 内的大整数乘法足够使用。

4. 提交

从 LeetCode Multiply Strings 一题的提交结果来看：

• FFT 的递归版本比 FFT 的迭代版本需要更多的运行时间和空间。这是因为 FFT 的递归版本会使用额外的栈空间用来辅助计算，开辟和回收这些空间也需要一定的时间，因此相较于 FFT 的迭代版本，FFT 的递归版本的运行时间和空间消耗更多
• NTT 的迭代版本比 FFT 的迭代版本需要的运行时间和空间都更少。这是因为 NTT 使用的是原根，相较于单位复根需要的空间更少，计算次数也更少。

总体来说，FFT 和 NTT 的时间复杂度是相同的，均为 $O(n\log n)$。NTT 在一定大小范围内的大整数乘法运算中不论是空间还是时间都好过于 FFT，计算的常数低于 FFT ，精度也高于 FFT ，唯一的问题是需要找到合适的原根保证可以覆盖运算需要的范围，当寻找到的原根可以覆盖需要的运算范围，可以选择 NTT ，否则就选择 FFT。

5. 补充

5.1 复数

形如 $z=a+bi$ 的数称为复数，其中：$a,b\in R,i^2=-1,i$ 虚数单位

复数的幅角表示形式：$z=r{e}^{i\theta }$，$r$ 为模长，$\theta$ 为幅角

加法法则：$(a+bi)+(c+di)=a+c+(b+d)i$

减法法则：$(a+bi)-(c+di)=a-c+(b-d)i$

乘法法则：$(a+bi)\times (c+di)=ac+adi+bci+bd{i}^2=ac-bd+(ad+bc)i$

除法法则：$(a+bi)÷(c+di)=\frac{a+bi}{c+di}=\frac{ac+bd+(bc-ad)i}{c^2-d^2}$

欧拉定理：$e^{i\theta }=\cos \theta +i\sin \theta$

任意一个复数 $a+bi$ 对应复平面上的一个点 $(a,b)$ ，横坐标为实部，纵坐标为虚部。

模长 $r$ 即复平面上对应的点到原点的距离，即 $\sqrt{a^2+b^2}$

幅角 $\theta$ 即复平面上对应的点到原点的连线与 $x$ 轴的夹角

5.2 单位根

$z=a+bi=r{e}^{i\theta }=r(\cos \theta +i\sin \theta )$

两个复数相乘的结果为：两个复数的模相乘，再乘以 $e^{幅角的和}$

单位根是方程 $z^n=1$ 在复数范围内的 $n$ 个根，因为是单位根，因此模长 $r$ 为 $1$，本质就是求解 $n$ 个不同的 $\theta$，使得每个 $\theta$ 对应的 $z$ 都有 $z^n=1$ 。

$$\begin{gathered} z^n=r^n{e}^{n\theta i}=1 \\ \Rightarrow \{\begin{array}{l}r=1\\ n\theta =2k\pi \end{array} \end{gathered}$$

所以 ${\omega }^n=1$ 的 $n$ 个根为：${\omega }_n^k=e^{i\frac{2k\pi }{n}}=\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n},(k=0,1,2,...,n-1)$

主 $n$ 次单位根 ${\omega }_n=e^{i\frac{2\pi }{n}}$

5.2.1 消去引理

${\omega }_{dn}^{dk}={\omega }_n^k$

证明：${\omega }_{dn}^{dk}=(e^{\frac{2\pi }{dn}i}{)}^{dk}=(e^{\frac{2\pi }{n}i}{)}^k={\omega }_n^k$

5.2.2 折半引理

$(w_n^{k+\frac{n}{2}}{)}^2=({\omega }_n^k{)}^2={\omega }_{\frac{n}{2}}^k$

证明：${\omega }_n^{k+\frac{n}{2}}={\omega }_n^k{\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$ ，其中 ${\omega }_n^{\frac{n}{2}}={\omega }_2=-1$ ，$({\omega }_n^k{)}^2={\omega }_n^{2k}={\omega }_{\frac{n}{2}}^k$

5.2.3 求和引理

$\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=0$

证明：$\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{1^k-1}{{\omega }_n^k-1}=0$

5.3 参考视频和博客

1. 课程PPT

2. 快速傅里叶变换(FFT)——有史以来最巧妙的算法？

3. 算法学习笔记(42): 快速数论变换

4. 快速数论变换（NTT）超详解