多项式 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 群找到。