多项式模复合的几乎线性算法, 支持多元多项式在线求值的数据结构

本文简要介绍对于有限域 \(\mathbb F_q\), 如何快速计算多项式模复合 \(f(g(X)) \bmod h(X)\), 其中 \(f,g,h\) 均是次数不超过 \(n\) 的多项式. 介绍的思想汇总于 2022 年 Bhargava, Ghosh, Guo, Kumar 和 Umans 的工作: Fast Multivariate Multipoint Evaluation Over All Finite Fields.

为了省略一些技术细节, 我们只讨论 \(\mathbb F_p\) 且 \(p\) 比较大的情况, 也适当略去一些复杂度的分析细节.

第一个想法

首先, 可能每个人遇到这个问题的时候都有一个直接的想法:

1. 对多项式 \(g(X)\) 多点求值, 得到 \(y_i = g(x_i)\).
2. 对多项式 \(f(X)\) 在 \(y_i\) 处多点求值, 得到 \(z_i=f(g(x_i))\).
3. 对 \(z_i\) 进行插值, 还原出多项式 \(f(g(X))\).

但这个想法有个比较严重的错误, 就是如果想要完整的还原出 \(f(g(X))\), 那么注意到它的次数是 \(n^2\) 级别的, 我们需要的不是 \(O(n)\) 个点而是 \(O(n^2)\) 个点, 这就没意义了.

对想法的修正: 变量数和次数

为了让这个想法在插值的这一步变得正确, 我们需要对做求值, 插值的对象做适当的变换. 方式就是通过 增加元数 来 降低次数.

把 \(f(X)\) 这个多项式的上标做 \(d\) 进制分解, 我们可以把多项式 \(f\) 写成另一个多元多项式 \(F(X_0,\dots,X_{m-1})\), 使得 \(F(X,X^d,\dots,X^{d^{m-1}}) = f(X)\).

我们来看看这导致了什么, 现在 \(F\) 关于每个 \(X_i\) 的次数不超过 \(d-1\), 那么总共次数就不超过 \(m(d-1)\), 我们想要计算多项式复合 \(f(g(X))\) 的话, 只需要先计算出 \(g_0 = g \bmod h, g_1 = g^d \bmod h, \dots, g_{m-1} = g^{d^{m-1}} \bmod h\) 的结果, 然后对之前的插值思想如法炮制:

1. 对 \(g_0,\dots, g_{m-1}\) 多点求值, 得到 \(\alpha_i = (g_0(x_i), g_1(x_i),\dots,g_{m-1}(x_i))\).
2. 对 \(F\) 多点求值, 得到 \(z_i = F(\alpha_i)\).
3. 对 \(z_i\) 多点插值, 还原出 \(F(g_0,\dots,g_{m-1})\).

我们看看现在我们需要多少个点值, \(F(g_0,\dots,g_{m-1})\) 的次数是 \(O(nmd)\) 的, 而 \(m = \log_d n\), 所以只要能够对于某个 \(d=n^{o(1)}\) 级别的 \(d\) 给出一个高效的多点求值算法, 我们就能够在 \(n^{1+o(1)}\) 的时间内完成多项式模复合了.

但后者其实不是一个平凡的问题, 大家一直没有找到 \(n^{1+o(1)}\) 的代数算法, 直到将目光放到 非代数算法 的路径上.

间奏: 代数算法的下界

我们知道, 多点求值有一个 \(O(n\log^2 n)\) 的优美的算法, 它的特点不仅在于简单而且真的可以跑, 还在于它是一个 代数算法.

我们不太精确地勾勒一下什么是代数算法. 一般来说, 代数算法考虑的不是一个图灵机, 而是一个代数电路, 电路的基本单元是无限域 \(K\) 上的四则运算.

让 \(K\) 无限是非常重要的, 如果 \(K\) 有限, 那么它的功能和布尔电路相当, 而布尔电路的下界是相当棘手的, 目前虽然有一些进展, 但看起来距离一些本质问题的答案有很远的距离. 而当 \(K\) 无限的时候, 我们似乎能够通过求助一些代数不变量来刻画问题的困难程度.

一个令人困惑的问题是这样的: 如果我们可以 \(O(n\log^2n)\) 一次性求出 \(n\) 次多项式的 \(n\) 个点值, 那我们可以多快在线地求出一个点值呢?

如果只考虑代数算法, 有一个令人沮丧的答案: 每次必须花 \(O(n)\) 的时间. 我们这里给出一个粗糙的勾勒, 说明为什么是这样:

我们的算法应当有两个 phase, 在第一个 phase 里, 我们接受多项式的所有系数, 然后对此进行任意有限多的预处理. 在第二个 phase 里, 我们接受一个输入 \(x\), 然后需要输出 \(f(x)\) 的值.

考虑我们建出了这个代数电路, 如果它确实能输出结果, 会怎么样呢? 我们从输出答案的那个单元往回追溯, 并且一路标记出所有 取值和 \(x\) 无关 的单元, 如果这样的单元有 \(m\) 个, 记它们的输出为 \(q_1,\dots,q_m\), 那本质上我们的答案可以表示成一个有理函数 \(Q(x, q_1,\dots,q_m)\), 其中每个 \(q_i\) 都是关于多项式的系数 \(a_0,\dots,a_n\) 的一个有理函数. 那么 \((a_0,\dots,a_n)\mapsto (q_1,\dots,q_m)\) 应该是个单射, 但这只有在 \(m\geq n+1\) 的时候才是可能的. (严格地论述它还需要一些功夫, 但至少这是符合直觉的: 代数运算不应该能够 "压缩信息").

你会发现, 如果我们对于有限域的情况的求值过程太机械, 就也会受到前述论证的下界限制. 后面我们会看到, 打破这个限制的一个原因在于 寻址, 这极大程度扩大了求值过程中输入信息的信息量.

算法 1: CRT 约化

回顾我们现在的问题: 给定 \(m\) 元, 关于每个元不超过 \(d-1\) 次的, \(\mathbb F_p\) 上的多项式 \(F\), 希望求出它在 \(N\) 个点处的点值.

首先我们应该知道在何种情况上我们已经知道多点求值是简单的, 那就是点集 \(= [p]^m\) 的情况, 这是可以每个变量做一次 FFT 就能求出来的. 但问题在于, 我们的 \(p\) 可以很大, 所以 \(p^m\) 实际上是很大的.

注意到, 我们现在的 \(d = N^{o(1)}\), 是低次多项式, 我们将再次用到这个好处: 低次 \(\implies\) 更少的插值. 这使得我们可以考虑先把 \(F\) 看做一个 \(\mathbb Z\) 上的东西, 求值的点也是 \(\mathbb Z\) 上的.

1. 我们选取最小的一系列素数 \(p_i\), 使得它们的乘积覆盖了可能结果的上界.
2. 将问题约化成 \(\bmod p_i\) 的各运算, 将它们各自解决
3. 通过 CRT 还原出作为 \(\mathbb Z\) 上求值的结果, 再对 \(p\) 取模.

我们看看 \(\bmod p_i\) 的 \(p_i\) 会多大. 假设一开始的模数是 \(M\), 由于 \(F\) 的总次数 \(< md\), 单点求值的结果不会超过 \(d^m p^{md}\).

一个基本的事实是, 数值在 \(O(\log N)\) 以内的素数之乘积大于 \(N\), 所以使用 CRT 的素数规模是 \(O(md\log p)\) 的. 注意, 如果 \(\log p\) 仍然过于大了, 我们这个过程可以继续递归进行, 直到所有叶子节点的 \(p_i\) 大小大概是 \(x = md\log x\) 的解的范围, 也即 \(O(md\log (md))\).

所以在 \(p_i = O(md\log (md))\) 的时候, 求所有点值的复杂度就变成了 \(O(md\log (md))^d\). 为了让它 \(= (m^d)^{1+o(1)}\), 需要令 \(m = d^{o(1)}\).

经过一定仔细的分析, 可以证明在 \(m = d^{o(1)}\) 的时候, 这个算法可以在 \((d^m+N)^{1+o(1)}\operatorname{poly}(m,d,\log p)\) 的时间内完成多元多项式的多点求值.

接下来我们完成我们的另一个承诺: 它能帮我们做在线多点求值. 首先注意到这个算法本身预处理了一个巨大的数表, 我们每次查询一个点值的时候只需要按照 CRT 约化的树递归下去, 然后不断 CRT, 取模就行了. 经过分析可以证明的是, 调用它本身的时间是 \(\operatorname{poly}(m,d,\log p)\) 的.

但我们有 \(m=d^{o(1)}\) 的限制, 对于一元多项式我们可以故技重施, 将它拆成多元低次的多项式, 这样就完成了在线求值.

comments

值得一提的是, CRT 转化还曾经用到过有限域上多项式乘法的 双线性复杂度 (bilinear complexity) 的研究中. 粗略地说, 双线性复杂度只关心把输入 \(a\) 和 \(b\) 之间进行乘法的次数. 这个复杂度足够有效刻画矩阵乘法了, 但是对多项式乘法来说比较奇怪: 在这个模型下, 如果 \(K\) 是无限域, 多项式乘法只需要 \(n\) 次双线性乘法: 求值, 乘法, 再做插值即可.

但对于有限域而言, 直接求值的点数有限制, 所以需要从小到大枚举不可约多项式, 对所有取模的结果做 CRT, 然后递归成更小的, 直到全部递归成一次项. 这样会得到一个非常接近线性的双线性复杂度.

算法 2: 模幂下 Taylor 展开

前述的算法 1 其实是 Kedlaya 和 Umans 在 2011 年的工作, 而 2022 年的一个进展是将 \(m=d^{o(1)}\) 这个限制改成了一个更弱的限制.

思路的起点在于如何设计一种策略, 让我们选取一个更小的集合 \(S\), 使得求出全部 \(S\) 内的点值是容易的, 同时对于 \(\alpha \notin S\), 也有办法比暴力求值更快地还原出 \(F(\alpha)\) 的值.

这件事如果不考虑 \(\bmod p\) 而是 \(\bmod r^s\) 的话则可以有如下的考虑. 我们选取 \(S=[r]^d\), 这种乘积集合的求值总是简单的, 接下来, 如果 \(\alpha - \overline{\alpha}\) 是 \(r\) 的倍数的话, 我们考虑 \(F(\overline{\alpha} + (\alpha - \overline{\alpha}))\) 做 Taylor 展开 (在原论文中称为考虑 Hasse 导数), 我们其实可以用 \(F\) 在 \(\overline{\alpha}\) 的 \(<s\) 阶导数的值来还原出 \(F(\alpha)\) 的值.

具体来说, 递归过程将算法 1 修改如下 (假设要求的计算结果 over \(\bmod r^s\), 且 \(s\leq m\)):

1. 选取最小的一些列素数 \(p_i\), 使得 \(\prod p_i \geq dr^d\).
2. 对于每个 \(e_1 + \cdots + e_m < s\), 计算出 Hasse 导数多项式 \(\frac{\partial_1^{e_1}\cdots \partial_m^{e_m}}{e_1!\cdots e_m!}F\).
3. 将问题约化成 \(\bmod p_i^m\) 的各运算, 且对每一个导数多项式都做在约化到 \(\overline{\alpha}\in[r]^m\) 上的多点求值.
4. 将计算的结果用 CRT 和 Taylor 展开式合并出答案.

接下来我们还是粗略地说明正确性和时间复杂度.

注意到约化到的点值 \(<r\), 而系数 \(< r^s\), 所以作为 \(\mathbb Z\) 上的答案 \(<d^m \cdot r^m \cdot r^{m(d-1)} \leq (dr^d)^m\), 这说明 \(\bmod p_i^m\) 的所有答案还原出的值域能覆盖答案.

相比于只用 CRT, 我们现在每次递归调用的 \(r\) 的规模从 \(r\) 变成 \(O(\log (dr^d)) = O(d\log r)\), 所以叶子节点有 \(r = O(d\log d)\).

另外值得注意的是我们递归的时候每次不只是分成几个素数就有几个岔, 还对每个导数多项式都递归下去了, 这样的导数是 \(\sum e_i < m\) 的解数, 也就是 \(\binom{2m-1}{m-1} = 2^{O(m)}\).

为了让这个不太爆炸, 需要对递归树的深度进行相当的限制, 注意到递归的目的大概是每次让 \(r\) 的大小取对数, 所以为了让递归的层数是 \(O(1)\), 需要 \(r \leq \underbrace{\exp (\exp (\cdots (\exp}_{O(1)} (d))))\), 或者写作 \(\operatorname{slog}(r) \leq \operatorname{slog}(d)+O(1)\).

此外, 这个算法也不能处理常数大小的 \(d\), 因为它的复杂度的对应项是 \(O(d\log d)^{m}\), 但只要是 任意缓慢增长 的 \(d\), 都有 \(O(d\log d)^{m} = (d^m)^{1+o(1)}\).

这种算法显然也是支持在线求值的, 但是复杂度项里会再带个 \(O(1)^m\).

算法 3: 有限域上的挂谷集合

另一个策略则仍然在有限域上考虑. 假设集合 \(K\) 满足这样的性质: 虽然 \(\alpha\notin K\), 但是我们可以将 \(\alpha\) 参数化到一个低次曲线 \(\alpha\in \Gamma = \{ C(t) | t\in \mathbb F \}\) 上, 并且 \(\Gamma\) 交 \(K\) 于充分多的点.

这有什么用? 注意到如果我们想求的是一个多项式 \(F\) 的值, 由于 \(F(C(t))\) 对于某个 \(t\) 的值, 而 \(F(C(t))\) 是一个低次多项式, 我们又知道很多其它位置的值, 那就可以插值出多项式 \(F(C(t))\), 然后带入得到 \(F(\alpha)\) 了.

挂谷集合原本的定义是在 \(\mathbb R^2\) 上一个能塞进任何一个方向的单位长度线段的集合, 这个大概科普文章网上满天飞了. 有限域上原本的挂谷集合定义则是每个方向都完整包含一条直线. 如果 \(C\) 是一次多项式并且要求 \(\Gamma-\alpha \subset K\), 这种集合被称作 Nikodym 集合, 这两种集合都被 Zeev Dvir 利用多项式方法 (polynomial method) 的出色工作证明了 \(|K| \geq c_n |\mathbb F|^n\). 所以我们的问题里不可能要求这个集合的交集有 \(q\) 这么大.

构造

这个构造来源于 Björklund, Kaski 和 Williams 在 2019 年的另一篇工作: Generalized Kakeya Sets for Polynomial Evaluation and Faster Computation of Fermionants. 他们利用这个技术把有界矩阵的积和式的计算从 \(2^{m - \Omega(\sqrt{m/\log m})}\) 优化到了 \(2^{m - \Omega(\sqrt{m/\log\log m})}\).

令设有限域的大小为 \(q\), 并且整数 \(u\) 满足 \(u+1\) 整除 \(q-1\),
考虑集合

\[S_\tau = \left\{ \left(\frac{x}{u+1}+\tau\right)^{u+1} - \tau^{u+1} \middle | x \in \mathbb F_q \right\}, \]

令

\[K = \bigcup_{\tau \in \mathbb F_q} S_\tau^m. \]

我们首先分析 \(K\) 的大小. 注意到 \(u+1\) 是 \(q-1\) 的一个因子, 所以一个数的 \(u+1\) 次幂本质上只有 \((q-1)/(u+1) + 1\) 种可能性. 所以我们有

\[|K| \leq \left(\frac{q-1}{u+1} + 1\right)^{m+1}. \]

接下来, 真正的做法并不是直接构造一条曲线穿过 \(\alpha\), 而是考虑展开定义式, 定义

\[G_\alpha(t) = \left(\left(\frac{\alpha_i}{u+1}+t\right)^{u+1} - t^{u+1}\right)_{i=1}^m, \]

显然, \(G_\alpha\) 完全经过 \(K\), 而且最高次项是 \(u\) 次, 系数是 \(\alpha_i t^u\). 我们称这种集合为广义挂谷集合 (generalized Kakeya set).

算法

接下来, 实际上进行值的还原的策略其实是这样: 先把多项式 \(F\) 拆成各个同次部分 \(F = F_0 + F_1 + \cdots + F_{m(d-1)}\), 而对于 \(F_i\), 多项式 \(F_i(G_\alpha)\) 的最高次项是 \(iu\) 次, 系数给出了 \(F_i(\alpha)\) 的值.

我们的目的是让 \((q-1)/(u+1)+1\leq d\), 设 \(\tilde d \leq d\), 且 \(\tilde d-1 \mid q-1\), 我们取 \(u = q/(\tilde d+1) - 1\). 这样我们求值的点的范围就是 \(|K| \leq \tilde d^{m+1} \leq d^{m+1}\) 了.

但下一个锅来了: \(u\) 的次数其实挺大的, \(iu\) 可以有多大呢? \(i \leq m(d-1)\), 所以 \(iu\) 可以有 \(2(d/\tilde d)mq\) 这么大, 我们可以要求 \(\tilde d = \Omega(d)\), 但依然插 \(q\) 个点值都不够用!

这个时候, 又轮到 Hasse 导数出马了, 我们其实需要求出 \(K\) 上每个点的 Hasse 导数直到 \(O(m)\) 阶, 这样就知道了 \(F_i(G(t))\) 在每个点处的更高阶导数信息, 通过 Hermite 插值还原出足足 \(O(mq)\) 项的 \(F_i(G(t))\).

这样一来, Hasse 导数让复杂度膨胀了 \(\binom{m+O(m)}{m} \leq O(1)^m\) 倍, 但我们拯救了求值的点集大小. 我们只需要把算法 1 的底层叶子节点换成这个算法, 复杂度就变成了 \((O(d)^m + N\cdot O(1)^m) \cdot \operatorname{poly}(m, d, q)\), 注意 \(\operatorname{poly}\) 项现在有了 \(q\), 现在底层的 \(q\) 虽然不需要太小, 但也不能太大了.

来自解析数论的修正

我们还没有解决一个问题, 就是需要让 \(\tilde d = \Omega(d)\) 是可能的. 这件事其实可以反过来考虑, 就是我们先选取 \(\tilde d\), 然后寻找满足 \(p\bmod (\tilde d-1) = 1\) 的素数.

我们知道 Dirichlet 定理指出, 在素数中, \(\bmod x=1\) 的素数占比趋近于 \(1/\varphi(x)\), 但可惜的是这是对每个单独的 \(x\) 给出的结果, 不能一致地控制 \(x\).

所以, 对于本算法的分析而言, 需要请出的是更强的 Bombieri–Vinogradov 定理, 它给出了一定范围的 \(x\) 内的一个一致控制.

这里的技术细节我们略去, 最后得出的结论大概是, 对于任何充分大的 \(d\), 选取一个在 \([0.8d, d]\) 内的 \(\tilde d\), 就可以找到一系列 \(p_i \bmod \tilde d=1\), 使得 \(\prod p_i > M\), 且 \(p_i \leq O((D\log M)^3)\).

遗留问题

现在对于 任意增长 的 \(d\), 挂谷集合的做法漂亮地做到了我们要求的复杂度, 进一步地, 我们可以对于任意 \(\epsilon > 0\), 存在一个常数 \(d\), 使得对于该 \(d\) 来说复杂度是 \((d^m + N) \cdot d^{\epsilon m} \cdot \operatorname{poly}(m, d, q)\).

目前的做法似乎还不能本质上将类似 \(d=2\) 的情况 (多重线性型) 的复杂度压到近线性, 或许是一个值得思考的方向.