整数乘法的长征

可以点开 pdf.

问题和计算模型

整数乘法无疑是再熟悉不过的问题了, 但是出于严谨性, 我们先在本文开头明确一下问题以及计算模型.

输入为两个非负整数 \(a, b\), 它们都是 \(n\) 位的二进制数. 我们的目标是计算它们的乘积 \(a \cdot b\), 也是按顺序输出它的各位数码. 原则上说, 也可以考虑固定某个常数 \(B\) 的 \(B\) 进制数, 将每一位的二进制编码输入, 这会给出不同的问题, 因为进制转换也需要一定时间. 但对于接下来给出的算法而言, 由于基本上没有涉及到二进制的特殊性质, 它们都可以直接推广到 \(B\) 进制数的情形, 给出相同的复杂度. 总而言之, 我们下文的讨论中都当做二进制数来处理.

选取的计算模型是衡量的一个一个位的计算时间, 也就是说像 word RAM 那样一次性可以操作一个字长的计算机是不被允许的. 下面所给出的算法基本上都可以应用于多带图灵机或者布尔电路. 也就是说, 接下来我们说一个算法的复杂度是 \(O(T(n))\) 的时候, 应该都能给出一个多带图灵机上 \(O(T(n))\) 时间的算法, 或者一个大小为 \(O(T(n))\) 的布尔电路来计算整数乘法.

Karatsuba 型算法

在人类漫长的历史中, 乘法毫无疑问是最基本的运算之一. 在上个世纪, 计算复杂度开始成为一个
严肃的研究问题之后, 大家开始审视算法所消耗的时间. 人们一直以来使用的朴素的竖式乘法
在衡量逐位进行的计算时, 时间复杂度是 \(O(n^2)\) 的, 其中 \(n\) 是乘数的位数.

于是, Kolmogorov 大胆地提出了如下猜想: ''会不会整数乘法的时间复杂度下界就是 \(O(n^2)\) 呢?'' 有趣的是, 没过多久, 他讨论班上的一个学生 Karatsuba 就证明了这个猜想是完全错误的.

Karatsuba 算法.
整数乘法可以在 \(O(n^{\log_2 3})\) 时间内完成. 注意 \(\log_2 3 < 1.585\).

这个 Karatsuba 算法的大概思路可以按照如下思路概括.

不妨设 \(n\) 是偶数, 将两个正整数写成 \(a = a_0 + a_1 \cdot 2^{n/2}\), \(b = b_0 + b_1 \cdot 2^{n/2}\) 的形式. 那么欲要求出 \(ab\), 也就是求出

\[ab = a_0b_0 + (a_0b_1 + a_1b_0)2^{n/2} + a_1b_1 2^n, \]

注意到

\[a_0b_1 + a_1b_0 = (a_0+a_1)(b_0+b_1) - a_0b_0 - a_1b_1, \]

而后面二者都是已经需要计算的, 所以我们将两个 \(n\) 维乘法问题转化为了三个 \(n/2\) 维乘法问题. 执行这个递归算法, 有递归式

\[T(n) = 3T(n/2) + O(n), \]

根据主定理解得 \(T(n) = O(n^{\log_2 3})\).

注意, 其实对于 \(a_0+a_1\) 和 \(b_0+b_1\), 我们实际上是 \(n+1\) 位的乘法, 但多出的位数可以先暴力计算出来, 这样上述递归式就是正确的了. 后续我们讨论更加复杂的算法的时候, 会因为技术细节稍微写出一个不那么正确的递归式, 但是大体上思路都会是对的.

上述算法用如下想法概括: 我们要计算多项式乘法 \((a_0 + a_1 X) (b_0 + b_1 X)\), 结果是一个 \(2\) 次多项式, 所以我们需要 \(3\) 个点值. 实际上, Karatsuba 所选取的点值就是 \(X = 0, 1, \infty\) 的三个点. 这里 \(\infty\) 点就对应于最高次项.

Toom--Cook 算法

显然, 按照多项式点值的理解思路, 我们可以继续将整数剖成更多的部分. 如果写成 \(a = a_0 + a_1 X + \cdots + a_k X^k\), 那么 \(ab\) 就是一个 \(2k\) 次多项式, 我们需要 \(2k+1\) 个点值. 对于 \(k=2\) 的情况, 根据递归式

\[T(n) = 5T(n/3) + O(n), \]

得到了 \(T(n) = O(n^{\log_3 5})\).

Toom--Cook 算法的一般形式是如下的递归式

\[T(n) = (2k+1)T(n/(k+1)) + k^{O(1)} n, \]

其中 \(k^{O(1)}\) 是进行 Lagrange 插值所需要支付的代价. 由于 \(\log_{k+1}(2k+1) \to 1\), 当 \(k\) 取任意大的整数, 可以得到复杂度 \(T(n) = O(n^{1+\varepsilon})\). 也可以进一步地分析 \(\varepsilon\) 的大小. 这个递归式的解是 \(T(n) = O(n^{\log_{k+1} (2k+1)} k^{O(1)})\). 写作

\[n^{\log_{k+1} (2k+1)} k^{O(1)} = n \cdot \exp \left\{ O\left(\frac {\log n}{\log k}\right) + O(\log k) \right\}, \]

取 \(\log k = \Theta \left(\sqrt{\log n}\right)\) 可以做到复杂度

\[n 2^{O\left(\sqrt{\log n}\right)}. \]

后续有一系列改进来优化这个 \(O(\cdot)\) 的大小.

在研究整数乘法的时候, 我们总是可以考虑类似的多项式乘法的问题. 如果只衡量四则运算的次数, 多项式乘法有一个好处, 就是它不会涉及进位,
很多问题的分析就会容易很多. 上面的 \(k^{O(1)}\) 如果考虑多项式乘法的时候, 注意如果使用暴力的 Lagrange 插值, 复杂度为 \(O(k^2)\)
插出多项式, 而我们将长为 \(n\) 的多项式分成了大小为 \(n/k\) 的块, 说明复杂度就是

\[T(n) = (2k+1)T(n/(k+1)) + O(kn). \]

这就可以看做一个标尺, 作为整数乘法用 Toom--Cook 算法的复杂度差不多能达到的最好的情况.

初步引入快速 Fourier 变换

现在一般认为, 早在 1876 年, C.F.Gauss 就已经发现了快速 Fourier 变换的算法.
而正式提出的快速 Fourier 变换被称作 Cooley--Tukey 算法.

快速 Fourier 变换, Cooley--Tukey 算法
令 \(n\) 是 \(2\) 的幂, 那么 \(n\) 个点的快速 Fourier 变换

\[b_k = \sum_j a_j \omega^{kj}, \]

给定序列 \(a\) 求 \(b\), 其中 \(\omega = e^{2\pi i / n}\), 可以在 \(O(n\log n)\) 次运算内完成.
其中运算是 \(\mathbb C\) 上的代数运算.

相信大家都已经很熟悉快速 Fourier 变换了, 这里不再赘述.

但是, 当我们严格讨论整数乘法的时候, 就有更多需要讨论的细节. 对于两个 \(n\) 位的整数乘法, 考虑拆解成多项式

\[a = \sum_{k=0}^{m-1} a_k X^k, \]

其中 \(m = n / \ell\), 那么每个 \(a_k\) 都是 \(\ell\) 位的整数. 这样多项式乘法

\[ \left( \sum_k c_k X^k \right) = \left(\sum_{k=0}^{m-1} a_k X^k\right) \left(\sum_{k=0}^{m-1} b_k X^k\right) \]

的各项系数不超过 \(n\cdot 4^\ell\), 如果取 \(\ell = O(\log n)\), 得数都是 \(O(\log n)\) 位的整数.

所以, 一种方法是将这些数都当做精度为 \(O(\log n)\) 位的复数 \(\mathbb C\) 来处理 (其实在 \(\mathbb C\) 中制备单位根以及精度分析都很需要进一步的论证, 不过暂时不谈), 另一种方式是将 \(O(\log n)\) 位的的数放在大小合适的有限域 \(\mathbb F_p\) 中处理. 二者都形如将问题规约成了 \(O(m\log m) = O(n)\) 个 \(O(\log n)\) 级别的复数或者有限域上的乘法. 这样这种算法的复杂度形如

\[T(n) = O(n) \cdot T(O(\log n)), \]

它的复杂度可以大概写作

\[\begin{align*} T(n) & \leq (K n) (K \log n) (K \log \log n) \cdots\\ &= O \left( K^{\log^* n} n (\log n) (\log \log n) \cdots \right), \end{align*} \]

其中这个 \(\cdots\) 是一直乘到 \(\log^{\circ k} n < 2\) 为止.

看起来这个 \(\log\) 迭代是不断规约子问题所难以避免的, 这也说明了接下来要介绍的 Schonhage--Strassen 算法的巧妙.

Schonhage--Strassen 算法

Schonhage--Strassen 算法.
整数乘法可以在 \(O(n\log n\log\log n)\) 时间内计算.

Schonhage 和 Strassen 的想法在直觉上大概可以解释如下: 快速 Fourier 变换的过程中, 涉及的运算只有某个环 \(R\) 上的元素的加减法, 以及元素乘以单位根的幂 \(\omega^i\). 然后, 整个卷积就通过

\[a * b = \mathcal F^{-1} \{ \mathcal F \{ a \} \cdot \mathcal F \{ b\} \} \]

完成. 也就是说, 真正的 \(R\) 中完成的任意两个元素的乘法只有 \(n\) 次, 而不是 \(n\log n\) 次, 如果能让环 \(R\) 的编码具有某些 更加美好的性质, 或许就能让其他部分的时间复杂度降低.

为了完成这一点, Schonhage--Strassen 算法的核心在于并不选取传统意义上单位根存在的域, 而是构造出合适的环.

具体而言, 他们首先考虑的是对于合适的 \(n\), 计算 \(\bmod (2^n + 1)\) 的乘法. 我们如果最后的目的是带入 \(T = 2^{n/m}\), 那么做拆解,

\[a = \sum_{i=0}^{m-1} a_i T^i, \]

就只需计算多项式乘法

\[\left(\sum_{i=0}^{m-1} a_i T^i\right)\left(\sum_{i=0}^{m-1} b_i T^i\right) \bmod (T^m + 1), \]

那么 \(T^m +1 = 0\) 需要 \(2m\) 次单位根. 我们注意系数 \(a_i\) 是 \(n/m\) 位的, 那么如果我们取 \(R = \mathbb Z / (2^{\ell}+1) \mathbb Z\), 其中 \(\ell \geq 2n / m + \log m\), 这样 \(R\) 中乘法的结果就是没有损失的.

最后, 选取的 \(\ell\) 满足了在 \(R\) 中, 令 \(\omega = 2\), 我们有 \(\omega^{2\ell} = 1\), 是一个 \(2\ell\) 次单位根, 所以这个环 \(R\) 最大可以支持 \(m=\ell\) 时候的多项式乘法. 这大概支持了 \(n\) 级别的问题规约到了 \(m \approx \sqrt{2n}\) 级别的问题.

现在我们发现, 在计算快速 Fourier 变换的时候, \(R\) 上的加减法都是线性的, 乘以 \(\omega^i\) 此时也变得十分简单, 就是移位之后进行取模. 最后, 逆 Fourier 变换需要除以 \(2^m\), 也是容易计算的. 唯一递归的部分是将两个经过了 Fourier 变换的数列进行点乘, 我们有递归式

\[T(n) = m T(\ell) + O(n\log n), \]

这里 \(m, l \approx \sqrt{2n}\), 我们不严谨地写成

\[T(n) = \sqrt{2n} T (\sqrt{2n}) + O(n\log n), \]

记 \(F(n) = T(2n)\), 有

\[F(n) = T(2n) = 2\sqrt n T(2\sqrt n) + O(n\log n) = 2\sqrt n F(\sqrt n) + O(n\log n), \]

这样就能看出, 第二层的复杂度是 \(2\sqrt n \cdot \sqrt n \log \sqrt n = n\log n\), 往下一直都是如此. 总共递归 \(\log \log n\) 层, 就得到了

\[T(n) = O(n\log n\log \log n). \]

Cantor--Kaltofen 算法

上面这个算法可以很轻易地改进成一个 \(\mathbb Q\) 上的多项式乘法的算法, 只进行 \(O(n\log n\log \log n)\) 次有理数的四则运算. 大概道理就是将 \(R\) 定为 \(\mathbb Q[T] / (T^\ell + 1)\).

如果能得到一个 \(\mathbb Z\) 上的算法就更好了, 这意味着彻底避免除法, 也就是说, 我们可以将这算法和任意环 \(A\) 做张量, 就能得到一个 \(A\) 上的多项式乘法算法. 这个 \(\mathbb Z\) 上的算法被称作 Cantor--Kaltofen 算法.

Cantor 和 Kaltofen 的思路其实也可以大致地概括. 他们注意到, 原来的 Schonhage--Strassen 算法直接改写到 \(\mathbb Q\) 上, 所涉及的除法, 所涉及的除法其实除以的都是 \(2\) 的整次幂.

类似地, 如果考虑 \(3^n\) 次本原单位根, 它满足的最小多项式是分圆多项式

\[\Phi_{3^n}(X) = (X^{3^n} - 1) / (X^{3^{n-1}} - 1) = X^{2\cdot 3^{n-1}} + X^{3^{n-1}} + 1, \]

也是稀疏的. 这意味着, 对于多项式乘法问题

\[a(X) \cdot b(X) \bmod \Phi_{3^n}(X), \]

有次数 \(N = 2\cdot 3^{n-1}\), 那么取环 \(R = \mathbb Z[T] / \Phi_{3^m}(T)\) 作为子问题, 有子问题规模 \(M = 2\cdot 3^{m-1}\), 我们分割为

\[a = \sum_{i=0}^{L - 1} a_i Y^{i}, \]

由于 \(a_i\) 次数是 \(3^{m}-1\), 最后替换为 \(Y = X^{3^{m-1}}\), 满足 \(\Phi_{3^{n-m+1}}(Y) = 0\), 说明计算的是 \(\bmod \Phi_{3^{n-m}}(Y)\) 的多项式乘法, 有 \(L = 2\cdot 3^{n-m+1}\).

不难描述 \(3\) 进制的快速 Fourier 变换, 可以得到复杂度递归式
[ T(N) = L T(M) + O(N \log N), ]
注意最优情况还是 \(M = L = \sqrt{2N}\) 的时候得到的, 递归式解出来复杂度还是 \(O(N\log N \log \log N)\).

以上两个算法有什么用处呢? 重点在于, \(2\) 进制的 Schonhage--Strassen 算法只会除 \(2\), \(3\) 进制的 Schonhage--Strassen 算法只会出现除 \(3\). 对于 \(\mathbb Z\) 上的运算, 如果我们知道 \(x = y / 2^k\), 又知道 \(x = z / 3^m\), 那么考虑 Bezout 系数 \(2^k a + 3^m b = 1\), 就有

\[y = 2^k x, z = 3^m x \implies x = (2^k a + 3^m b)x = ay + bz. \]

所以我们只需要在两个算法的执行过程中都不做除法, 最后把两个算法的结果用 Bezout 系数合并就可以了.

两种思想的结合 --- Furer 算法

Schonhage--Strassen 算法维持了几十年的时间复杂度的最优, 直到 Furer 的突破重新让整数乘法这个问题回到了人们的视线里.

Furer 算法. 整数乘法可以在 \(O(n\log n 2^{O(\log^* n)})\) 时间内完成.

我们现在可以大概感受到, 前面有两种不同的思路有着各自的好处和看似不可弥合的矛盾. 如果使用精度够高的复数或有限域运算, 只需要 \(O(\log n)\) 位就能够提供 \(n\) 次单位根, 但是这样的话, 做快速 Fourier 变换的过程就会需要反复递归调用乘法. 如果造出抽象的单位根, 可以使得快速 Fourier 变换过程中的乘法次数降低, 但这个环本身需要 \(\sqrt n\) 量级的复杂度来构造,
这使得递归的深度更大.

为此, Furer 的思路仍然是考虑精度 \(p = \Theta(\log n)\) 的 \(\mathbb C\), 但在过程中, 他改进了做快速 Fourier 变换的过程.

这里的关键在于, 做 \(p\) 进制的快速 Fourier 变换. 考虑

\[A(X) = \sum_i A_i(X^p) X^i, \]

那么记此时长度为 \(m = \Theta(n/\log n)\), 有

\[A(\omega_m^{j(m/p) + k}) = \sum_i A_i(\omega_{p}^k) \omega_m^{(j(m/p)+k)i}, \]

那么固定 \(k\) 的时候,

\[ \left\{ A(\omega_m^{j(m/p) + k}) \right\}_j = \left\{\sum_i A_i(\omega_{p}^k) \omega_m^{ik} \omega_{p}^{ij} \right\}_j = \mathcal F \left\{ A_i(\omega_{p}^k) \omega_m^{ik} \right\}_i, \]

也就是说, 我们可以不断递归地做长为 \(p\) 的 Fourier 变换.

回忆 Chirp-Z 变换, 或者叫 Bluenstein 算法, 做 Fourier 变换可以转化为 \(O(p)\) 次 \(\mathbb C\) 上的乘法, 以及一次长为 \(2p\) 的多项式乘法. 将它拆包, 就是 \(O(p^2)\) 长度的整数乘法.

显然, 上面过程的复杂度被 \(O(p^2)\) 长度的整数乘法控制, 那么这个做了多少次呢? DFT 的每层我们做了 \(m/p\) 次, 总共有 \(\log_p m = \log m / \log p\) 层, 所以次数是

\[\frac{m \log m}{p \log p} = O\left( \frac{n}{\log n\log \log n} \right). \]

注意, 如果我们直接选取 \(p = O(\log^2 n)\) 精度的复数, 那么得到的调用次数是 \(O(n/\log n)\) 的, 这就是 Furer 算法具体改进了的地方.

接下来, 递归式形如

\[T(n) = O\left( \frac{n}{\log n\log \log n} \right) \cdot T(O(\log^2 n)). \]

记 \(n' = O(\log^2 n)\), 整理一下可以写成

\[\frac{T(n)}{n\log n} = O\left(\frac{T(n')}{n'\log n'}\right). \]

注意 \(f(n) = O(\log^2 n)\) 这个函数迭代 \(k\) 次形如 \(O((\log^{\circ k} n)^2)\), 所以最后迭代到 \(k = \log^* n\), 复杂度就是

\[T(n) = O\left(n\log n 2^{O(\log^* n)}\right). \]

Harvey--van der Hoeven 的高维 Fourier 变换

整数乘法可以在 \(O(n\log n)\) 时间内完成.

Harvey 和 van der Hoeven 的观察来自如下思路: Schonhage--Strassen 算法如果造了一个 \(k\) 次的单位根, 这只能支持做长为 \(k\) 的快速 Fourier 变换. 这导致了只能取 \(k = \sqrt n\). 但是, 如果我们想做高维 Fourier 变换 \(R[X_1,\dots,X_d] / (X_1^k - 1, \dots, X_d^k - 1)\) 的话, 仍然只是需要 \(k\) 次单位根的. 这样就可以取 \(k = n^{1/(d+1)}\) 了.

但遗憾的是, 这样的高维 Fourier 变换是没有用的. 我们知道, 对于 \(n = n_1 \cdots n_d\), 其中 \(n_1,\dots,n_d\) 是互素的整数, 那么中国余数定理告诉我们, 有同构

\[R[X_1,\dots,X_d] / (X_1^k - 1, \dots, X_d^k - 1) \cong R[X]/(X^n - 1). \]

这就意味着, 一维 Fourier 变换可以拆成高维 Fourier 变换, 但是这里有一个关键条件在于 \(n_1,\dots,n_d\) 互素.

在介绍具体算法之前, 让我们先来想象一下, 如果能够弥合上面两个矛盾, 我们可以做什么呢? 我们只需要在 \(O(n\log n)\) 的时间就能够完成快速 Fourier 变换的部分, 然后最后的乘法要做 \(n^{1/(d+1)}\) 次数大小的乘法, 这样的递归式就形如

\[T(n) = O(n^{1 - 1/(d+1)})T(n^{1/(d+1)}) + O(n\log n). \]

注意, 这里的 \(O(n^{1 - 1/(d+1)})\) 里隐藏的常数非常重要. 如果这里常数刚好是 \(d + 1\), 也即

\[T(n) = (d+1) n^{1 - 1/(d+1)}T(n^{1/(d+1)}) + O(n\log n), \]

那么递归每一层的复杂度都是 \(O(n\log n)\), 那么总的复杂度就是 \(O(n\log n)\).

但是假设有正常数 \(\epsilon > 0\), 递归式的常数是

\[T(n) = (d+1 - \epsilon) n^{1 - 1/(d+1)}T(n^{1/(d+1)}) + O(n\log n), \]

那么复杂度的求和第 \(k\) 层就会是 \((1 - \epsilon / (d+1))^{k} \cdot O(n\log n)\), 这个最后还是 \(O(n\log n)\).

在原来的 Schonhage--Strassen 算法里, 递归的常数就是 \(2\), 但相应的递归形式是 \(2\sqrt n T(\sqrt n)\). 接下来的算法中, 我们希望的是确实存在一个类似于 \(T(n) = O(n^{1 - 1/(d+1)})T(n^{1/(d+1)}) + O(n\log n)\) 的递归式成立, 这里的 \(O(n^{1 - 1/(d+1)})\) 内涵的常数和 \(d\) 的选取无关. 这样的话, 只需要固定选取一个充分大的 \(d\), 就能够得到 \(O(n\log n)\) 的复杂度了.

为此, Harvey 和 van der Hoeven 的思路仍然是考虑高维 Fourier 变换的转化, 但是用一些方法重新让低阶的单位根发挥作用.

基于数论假设的算法

在 Furer 算法中, 我们将快速 Fourier 变换用 Bluenstein 算法转化成了卷积, 但这导致卷积的长度更长了. 注意 Bluenstein 算法是对任何等比数列多点求值都转化成卷积的技巧, 对于单位根来说, 我们还有更有性价比的 Rader 变换.

令 \(p\) 是个素数, 考虑 \(p\) 阶快速 Fourier 变换, 也就是说我们要计算

\[\tilde a_j = \sum_{k\in \mathbb Z / p \mathbb Z} a_k \omega^{jk}, \]

我们拆解成

\[\tilde a_j = a_0 + \sum_{k \in (\mathbb Z / p \mathbb Z)^*} a_k \omega^{jk}, \]

由于我们知道 \((\mathbb Z / p \mathbb Z)^*\) 是一个循环群, 记 \(g\) 是一个原根, 我们可以写成

\[\tilde a_{g^j} = a_0 + \sum_{k=1}^{p-1} a_{g^k} \omega^{g^{j+k}}, \]

这就变成了一个 \(p-1\) 阶的循环卷积.

接下来让我们看看如何将已有的东西拼起来. 我们将问题转换成环

\[\mathcal R = \mathbb A[X_1,\dots,X_d] / (X_1^{p_1} - 1, \dots, X_d^{p_d} - 1) \]

上的乘法,
其中

\[\mathbb A = (\mathbb Z / m \mathbb Z)[U] / (U^s + 1). \]

其中 \(s\) 是 \(2\) 的幂.

这些 \(p_i\) 如何选择呢? 答案当然是我们希望每一维的 Fourier 变换都能做. 我们希望协调 \(p_i\), \(m\), \(s\) 的选取使得

• 在 \(\mathbb Z / m \mathbb Z\) 中存在任何 \(p_i\) 阶单位根.
• 但是因为素数阶 Fourier 变换不好做, 我们会需要用 Rader 变换转换为 \(p_i - 1\) 阶循环卷积, 这就需要 \(U\) 能帮助给出 \(p_i-1\) 阶单位根.

接下来的选取是一些技术细节. 首先令 \(\ell\) 也是一个 \(2\) 的幂, 且满足 \(s^{1/3} \leq \ell \leq s^{2/3}\). 让 \(p_i\) 是等差数列 \(\ell + 1, 3\ell+1, 5\ell + 1,\dots\) 上顺序出现的素数, 然后取定

\[v = \operatorname{lcm}(\ell^3, p_1 - 1, \dots, p_d - 1) \cdot p_1 \cdots p_d. \]

我们需要 \(m\) 是一个使得 \(\mathbb Z / m \mathbb Z\) 中存在 \(v\) 阶本原单位根的数.

那么我们现在考虑如何做 Rader 变换.

有人可能要问了, 我们现在是不是每维用 Bluenstein 算法 Fourier 变换就可以了呢? 并不是, 因为每次调用 Bluenstein 算法
都会调用作为子问题的多项式乘法, 这样递归下去就会带一个常数 \(d\), 就失去意义了. 问题的关键在于如何让我们只调用''一次'' Bluenstein 算法.

当做 \(\mathbb A[X_1]/(X_1^{p_1} - 1)\) 上的 Fourier 变换的时候, 我们将它拆成了一个 \(a_0\) 单独的贡献和剩下部分的形如 \(\mathbb A[Z_1]/(Z_1^{p_1-1} - 1)\) 上的循环卷积. 那么对于 \(\mathbb A[X_1,X_2]/(X_1^{p_1} - 1,X_2^{p_2} - 1)\) 而言, 就会有一个 \(a_{00}\), 一个 \(\mathbb A[Z_1]/(Z_1^{p_1-1} - 1)\), 一个 \(\mathbb A[Z_2]/(Z_2^{p_1-1} - 1)\) 和一个 \(\mathbb A[Z_1, Z_2]/(Z_1^{p_1-1} - 1,Z_2^{p_2-1} - 1)\). 总的来说, 在整个 \(d\) 维上, 就有一个大头是

\[\mathcal S = \mathbb A[Z_1,\dots,Z_d] / (Z_1^{p_1-1} - 1, \dots, Z_d^{p_d-1} - 1) \]

上的循环卷积占据了主要的复杂度部分.

设 \(p_i = q_i \ell + 1\), 我们已经预设了 \(q_i\) 是奇数, \(\ell\) 是 \(2\) 的幂, 所以它们互素,
我们又可以拆成

\[\begin{align*} \mathcal S &\cong \mathbb B[Y_1,\dots,Y_d] / \left(Y_1^{\ell}-1,\dots,Y_d^\ell-1\right),\\ \mathbb B &= \mathbb A[V_1,\dots,V_d] / \left(V_1^{q_1}-1, \dots, V_d^{q_d}-1\right). \end{align*} \]

我们的目的是计算 \(\mathcal S\) 上的乘法, 首先它可以在 \(\mathbb B\) 上做 \(\ell^d\) 的快速 Fourier 变换, 由于 \(\mathbb A\) 中的 \(U\) 可以直接提供单位根, 这个过程不需要做任何真正的乘法, 就像 Schonhage--Strassen 所做的那样.

Linnik 常数

可以说, 上面这个算法最后我们不得不面对 \(\mathbb B\) 上的乘法, 而 \(\mathbb B\) 的位数大概是

\[ \underbrace{q_1 \cdots q_d}_{\mathbb B} \cdot \underbrace{s}_{U} \cdot \underbrace{\log m}_{\mathbb Z / m\mathbb Z}. \]

我们需要将它镇压, 使得它可以最后有形如 \(n^{1/d}\) 量级的级别进行递归.

这需要我们在正整数中能找到足够的资料来做这件事. 遗憾的是, van der Hoeven 和 Harvey 并没能无条件地找到这样的资料. 他们的算法需要假设 Linnik 常数充分地好. 我们这里不介绍推导的技术细节, 仅介绍什么是 Linnik 常数.

对于正整数 \(k\) 和与 \(k\) 互素的 \(a\), 我们知道等差数列 \(\{a + kn\}_n\) 中总是存在素数. 记 \(P(k, a)\) 是最小的这样的素数, 然后记 \(P(k)\) 为最大的 \(P(k, a)\).

根据素数分布按照概率的理解, 大家相信实际上有 \(P(k) = O(\varphi(k) \log^2 k)\), 但远远不能证明. 相应地, Linnik 证明了弱一些的事实:

存在常数 \(L > 1\) 使得 \(P(k) = O(k^L)\).

这个常数 \(L\) 就被成为 Linnik 常数. 现在最好的无条件估计是 Xylouris 给出的 \(L=5.18\).

但是我们需要多好的常数呢? Harvey 和 van der Hoeven 的算法需要的是这样的:

如果对于某个 \(L < 1 + \frac 1{303}\) 的 Linnik 定理成立, 那么整数乘法存在 \(O(n\log n)\) 算法.

如果对于某个 \(L < 1 + 2^{-1162}\) 的 Linnik 定理成立, 那么 \(\mathbb F_q\) 上的多项式乘法存在 \(O(n\log q \log (n\log q))\) 算法.

注意, 原则上上面给出的多项式乘法算法还是一个代数算法, 可以理解为, \(n\) 次多项式乘法可以在 \(O(n\log n)\) 的代数运算内完成.
由于 Linnik 定理的强度远远没有达到现在的需求, 所以能否无条件地得到多项式乘法的快速算法还是个开放问题.

无条件整数乘法 --- 数值方法

对于整数乘法, 我们可以自己选取是在什么东西上的多项式, 它既可以是给定精度的复数, 也可以是有限域. 对于有限域, 我们就把问题变成了前述情况, 但是对于复数来说, 我们还有一些潜在的工具没有动用: 那就是数值算法.

记单位环 \(\mathbb T = \mathbb R / \mathbb Z\), 如果我们把 Fourier 变换看做多项式在 \(f(j / n) = F(e^{2\pi i j / n})\), 这可以很自然地将 \(f\) 看做一个定义在 \(\mathbb T\) 上的函数 \(f \colon \mathbb T\to \mathbb C\), 其中 \(f(t) = F(e^{2\pi i t})\).

那么如果说对于稍大的整数 \(m\), 我们求出了所有 \(f(j / m)\), 那对于所有 \(f(j / n)\), 看起来我们有可能根据附近的信息还原出 \(f(j / n)\).

所以, 接下来的想法就是, 先选取差不多大小的素数 \(p_1,\dots, p_d\), 取 \(\mathbb A = \mathbb C[U]/(U^s+1)\), 其中 \(\mathbb C\) 取到 \(p =\Theta(\log n)\) 级别的精度.
然后计算

\[\mathbb A[X_1,\dots,X_d]/(X_1^{p_1} - 1, \dots, X_d^{p_d} - 1) \]

上的乘法. 为了计算利用 \(\mathbb C\) 中单位根的快速 Fourier 变换, 我们用 ''某种方法'' 将大小为 \((s_1,\dots,s_d)\) 的 Fourier 变换规约为稍大的 \((n_1,\dots,n_d)\) 上的 Fourier 变换, 其中 \(n_i\) 皆是 \(2\) 的幂.

这和之前面临的问题一样, 我们还是需要只做一次 Bluenstein 算法, 这只需要我们对 Bluenstein 算法做一点细致的观察. 对 \(n\) 阶 Fourier 变换

\[\tilde a_i = \sum_j a_j \omega^{ij}, \]

为了使用 \(2ij = (i+j)^2 - i^2 - j^2\), 我们需要 \(2n\) 次单位根

\[\begin{align*} \tilde a_i &= \sum_j a_j \omega_{2n}^{(i+j)^2 - i^2 - j^2}\\ &= \omega_{2n}^{-i^2} \sum_j \left(a_j \omega_{2n}^{- j^2}\right) \omega_{2n}^{(i+j)^2}. \end{align*} \]

当 \(n\) 是偶数的时候, 可以发现 \((i+n)^2 \equiv i^2 \pmod {2n}\), 这是可以用 \(n\) 长度的循环卷积计算的. (我们平时很难注意到这一点, 因为平时谁会在能做 FFT 的时候还重复套用一次 Bluenstein 呢?)

所以, 我们就可以将 \((n_1,\dots,n_d)\) 长度的 Fourier 转化为一次

\[\mathbb A[X_1,\dots,X_d]/(X_1^{n_1} - 1, \dots, X_d^{n_d} - 1) \]

的循环卷积, 这样就可以利用 \(\mathbb A\) 本身含有的单位根来计算. 进行真正卷积的规模都只是 \(s\) 级别的.

这样一来, 递归式就确实形如 \(T(n) = O(n^{1-1/(d+1)}) T(n^{1/(d+1)}) + O (n\log n)\) 了. 不过我们还需要确定一个非常关键的问题, 就是如何将 \((p_1,\dots,p_d)\) 的 Fourier 变换规约为 \((n_1,\dots,n_d)\) 的 Fourier 变换. 这便是 Harvey 和 van der Hoeven 的杰作: 我们这里需要用到称为 ''Gauss 重采样'' 的数值算法.

Gauss 重采样

我们推导的起点是 Gauss 函数

\[\rho_s(z) = e^{-\pi z^2 / s} \]

的 Fourier 变换

\[\widehat {\rho_s} (t) = \int_{\mathbb R} \rho_s(z) e^{-2\pi i zt} dz = s^{1/2} \rho_{1/s}(t), \]

而满足适当条件的函数和自身的 Fourier 变换满足 Poisson 求和公式

\[\sum_{k\in \mathbb Z} F(k) = \sum_{k\in \mathbb Z} \widehat F(k), \]

这是下文所作的变换的基础.

这里我们考虑这样的一般情况: 对于互素的 \(s < t\), 如何将 \(s\) 阶的 Fourier 变换规约到 \(t\) 阶的 Fourier 变换上. 注意在实际应用中, \(s\) 是大素数, \(t\) 是 \(2\) 的幂, 所以这个性质显然是满足的.

这里的一个直观是将序列看做 \(\mathbb T\) 上位于 \(j/s\) 和 \(k/t\) 这些分数位置的点.

具体地说, 考虑循环序列 \(\{u_j \in \mathbb C\}_{j\in \mathbb Z / s\mathbb Z}\) (也就是说 \(u_j = u_{j\bmod s}\)), 选定某个实数 \(\alpha > 0\), 我们定义 重采样映射 为

\[ (\mathcal S u)_k = \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2\left(\frac k t - \frac j s\right)^2} u_j, \]

注意将 \(k\) 和 \(k+t\) 对比可以得到

\[\begin{align*} (\mathcal S u)_{k+t} &= \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2\left(\frac {k+t} t - \frac j s\right)^2} u_j\\ &= \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2\left(\frac {k} t - \frac {j+s} s\right)^2} u_j\\ &= \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2\left(\frac {k} t - \frac {j} s\right)^2} u_{j}, \end{align*} \]

所以 \(\mathcal S u\) 是定义在 \(\mathbb Z / t\mathbb Z\) 上的.

那么让我们看看做完重采样之后做 Fourier 变换会怎么样呢?

\[\begin{align*} (\mathcal F\mathcal S u)_l &= \sum_{k} e^{2\pi i kl / t} (\mathcal S u)_k \\ &= \sum_{k \in \mathbb Z / t \mathbb Z} e^{2\pi i kl / t} \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2 \left(\frac k t - \frac j s\right)^2} u_j \\ &= \sum_{k \in \mathbb Z} e^{2\pi i kl / t} \alpha^{-1} \sum_{j\in \mathbb Z / s\mathbb Z} e^{-\pi \alpha^{-2} s^2 \left(\frac k t - \frac j s\right)^2} u_j\\ &= \sum_{j\in \mathbb Z / s \mathbb Z} u_j \alpha^{-1} \sum_{k \in \mathbb Z} e^{2\pi i kl / t} e^{-\pi \alpha^{-2} s^2 \left(\frac k t - \frac j s\right)^2}, \\ \end{align*} \]

我们需要考虑

\[\begin{align*} F(z) &= \alpha^{-1} e^{(2\pi i l / t)z} e^{-\pi \alpha^{-2} s^2\left(\frac z t - \frac j s\right)^2}\\ &= \alpha^{-1} \rho_{\alpha^2} \left( \frac{s z}{t} - j \right) e^{2\pi i l z / t} \end{align*} \]

的 Fourier 变换, 也就有

\[\begin{align*} \widehat F(w) &= \int_{\mathbb R} F(z) e^{-2\pi i w z} dz\\ &= \int_{\mathbb R} \alpha^{-1} \rho_{\alpha^2} \left( \frac{sz}{t} - j \right) e^{2\pi i l z / t} e^{-2\pi i w z} dz\\ &= \frac t s \int_{\mathbb R} \alpha^{-1} \rho_{\alpha^2} (z) e^{2\pi i l z/s} e^{-2\pi i w tz / s} dz\\ &= \frac t s \int_{\mathbb R} \alpha^{-1} \rho_{\alpha^2} (z) e^{2\pi i (l-wt) (z+j)/s} dz \\ &= \frac t s e^{2\pi i (l - wt) j/s} \rho_{\alpha^{-2}}\left(\frac {l-wt} s\right)\\ &= \frac t s e^{2\pi i (l - wt) j/s} e^{-\pi \alpha^{2} t^2 \left( \frac w s - \frac{l/s}{t} \right)^2 } \end{align*} \]

代回原式, 得到

\[\begin{align*} (\mathcal F\mathcal S u)_l &= \sum_{j\in \mathbb Z / s \mathbb Z} u_j \sum_{k \in \mathbb Z} \frac t s e^{2\pi i (l - kt) j/s} e^{-\pi \alpha^{2} t^2 \left( \frac k s - \frac{l/s}{t} \right)^2 }\\ &= \frac t s \sum_{k \in \mathbb Z} \tilde u_{l - kt} e^{-\pi \alpha^{2} t^2 \left( \frac k s - \frac{l/s}{t} \right)^2 }, \end{align*} \]

根据 \(s, t\) 互素, 设 Bezout 系数 \(at+bs=1\), 我们有

\[\begin{align*} &= \frac t s \sum_{k \in \mathbb Z} \tilde u_{t(la -k)} e^{-\pi \alpha^{2} t^2 \left( \frac k s - \frac{l/s}{t} \right)^2 }\\ &= \frac t s \sum_{k \in \mathbb Z} \tilde u_{-tk} e^{-\pi \alpha^{2} t^2 \left( \frac {k+la} s - \frac{l/s}{t} \right)^2 }\\ &= \frac t s \sum_{k \in \mathbb Z} \tilde u_{-tk} e^{-\pi \alpha^{2} t^2 \left( \frac {k} s - \frac{bl}{t} \right)^2 }, \end{align*} \]

这说明, 设另一个 Gau\ss 重采样变换为

\[(\mathcal T u)_k = \sum_{j\in \mathbb Z} e^{-\pi \alpha^2 t^2 \left(\frac k t - \frac j s\right)^2} u_j, \]

我们可以看到,

\[ (\mathcal T \{ (\mathcal F u)_{-tk} \}_k)_{bl} = \sum_{k\in \mathbb Z} \tilde u_{-tk} e^{-\pi \alpha^{2} t^2 \left( \frac {k} s - \frac{bl}{t} \right)^2 }\\ = \frac s t (\mathcal F \mathcal S u)_l. \]

这说明, 先做 Gauss 重采样 \(\mathcal S\), 再做 \(t\) 阶 Fourier 变换 \(\mathcal F\), 和先做 \(s\) 阶 Fourier 变换
再做 Gauss 重采样 \(\mathcal T\), 只相差一个坐标变换.

所以, 我们就把问题拆解成了如下两个:

1. 如何从 \(u\) 出发 (近似地) 快速计算 \(\mathcal S u\)?
2. 如何通过 \(\mathcal T u\) 的结果 (近似地) 快速还原出 \(u\)?

当然, 涉及误差分析的部分实在是非常地技术性, 而且前面的算法我们也都没有做误差分析, 所以我们只能做一个大概的介绍. Harvey 和 van der Hoeven 经过计算得到的结论是, 选取 \(\alpha = \Theta(p^{1/4})\) 的情况下, 可以同时满足高效和足够的精度.

首先是快速计算 \(\mathcal S u\), 注意到求和

\[ (\mathcal S u)_k = \alpha^{-1} \sum_{j\in \mathbb Z} e^{-\pi \alpha^{-2} s^2\left(\frac k t - \frac j s\right)^2} u_j \]

在 \(\left|j - \frac {ks}{t}\right| > L\) 部分的 \(j\) 的贡献只有 \(e^{-\Omega(L^2 \alpha^{-2})}\) 量级,
所以 \(L = \Theta( \alpha \sqrt p )\) 可以得到 \(\Theta(p)\) 位的精度. 这个时候只需要求 \(O(\alpha\sqrt{p})\)
项就达到了误差精度, 所以只花了 \(O(t (\log n)^{3/4})\) 次复数运算, 这比 FFT 直接低一阶, 所以可以接受.

接下来是通过 \(\mathcal T u\) 恢复 \(u\). 注意我们现在持有的实际上是 \(u\) 的一个近似结果, 而且 \(t\) 比 \(s\) 大, 所以我们解方程的过程看起来也是一个数值算法. 首先, 我们可以直接将 \(t\) 个结果保留 \(s\) 个结果, 然后尝试求解. 如何保留呢? 现在观察我们需要求解的形式

\[(\mathcal T u)_k = \sum_{j\in \mathbb Z} e^{-\pi \alpha^2 t^2 \left(\frac k t - \frac j s\right)^2} u_j, \]

由于我们现在分母上是 \(\alpha^2\) 了, 在 \(\left|j - \frac {ks}{t}\right| > L\) 部分的 \(j\) 的贡献只有 \(e^{-\Omega(L^2 \alpha^{2})}\) 量级, 所以 \(L = \Theta(\sqrt p / \alpha)\) 可以得到 \(\Theta(p)\) 位的精度.

直观上看, \(u_j\) 这项贡献到的最大的项是让 \((k/t - j/s)^2\) 最小的 \(k\), 也即

\[k= \left\lfloor \frac{jt}{s} \right\rceil, \]

将对应的最大的那项贡献的 \(k\) 记为 \(k'\), 设 \(\tilde T\) 是选取的这些 \(k'\) 给出的矩阵, 我们将变换写作

\[ (\mathcal T u)_{k'} = e^{-\pi \alpha^2 t^2 \left(\frac{k'}{t} - \frac{j_0}{s}\right)^2} \left( u_{j_0} + \sum_{j \neq 0} e^{-\pi \alpha^2 t^2 \left[\left(\frac{k'}{t} - \frac{j}{s}\right) - \left(\frac{k'}{t} - \frac{j_0}{s}\right)^2\right]} u_j \right), \]

这写成了矩阵分解 \(\tilde T = D (I + \mathcal E)\), 其中 \(D\) 是对角阵, \(\mathcal E\) 是误差项.
如果直接使用展开

\[\tilde T^{-1} = (I - {\mathcal E} + {\mathcal E}^2 - \cdots) D^{-1} \]

来计算会怎么样呢?

Harvey 和 van der Hoeven 计算的结果是, 记 \(\theta = t/s - 1\), 那么 \(\mathcal E\) 的算子范数是 \(\leq 2^{-\Omega(\alpha^2 \theta)}\) 的, 所以需要迭代 \(O(p / \alpha^2 \theta)\) 次才能达到 \(\Theta(p)\) 位的精度. 每次迭代又需要 \(\sqrt p / \alpha\) 次计算, 所以总共的迭代次数是 \(O(p^{3/2} / \alpha^3 \theta)\). 设定参数 \(\theta \geq \Omega(p / \alpha^4)\), 上面这个数就是 \(O(p^{1/2}\alpha)\), 和前述计算 \(\mathcal S\) 的代价一致. 同时, 这个参数并没有让 \(t\) 比 \(s\) 大太多, 因此对于 \(t_1\cdots t_d\) 比 \(s_1 \cdots s_d\) 的估计仍然可以进行. 具体来说, \(\theta\) 可以选定一个 \(1 + O(1/d)\) 级别的常数, 这导致最后 \(t_1\cdots t_d\) 只比 \(s_1\cdots s_d\) 大 \(O(1)\) 倍, 这足够完成递归的复杂度分析.

递归的大头仍然是 \(O(n^{1 - 1/(d+1)}) T(n^{1/(d+1)})\), 上面的数值算法调用的是 \(O\left(n / (\log n)^{1/4}\right)\) 次 \(T(O(\log n))\), 复杂度仍然被 \(O(n\log n)\) 控制.

这实在是一个让人感到惊叹的思想, 一直在之前的研究中显得鸡肋的复数 (之前的算法里复数几乎都可以用有限域替代!), 却帮助我们在临门一脚的时候得到了无条件的结果, 原因是我们真的需要数值分析.

当然, 我还是希望能够看到一个不借助数值分析的算法. 以及我本人可能比较关心的是, 我们有没有可能超越 Cantor 和 Kaltofen 的思路, 得到比 \(O(n\log n\log \log n)\) 更优的 \(\mathbb Z\) 上的多项式乘法?