利用Karatsuba算法、Toom-Cook算法和傅里叶变换进行快速乘法的奇特方法

介绍

算法的适用性和性能取决于某些常规计算的执行速度。 例如，在椭圆曲线密码学中，需要计算公钥为 $k \times g = g + g + g + \dots + g$，其中 $k$ 是一个非常大的整数（通常是一个数百位的数字），而 $g$ 是椭圆曲线上的一个点 $(x, y)$，被称为生成元。 如果以一种简单的方式完成，即重复将 $g$ 加到自身，则需要大约 $10^{100}$ 次运算（我们说算法以 $O(n)$ 运行，表明运算次数 - 直到某个常数因子 - 随着 $k$ 以线性方式增加）。 最快的超级计算机每秒可以执行少于 $10^{18}$ 次浮点运算； 光做一次计算就需要永远的时间。 尽管如此，由于开发了更快的算法，我们每天都在执行此类计算或相关计算（例如，我们可以通过重复添加 $g + g = 2g$，$2g + 2g = 4g$ 等来减少计算时间，从而将计算次数减少到 $O(\log(n))$ - 比较 $10^{100}$ 到类似 $100 \times \log_2 10$ 的东西）。

zk-SNARKs（零知识简洁非交互式知识论证）是重要的密码学原语，允许一方（证明者）说服另一方（验证者）某个陈述是真实的，而除了该陈述的有效性之外不透露任何其他内容。 zk-SNARKs 的应用范围非常广泛，因为它们有潜力成为新型治理、数据共享和金融系统的基础。 例如，你可以将一项困难的计算委托给不受信任的第三方，并获得一个证明，使你可以验证计算的完整性，而无需重新运行所有内容。 关键是证明是简洁的，因此可以在数百毫秒内验证，而不是执行整个计算。 构造依赖于将计算转换为多项式并检查多项式上的条件。 多项式乘法可以通过快速傅里叶变换以非常有效的方式完成 - 这是人类有史以来设计的最重要的算法之一。 此外，这种计算可以并行化：多个处理器可以运行算法的各个部分，使其更快。

即使是像整数乘法这样的简单计算（几乎发生在所有地方）也可以比学校规则更快地完成，只要我们尝试相乘的数字足够大。

如果你想学习如何加速一些普通计算并使你的算法运行得更快，那么接下来的章节适合你。

分而治之：Karatsuba

我们都在小学学过如何将两个数字相乘：我们将一个写在另一个下面，然后将上面的每个数字乘以下面的数字的每一位，然后我们将所有数字加起来：

          1234
×          152
———————————————
      2468 ( =  1234 ×     2)
     6170  ( =  1234 ×    50)
    1234   ( =  1234 ×   100)
———————————————
    187568 ( = 187568)

此算法具有 $O(n^2)$。 1960 年， Kolmogorov 推测这代表了乘法的渐近界限（即，两个数字的乘法不能少于 $O(n^2)$ 次运算）。 他就此主题发表演讲，其中一位学生 Karatsuba 当时 23 岁，提出了一个以 $O(n^{\log_2(3)})$ 运行的解决方案，从而推翻了 Kolmogorov 的猜想。 Karatsuba 算法的基本思想如下：假设我们想将 $x$ 和 $y$ 相乘； 我们可以将它们分解为更小的数字：

$x = x_1 \times 10^m + x_0$

$y = y_1 \times 10^m + y_0$

其中 $x_0$ 和 $y_0$ 都是小于 $10^m$ 的数字。 产品 $x \times y$ 仅仅是：

$x \times y = x_1 \times y_1 \times 10^{2m} + (x_1 \times y_0 + y_1 \times x_0) \times 10^m + x_0 y_0$

Karatsuba 发现可以有效地计算 $x_1 y_0 + y_1 x_0$，但需要进行一些加法运算：

$x_1 \times y_0 + y_1 \times x_0 = (x_1 + x_0) \times (y_1 + y_0) - x_1 \times y_1 - x_0 \times y_0$。

即使有一些额外的计算，这些计算也是在较小的数字上进行的，从而导致大数字的整体成本更低。

Toom-Cook 算法

可以进一步采用分而治之的策略，从而降低乘法算法的复杂性。 Toom 和 Cook 开发了几种方法（称为 Toom-X，其中 X 是一个数字），包括以下阶段：

1. 分割
2. 评估
3. 逐点乘法
4. 插值
5. 重组

该算法的几种变体在 GNU 多精度算术库 中实现。 Toom-2 与 Karatsuba 算法相同。 Toom-X 首先将数字 $x$ 和 $y$ 分成 X 个长度相等的部分(1)，这些部分被视为某些多项式的系数（我们专注于 Toom-3，但你可以在此处 中查看更多详细信息）(2)：

$x(t) = x_2 t^2 + x_1 t + x_0$

$y(t) = y_2 t^2 + y_1 t + y_0$

如果我们在 $t = b$ 处评估 $x$，$y$，我们会得到数字。 两个数字的乘积等于一个 $2(X-1)$ 次多项式，

$w(t) = w_4 t^4 + w_3 t^3 + w_2 t^2 + w_1 t + w_0$ 我们可以在 5 个不同的点评估多项式，这足以根据插值定理唯一地确定多项式 $w$。 我们可以选择 5 个方便的点，这使得多项式的评估和重构变得容易。 常见的点有 $0, 1, -1, 2$ 和 $\infty$（最后一个只是主要系数的乘积）。 让我们看看每个值的形式：

$w(0) = x(0) y(0) = x_0 y_0$

$w(1) = x(1) y(1) = (x_0 + x_1 + x_2)(y_0 + y_1 + y_2)$

$w(-1) = x(-1) y(-1) = (x_0 - x_1 + x_2)(y_0 - y_1 + y_2)$

$w(2) = x(2) y(2) = (x_0 + 2x_1 + 4x_2)(y_0 + 2y_1 + 4y_2)$

$w(\infty) = x(\infty) y(\infty) = x_2 y_2$

如果我们从 $w$ 及其系数来看待事物，我们会得到：

$w(0) = w_0$

$w(1) = w_4 + w_3 + w_2 + w_1 + w_0$

$w(-1) = w_4 - w_3 + w_2 - w_1 + w_0$

$w(2) = 16w_4 + 8w_3 + 4w_2 + 2w_1 + w_0$

$w(\infty) = w_4$

这只是求解一个线性系统（其中 2 个系数很简单）。 一旦系数已知，剩下的就是评估 $t = b$ 处的 $w$ 并相加。 Toom-3 比 Karatsuba 方法（$\mathcal{O}(n^{1.58})$）具有较低的阶数（$\mathcal{O}(n^{\log(5)/\log(3)}) = \mathcal{O}(n^{1.46})$），因此对于足够大的整数，它的运行速度更快。

对于更大的整数（大约 10,000 到 40,000 位），我们可以通过 Schönhage-Strassen 算法更快地运行，该算法使用快速傅里叶变换 (FFT) 来实现复杂度 $O(n \log(n) \log(\log(n)))$。 在我们可以解释该算法之前，我们需要先介绍 FFT。 阶数可以进一步降低到 $O(n \log(n))$，但该算法仅适用于（超超）非常大的数字，并且是银河算法的一个例子。

快速傅里叶变换

FFT 是许多重要算法的关键构建块之一，例如非常大的数字的快速乘法、多项式乘法、求解有限差分方程、纠错码（Reed-Solomon 码）和数字信号处理。 当 Gauss 试图插值小行星 Pallas 和 Juno 的轨道时，他在 19 世纪初使用了它。 一个简单的实现需要 $O(n^2)$ 操作。 1965 年，Cooley 和 Tukey 意识到该算法可以更有效地实现，将其降低到 $O(n \log(n))$，从而导致其广泛使用。 几乎每种语言和数值计算库都实现了它。 在 Rust 中，你可以查看此链接。

为了了解对朴素算法的巨大改进，让我们看一下不同样本的计算次数：

我们看到，对于具有 1000 个或更多元素的样本，计算量减少了两个数量级以上！

复数上的 FFT

傅里叶变换将函数从其原始域（空间或时间）映射到另一个取决于（空间或时间）频率的函数。 换句话说，它将一个函数分解为一组具有不同频率和幅度的正弦波，这些正弦波可用于分析给定系统的行为。 我们还可以执行反转，添加所有这些波以恢复原始函数。 即使（连续）傅里叶变换有许多应用，我们也将对离散傅里叶变换 (DFT) 感兴趣，其中我们有一个有限的数据集合。 给定数据 $x_0$、$x1$、...$x{N-1}$，DFT 给出一个序列 $X_0, X1, ... X{N-1}$，其中

$X = \sum_{k=0}^{N-1} x_k \exp(-2 \pi i k / N)$

其中 $i^2 = -1$ 是虚数单位。 DFT 的反演由下式给出

$x = \frac{1}{N} \sum_{k=0}^{N-1} X_k \exp(2 \pi i k / N)$。

DFT 可以表示为矩阵向量乘积的形式，$X = Mx$，其中 $M$ 是 $N \times N$ DFT 矩阵：

$M_{ij} = \omega^{(i-1) \times (j-1)}$

其中 $\omega = \exp(2 \pi i / N)$ 并且 $i$ 和 $j$ 取值 $1, 2, 3, ... N$

以这种方式实现，DFT 需要 $N^2$ 次运算，这是由向量矩阵乘法产生的。 FFT 将通过利用结构并使用分而治之的策略使此计算更有效。

我们也可以将 DFT 视为在单位根上评估系数为 $x_k$ 的多项式。 这在讨论快速多项式乘法时会很有用。

关键是计算具有 $N$ 个点的 DFT 可以减少到计算具有 $N/2$ 个点的两个 DFT。 我们可以递归地应用它，将一个非常大的问题分解为一组更小的、更容易解决的子问题，然后将这些结果重新组合以获得 DFT。

该算法还利用了复平面中 $n$ 次单位根的性质。 如果 $z^n = 1$，则数字 $z$ 称为 $n$ 次单位根。 这些是形式

$z_k = \exp(2 \pi i k / n)$ 对于 $k = 0, 1, 2, ..., n-1$。 一个有趣的点是这些根以共轭对出现：对于每个根 $r$，我们都有相应的 $\bar{r}$（事实上，它们在乘法下形成一个 $n$ 阶的有限群）。 例如，四次单位根是：$1, i, -1, -i$。 很容易看出哪些是对。

要了解所有工作原理，假设我们有一个向量 $x = (x_0, x_1, x2, ... x{n-1})$，并且我们想要计算 FFT。 我们可以将偶数项和奇数项分开：

$X = \sum{k=0}^{n/2-1} x{2k} \exp(2 \pi i 2k / n) + \sum{k=0}^{n/2-1} x{2k+1} \exp(2 \pi i (2k+1) / n)$

我们可以通过提出一个 $\exp(2 \pi i / n)$ 的因子，以不同的方式表达奇数项，

$X = \sum{k=0}^{n/2-1} x{2k} \exp(2 \pi i 2k / n) + \exp(2 \pi i / n) \sum{k=0}^{n/2-1} x{2k+1} \exp(2 \pi i (2k) / n)$

我们现在可以看到，只要 $k$ 大于 $n/2$，对应于 $n$ 次单位根的因子就会重复自身。 另一种看待这一点的方式是通过从指数的分子中取出 2 并将其发送到分母来重新排列这些项：

$X = \sum{k=0}^{n/2-1} x{2k} \exp(2 \pi i k / (n/2)) + \exp(2 \pi i / n) \sum{k=0}^{n/2-1} x{2k+1} \exp(2 \pi i (k) / (n/2))$

我们现在发现 $\sum{k=0}^{n/2-1} x{2k} \exp(2 \pi i k / (n/2)) = DFT(x{2k})$ 只是包含 $n/2$ 个点的偶数项的 DFT。 类似地，$\sum{k=0}^{n/2-1} x_{2k+1} \exp(2 \pi i (k) / (n/2))$ 是包含 $n/2$ 个点的奇数项的 DFT。 这样，我们将 $n$ 点 DFT 分解为两个较小的 $n/2$ 点 DFT，它们可以组合起来产生原始 DFT。 现在，每个 $n/2$ DFT 都可以分解为两个较小的 DFT，因此我们可以通过处理较小的样本来递归地减少计算量（这样，我们就可以省去大量的向量矩阵乘积）。

将 FFT 扩展到任意环

FFT 可以从复数或实数扩展到任意环，例如整数或多项式（查看我们的开发者数学生存工具包）。 特别是，我们可以使用数论变换，它将 FFT 专门用于 $\mathbb{Z}/p\mathbb{Z}$，即模 $p$ 的整数（一个素数）。 在这里，我们也有 $n$ 次单位根，由下式给出

$\alpha^n \equiv 1 \pmod{p}$

重要的是，我们将自己限制为素数：在这种情况下，我们得到 1 的平方根只是 1 和 $-1$。 例如，如果我们取 $p = 5$，$1^2 \equiv 1 \pmod{5}$ 和 $-1 \equiv 4$，$4^2 = 16 \equiv 1 \pmod{5}$。 这对于 8 不成立，因为 $1^2 \equiv 3^2 \equiv 5^2 \equiv 7^2 \equiv 1 \pmod{8}$ 并且我们将有 4 个平方根！

在有限域中使用 FFT 的问题在于，我们不能随意选择域和域。 我们需要选择一个阶数为 $2^n$ 的乘法子群（换句话说，我们需要选择一个由元素 $g$ 生成的群，并且该群包含其高达 $2^n$ 的幂）。 例如，如果我们取 $p = 5$，我们有一个阶数为 $4 = 2^2$ 的群，它由 2 生成：${2^1 = 2, 2^2 = 4, 2^3 \equiv 3, 2^4 \equiv 1}$； 尽管如此，它不需要跨越该域的所有元素。

FFT 乘法算法

该算法遵循与 Karatsuba 和 Toom 相同的路线：

1. 分割
2. 评估
3. 逐点乘法
4. 插值
5. 组合

关键区别在于使用 FFT 加速计算。

多项式乘法

让我们从多项式乘法开始。 给定两个多项式，$p(x) = pd x^d + p{d-1} x^{d-1} + \dots + p_0$ 和 $q(x) = qd x^d + q{d-1} x^{d-1} + \dots + q_0$，我们想找到它们的乘积，$w(x) = p(x)q(x)$。 最简单的算法是重复应用分配律，执行乘法并重新排列所有内容。 两个 $d$ 次多项式的乘积是一个 $2d$ 次多项式。 我们可以看到，这种策略涉及 $\mathcal{O}(d^2)$ 阶的运算，即运算随着所涉及多项式的次数呈二次方增长。 我们可以利用多项式的结构和插值定理。 我们至少有两种形式来描述同一个多项式：

• 给出 $d+1$ 个系数。
• 在精确的 $d+1$ 个点上指定多项式的值(3)。

第二种选择有什么优势？ 我们可以自由选择点并减少计算次数。 例如，如果我们有一个偶函数 $f(x) = f(-x)$，我们可以评估较少的点。 同样，如果函数是奇函数，$f(-x) = -f(x)$ 并且我们必须更改符号才能获得 $-x$ 的值。 因此，选择对 $x$ 和 $-x$，我们将评估次数减少一半（例如，除非我们选择 0）。 我们可以将多项式分成两个多项式：一个有奇数项，另一个有偶数项：

$p(x) = p_e(x) + x p_o(x)$。

例如，如果 $p = x^5 + 3x^4 + 5x^3 + 2x^2 + 6x + 3$，我们将其拆分：$p(x) = (3x^4 + 2x^2 + 3) + x(x^4 + 5x^2 + 6)$ 然后我们有：$p_e = (3x^4 + 2x^2 + 3)$ 和 $p_o = (x^4 + 5x^2 + 6)$，其中两个多项式都是偶函数！ 这样，我们可以很容易地看到：$p(-x) = p_e(x) - x p_o(x)$ 如果我们有对 $(x_k, p(x_k))$ 和 $(x_k, q(x_k))$，则在 $x_k$ 处评估的乘积多项式是 $(x_k, p(x_k) q(x_k))$。

为了确定乘积多项式，我们需要 $2d+1$ 个点; 利用上述策略，我们需要更少的点评估。 如果我们可以轻松地从系数形式转换为点评估，以该形式执行乘法，然后转换回系数形式，我们可以实现较低的复杂性。 我们可以将多项式 $p_e(x^2)$ 和 $p_o(x^2)$ 递归地分解为更小的多项式。

我们可以选择 $n$ 次单位根作为评估点，它们成对出现：$\exp(2 \pi i k / n)$，其中 $k = 0, 1, 2 ... n-1$。 换句话说，我们可以快速计算多项式的 DFT，将系数相乘并在找到乘积后反转 DFT。 这导致 $\mathcal{O}(d \log(d))$ 阶的运算。

整数乘法

要将 FFT 应用于整数乘法，我们需要将数字转换为多项式的系数，执行 FFT 乘法并最终重建结果。 总的来说，这将花费 $O(n \log(n) \log(\log(n)))$。 存在大量的开销，这将使该算法仅适用于非常大的整数。 例如，如果我们想将 3578 和 2457 相乘，我们可以定义向量 $(8, 7, 5, 3, 0, 0, 0, 0)$ 和 $(7, 5, 4, 2, 0, 0, 0, 0)$，其中我们方便地用零填充数字。

通常，运算是模 $2N+1$ 执行的，其中 $N$ 大于整数 $x$ 和 $y$ 的组合位数，以确保结果永远不会环绕。

傅里叶变换的优点是，诸如 $x$ 和 $y$ 的Rollup之类的运算可以从变换 $X$ 和 $Y$ 的乘积计算得出并变换回来：

$\sum_{k=0}^{N} xk y{N-k} = IFFT(FFT(y) \times FFT(x))$

Schönhage-Strassen 算法利用了负循环Rollup。 给定长度为 $n$ 的向量 $x$ 和 $y$，以及 $r$ 是一个 $2n$ 次（本原）单位根（即，$r^{2n} \equiv 1 \pmod{p}$ 且如果 $0 < k < 2n$，则 $r^k \not\equiv 1$），我们可以定义以下权向量：

对于 $0 \leq j < n$，$W_j = r^j$

对于 $0 \leq j < n$，$W_j^{-1} = r^{-j}$ $x$ 和 $y$ 的负循环Rollup (NCC) 可以计算为： $NCC(x, y) = W^{-1} IFFT(FFT(Wx) \times FFT(Wy))$

GNU 多精度算术库中实现的各种方法的比较显示在此链接中。

总结

选择正确的算法来执行常规计算（例如整数或多项式乘法）可能会对软件的性能产生巨大影响。 根据整数的大小，可以通过采用分而治之的方法来加快速度（减少计算次数）：我们将计算分解为较小的计算，这些计算很容易处理，或者继续分解它们直到它们可管理。 我们提出的所有快速算法都使用了这种方法，从而大大节省了计算量。 FFT 凭借其复杂度 $O(n \log(n))$ 可以成为加速计算的宝贵工具，即使它乍一看可能显得很奇怪或牵强！

笔记

(1) 如果这不可能，则最高有效部分可以比其余部分短。

(2) 为了方便起见，我们将删除乘法符号。

(3) 我们之前在 Toom-Cook 方法中提到过这一点。 例如，我们从几何学中知道我们需要给出两点来确定一条直线，这是一次多项式。

• 原文链接： blog.lambdaclass.com/wei...
• 登链社区 AI 助手，为大家转译优秀英文文章，如有翻译不通的地方，还请包涵～