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

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

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

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

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

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

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

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

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

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

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

我们引入记号

于是无非是

于是我们只需要算出前面那堆东西乘起来的结果 \(\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}\)，那么必须成立：

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

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

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

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

然后计算 \(\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=2m\)，那么也就是说我们只需要按照 \(m\) 做一次这种分治就可以得到 \(\deg A'=\frac 32m+o(m)\)，然后我们再按照 \(\frac 12 m\) 做一次这个分治，就可以得到 \(\deg A''=m+o(m)\)，也就是达成了我们的目标，而我们付出的代价仅有两次规模减半的调用，于是复杂度就对了！

下面叙述证明的细节。

首先我们介绍一个概念：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\) 要么是单位矩阵，要么满足：

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

Definition. 我们称一个 regular 矩阵 \(\mathbf M\) 可以将 \(\begin{bmatrix}A\\B\end{bmatrix}\) reduce 成 \(\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 算法的伪代码：

然后分析它的正确性：

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

那么如果满足

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

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

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\leq 2\deg A_0'\)，所以有

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

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