多项式 gcd 的正确姿势：Half-GCD 算法

很久以前就想学这个算法了，但当时迫于知识水平没法理解，就一直咕到现在了，最近重新深刻学习了一下。

本文的大部分内容来自文章 [3] 和 Chee Keng Yap 在 1999 年出版的书 [4]。

我们解决了什么

我们得到了如下问题的 \(\mathcal O(n\log^2 n)\) 做法：

• 计算两个多项式的 gcd
• 给出 \(f(x),h(x)\)，找到 \(g(x)\) 满足 \(f(x)g(x)\equiv 1\pmod{h(x)}\)
• 给出 \(P(x),Q(x)\)，找出两个多项式 \(A(x),B(x)\) 满足 \(A(x)P(x)+B(x)Q(x)=\gcd(P(x),Q(x))\)

换句话说，我们得到了大部分需要用到多项式 gcd 的问题的优秀算法。

通过修改这个算法可以得到最小线性递推式问题的 \(\mathcal O(n\log^2 n)\) 算法，可以参见 [6]。

算法原理

这里以多项式 gcd 为例阐述做法。考虑一下常见的多项式 gcd 的实现：

• 计算 \(R(x)=P(x)\bmod Q(x)\)
• 将问题递归到 \(\gcd(Q(x),R(x))\)。

由于每次递归子问题的时候最劣情况多项式的度数只会减少 \(1\)，这导致我们的复杂度寄了。实际上，不需要精心构造，大部分情况下随机数据都会卡满这个上界。

考虑一下我们的复杂度都浪费在哪里了：计算取模时我们先计算了除法，而计算除法可以只使用最高的若干次项得到答案，这部分的复杂度是对的，而真正破坏复杂度分析的部分是后面的计算取模。纵使除法的答案可能非常小，我们还是使用了 \(\mathcal O(\max(\deg P,\deg Q)\log)\) 的代价计算了取模。

一个简单的想法是，如果能规避掉取模，而是想办法直接得到除法的结果并计算答案就能让复杂度对起来。

什么是 Half-GCD

重新考虑一下多项式 gcd 的过程，实际上是构造了一个多项式序列 \(r_{-2},r_{-1},r_0,r_1,\cdots r_k\)，其中：

• \(r_{-2}=A(x)\)
• \(r_{-1}=B(x)\)
• \(r_i=r_{i-2}-q_ir_{i-1}\)，其中 \(q_i=r_{i-2}/r_{i-1}\)。
• \(r_k=0\)

我们可以不用考虑取模，而是把它写成一个线性变换：

\[\begin{bmatrix} 0 & 1\\ 1 & -q_i \end{bmatrix} \begin{bmatrix} r_{i-2}\\ r_{i-1} \end{bmatrix} = \begin{bmatrix} r_{i-1}\\ r_{i} \end{bmatrix} \]

我们引入记号

\[\lang Q\rang=\begin{bmatrix} Q & 1\\ 1 & 0 \end{bmatrix} \implies \lang Q\rang^{-1}=\begin{bmatrix} 0 & 1\\ 1 & -Q \end{bmatrix} \]

于是无非是

\[\lang q_k\rang^{-1}\cdots \lang q_1\rang^{-1}\lang q_0\rang^{-1} \begin{bmatrix} A(x)\\ B(x) \end{bmatrix} = \begin{bmatrix} \gcd(A(x),B(x)\\ 0 \end{bmatrix} \]

于是我们只需要算出前面那堆东西乘起来的结果 \(\mathbf M^{-1}\)，然后把它作用到我们的初始值上就可以大力算出答案，这样就规避了中间的取模过程，这就是 Half-GCD 算法的基本思想。

具体地，Half-GCD 没有直接算最后的矩阵，它算出的是这样一个矩阵 \(\mathbf M\)：假设 \(\deg A>\deg B\) 且 \(\mathbf M^{-1}\begin{bmatrix}A\\ B\end{bmatrix}=\begin{bmatrix}C\\ D\end{bmatrix}\)，那么必须成立：

\[\deg C\geq \lceil \frac {\deg A}2\rceil>\deg D \]

相信算出了这个矩阵之后该怎么做大家都会：再算一步欧几里得 \(E=C\bmod D\)，于是必然有 \(\lceil \frac {\deg A}2\rceil>\deg D>\deg E\)，于是我们就把问题递归到了一个规模减半的问题，递归调用就可以算出答案。

现在的问题就是如何实现 Half-GCD。

Outline. 这里给出一个 Half-GCD 的形象理解：

假设我们现在要计算 \(\mathrm{hgcd}(A,B)\)，我们设定某个阈值 \(m\)，然后把多项式拆开：

\[A=x^mA_0+A_1 \]

然后计算 \(\mathbf M=\mathrm{hgcd}(A_0,B_0)\)，并直接计算 \(\begin{bmatrix}A'\\B'\end{bmatrix}=\mathbf M^{-1} \begin{bmatrix}A\\B\end{bmatrix}\)。由于前半部分的叙述「除法的结果只和最高次项有关」，我们可以大胆猜测这个变换应用于整体多项式的时候，低次项的影响可以不用考虑，即 \(\deg A'\) 只和 \(\deg A_0'\) 有关，而和 \(A_1\) 无关。更进一步地，在大部分情况下都应该有

\[\deg A'=\deg A_0'+m \]

如果真的有这个性质，那么接下来的思路很简单了：假设原本 \(\deg A=2m\)，那么也就是说我们只需要按照 \(m\) 做一次这种分治就可以得到 \(\deg A'=\frac 32m+o(m)\)，然后我们再按照 \(\frac 12 m\) 做一次这个分治，就可以得到 \(\deg A''=m+o(m)\)，也就是达成了我们的目标，而我们付出的代价仅有两次规模减半的调用，于是复杂度就对了！

下面叙述证明的细节。

Half-GCD 原理

首先我们介绍一个概念：regular 矩阵。

Definition. 称矩阵 \(\mathbf M\) 是 regular 的，当且仅当下面三者中至少一个成立：

• \(\mathbf M\) 是单位矩阵

• 存在某个多项式 \(Q\)，满足 \(\mathbf M=\lang Q\rang\)

• 存在两个 regular 的矩阵 \(\mathbf M_1,\mathbf M_2\)，满足 \(\mathbf M=\mathbf M_1\mathbf M_2\)

换句话说，\(\mathbf M\) 必须是 \(\lang q_1\rang\lang q_2\rang\cdots\) 的形式。

然后有个引理

Lemma 1(ordering property). 设 \(\mathbf M=\begin{bmatrix}p & q\\r &s\end{bmatrix}\) 是一个 regular 的矩阵，则 \(\mathbf M\) 要么是单位矩阵，要么满足：

\[\deg p\geq \max(\deg q,\deg r)\geq \min(\deg q,\deg r)\geq \deg s\\ \deg p>\deg s \]

Proof. 直接归纳即可。\(\blacksquare\)

Definition. 我们称一个 regular 矩阵 \(\mathbf M\) 可以将 \(\begin{bmatrix}A\\B\end{bmatrix}\) reduce 成 \(\begin{bmatrix}C\\D\end{bmatrix}\)，当且仅当

\[\mathbf M^{-1}\begin{bmatrix}A\\B\end{bmatrix}=\begin{bmatrix}C\\D\end{bmatrix} \]

并记做

\[\begin{bmatrix}A\\B\end{bmatrix}\stackrel{\mathbf M}{\longrightarrow}\begin{bmatrix}C\\D\end{bmatrix} \]

于是欧几里得算法的实质就是我们从 \(\begin{bmatrix}A\\B\end{bmatrix}\) 开始，不断用 \(\lang q_i\rang\) reduce 它直到变成 \(\begin{bmatrix}\gcd(A,B)\\0\end{bmatrix}\)。

然后我们直接给出 Half-GCD 算法的伪代码：

\[\begin{array}{ll} &\textbf{Algorithm Polynomial }\operatorname{hgcd}(A,B)\text{:}\\ &\textbf{Input}\text{: }A,B\in\mathbb{F} _ p\lbrack x\rbrack \text{, }\deg(A)\gt \deg(B)\geq 0 \text{.}\\ &\textbf{Output}\text{: a regular matrix }\mathbf{M} \text{ which reduces }(A,B)\text{ to }(C,D)\text{ where}\deg C\geq \lceil\frac{\deg A}{2}\rceil>\deg D.\\ 1&m\gets \left\lceil\frac{\deg(A)}{2}\right\rceil \text{;}\\ 2& \textbf{if }\deg(B)\lt m\textbf{ then return} \left(\begin{bmatrix} 1&0\\ 0&1 \end{bmatrix}\right) \text{;}\\ 3&\mathbf{R}\gets \operatorname{hgcd}(A\operatorname{div} x^m,B\operatorname{div} x^m) \text{;}\\ 4& \begin{bmatrix} C\\ D \end{bmatrix} \gets \mathbf{R}^{-1} \begin{bmatrix} A\\ B \end{bmatrix} \text{;}\\ 5&\textbf{if }\deg(D)\lt m\textbf{ then return }\left(\mathbf{R}\right) \text{;}\\ 6& \begin{bmatrix} Q\\ E \end{bmatrix} \gets \begin{bmatrix} C\operatorname{div}D\\ C\bmod D \end{bmatrix} \text{;}\\ 7&\textbf{if }\deg(E)\lt m\textbf{ then return }(\mathbf{R}\langle Q\rangle )\text{;} \\ 8&k\gets 2m-\deg(D) \text{;}\\ 9&\mathbf{S}\gets \operatorname{hgcd}(D\operatorname{div}x^k, E\operatorname{div}x^k) \text{;}\\ 10&\textbf{return }\left(\mathbf{R}\langle Q\rangle\mathbf{S}\right)\text{;} \end{array} \]

然后分析它的正确性：

Lemma 2(Correctness Criteria). 设 \(A,B\) 为两个多项式且 \(\deg A>\deg B\)，然后设 \(m\) 是某个阈值，\(\mathbf M\) 是一个 regular 矩阵，然后定义如下几个东西：

\[A=x^mA_0+A_1\\ B=x^mB_0+B_1\\ \mathbf M^{-1} \begin{bmatrix} A_0\\ B_0 \end{bmatrix} = \begin{bmatrix} A_0'\\ B_0' \end{bmatrix}\\ \mathbf M^{-1} \begin{bmatrix} A_1\\ B_1 \end{bmatrix} = \begin{bmatrix} A_1'\\ B_1' \end{bmatrix}\\ A'=x^mA_0'+A_1'\\ B'=x^mB_0'+B_1' \]

那么如果满足

\[\deg A_0'>\deg B_0'\\ \deg A_0\leq 2\deg A_0' \]

就一定有如下两个柿子成立：

\[\deg A'=m+\deg A_0'\\ \deg B'\leq m+\max(\deg B_0',\deg A_0-\deg A_0'-1)\\ \]

特别地，上面两个柿子蕴含了

\[\deg A'>\deg B' \]

Proof.

建议先重新读一遍整个引理。

前半部分的一堆定义无非是我们把 \(A,B\) 截取出来，然后对它们分别做变换再合并回去。中间的条件部分可以理解为 \(\mathbf M\) 实际上是 \(\mathrm{hgcd}(A_0,B_0)\) 的产物，然后后面的结论是关于两个多项式的度数的性质。

我们来证明它。设 \(\mathbf M=\begin{bmatrix}P & Q\\R & S\end{bmatrix}\)，首先注意到 \(\deg A_0'>\deg B_0'\) 并且 \(A_0=A_0'P+B_0'Q\)，那么通过 Lemma 1 可以推知

\[\deg A_0=\deg A_0'+\deg P \]

然后又由于 \(\deg A_0\leq 2\deg A_0'\)，所以有

\[\deg P\leq \deg A_0' \]

注意到 \(\mathbf M\) 是 regular 矩阵，这意味着 \(\det \mathbf M=\pm 1\)，于是

\[\mathbf M^{-1}=\pm \begin{bmatrix} S & -Q\\ -R & P \end{bmatrix} \]

于是有 \(A_1'=\pm(A_1S-B_1Q)\)，那么就有

\[\begin{aligned} \deg A_1'&\leq \max(\deg A_1+\deg S,\deg B_1+\deg Q)\\ &\leq \max(\deg A_1,\deg B_1)+\max(\deg S,\deg Q)\\ &<m+\deg P\\ &\leq m+\deg A_0'\\ &=\deg (x^mA_0') \end{aligned} \]

于是我们现在已经得到了第一个子式：

\[\deg A'=m+\deg A_0' \]

第二个子式的证明是类似的。\(\blacksquare\)

这个引理告诉了我们什么：我们取高次项分治，然后把这个线性变换作用于整体多项式，得到的多项式的度数和低次项无关。也就是说，我们在上一节末尾的猜测确实是正确的！接下来只需要应用这个东西两次就可以得到真实答案。

我们来阐述一下具体的细节：

Theorem 1(HGCD Correctness). Algorithm HGCD is correct.

Proof.

我们重新写出我们的 reduce 的路线：

\[\begin{bmatrix}A\\B\end{bmatrix}\stackrel{\mathbf R}{\longrightarrow}\begin{bmatrix}C\\D\end{bmatrix}\stackrel{\lang Q\rang}{\longrightarrow}\begin{bmatrix}D\\E\end{bmatrix}\stackrel{\mathbf S}{\longrightarrow}\begin{bmatrix}D'\\E'\end{bmatrix} \]

返回的矩阵是 regular 是显然的，我们之后只用证明返回的多项式的次数是合法的。

• 第一处 return 语句（第 2 行）

显然度数是合法的

• 第二处 return 语句（第 5 行）

应用 Lemma 2，可以得到 \(\deg C=m+\deg A_0'\geq m\)，于是也是合法的。

• 第三处 return 语句（第 7 行）

由于在第二处没有 return，于是一定有 \(\deg D\geq m>\deg E\)，于是也是合法的。

• 第四处 return 语句（第 10 行）

设 \(m'=\deg A_0=\deg A-m,l=\deg D\)，于是 \(k=2m-l\)。

有个结论：

\[l\geq k\geq 0 \]

这是因为首先有 \(l\geq m\geq m+(m-l)=k\)，然后应用 Lemma 2 还有 \(l=\deg D\leq m+\max(\deg B_0',\deg A_0-\deg A_0'+1)\)，于是 \(l\leq m+\max(\lceil\frac{m'}{2}\rceil-1,\lfloor\frac{m'}{2}\rfloor+1)\leq m+\lfloor\frac{m'}{2}\rfloor+1\)，于是 \(l-m\leq \lfloor\frac{m'}{2}\rfloor+1\leq m\)，那么 \(k=2m-l\geq 0\)。

有了这个结论之后再来考虑一下第二次调用 hgcd 的过程，一定有

\[\deg D_0'\geq \lceil\frac{\deg D_0}{2}\rceil>\deg E_0' \]

而 \(\deg D_0=l-k=2(l-m)\)，于是

\[\deg D_0'\geq l-m>\deg E_0' \]

再次应用 Lemma 2，有

\[\deg D'=k+\deg D_0'\geq k+l-m=m \]

而

\[\begin{aligned} \deg E'&\leq k+\max(\deg E_0',\deg D_0-\deg D_0'-1)\\ &\leq k+\max(l-m-1,l-m-1)\\ &=k+l-m-1=m-1 \end{aligned} \]

于是得证。\(\blacksquare\)

关于复杂度：

Half-GCD 的复杂度显然是 \(T(n)=2T(\frac n2)+\mathcal O(n\log n)=\mathcal O(n\log^2 n)\)

于是 gcd 算法的复杂度就是 \(T(n)=T(\frac n2)+\mathcal O(n\log^2 n)=\mathcal O(n\log^2 n)\)

这是一份参考代码，实现的时候需要尤其注意多项式最高次项为 \(0\) 的 border case。

参考文献

[1] Half-GCD algorithm

[2] HALF-GCD 算法的阐述

[3] HGCD 算法的感性理解 & 伪代码（重置版）

[4] Chee Keng Yap. Fundamental Problems of Algorithmic Algebra.

[5] 彭雨翔: Introduction to Polynomials.

[6] 最短线性递推式求解与有理函数重建

注：[4] 和 [5] 可以在 U 群找到。