多项式欧几里得算法解析
warzone
·
2025-07-20 13:42:38
·
算法·理论
给定两个有限域上的多项式 f,g\in\mathbb{F}[x] （|\mathbb{F}| 为素数），即便使用 FFT，\Theta(n^2\log n) ，\Theta(n\log^2 n) 。
不论是 EI 还是 hly1204，都只对折半欧几里得算法作了大致上的理解，
本文章是我自己对折半欧几里得算法原理及实现细节的分析。
朴素的扩展欧几里得与多项式除法过程分析
在讲解折半欧几里得算法之前，我们先回顾一下朴素的扩展欧几里得算法与多项式除法的执行过程。
给定 f(x),g(x) ，欲求解 \gcd(f,g) ，则有迭代式
p_0(x)=f(x),p_1(x)=g(x),p_{k+1}(x)=p_{k-1}(x)\bmod p_k(x)
经过若干次迭代后，即可得 p_k(x)=\gcd(f,g) 且 p_{k+1}(x)=0 。这就是最基本的欧几里得算法。
扩展欧几里得算法则不止于此，其还希望给出 a(x),b(x) 使得
a(x)f(x)+b(x)g(x)=\gcd(f,g)
其计算过程是可以用矩阵表示的，具体地。设 q_k(x)=\left\lfloor\dfrac{p_{k-1}(x)}{p_k(x)}\right\rfloor ，并记
\left<q_k\right>=\begin{bmatrix}
0&1\\
1&-q_k
\end{bmatrix}
则
{p_k\brack p_{k+1}}=\left<q_k\right>{p_{k-1}\brack p_k}=\left<q_k\right>\left<q_{k-1}\right>\cdots\left<q_1\right>{p_1\brack p_0}
故若 p_k=\gcd(f,g) ，则矩阵 \textbf{M}(x)=\displaystyle\left<q_k\right>\left<q_{k-1}\right>\cdots\left<q_1\right> 就是扩展欧几里得算法的执行结果。
不难发现，虽然欧几里得算法执行过程中不断对多项式取余数，获得扩展欧几里得算法的执行结果只需要知道商而无需知道余数 。
此外，若 \deg f\ge \deg g ，则以下两式保证了扩展欧几里得算法的中间结果和最终结果都不会太长：
\deg q_k=\deg p_{k-1}-\deg p_k
\deg\textbf{M}=\sum_{i=1}^{k}\deg q_k=\deg p_0-\deg p_k=\deg f-\deg\gcd(f,g)
接下来回顾多项式除法的执行过程。n 次多项式 f 和 m 次多项式 g\ (n\ge m) ，q(x)=\left\lfloor\dfrac{f(x)}{g(x)}\right\rfloor ，f(x),g(x),q(x) 系数翻转后的结果分别为 f_R(x),g_R(x),q_R(x) ，则
f_R(x)\equiv g_R(x)q_R(x)\pmod{x^{n-m+1}}
也就是说，商 q(x) 只与 f(x),g(x) 最高的 n-m+1 个项有关 。f(x),g(x) ，f(x),g(x) 的高次项来计算一部分迭代。这就是折半欧几里得算法的核心思想。
截断对扩展欧几里得的影响
截断之后，我们并不能完全计算出所有的 q_k(x) ，而只能计算出前面的一部分。
F(x)=f(x)x^N+u(x)\quad(\deg u<N)
G(x)=g(x)x^N+v(x)\quad(\deg v<N)
其中 f(x),g(x) 已知，N,u(x),v(x) 未知，且 \deg f\ge \deg g ，F(x),G(x) 扩展欧几里得的结果（即每一次辗转相除的商）。设
P_0=F,P_1=G,P_{k+1}=P_{k-1}\bmod P_k,Q_k=\left\lfloor\dfrac{P_{k-1}}{P_k}\right\rfloor
P_k(x)=p_k(x)x^{N_k}+r_k(x)\quad(\deg r_k<N_k)
其中 p_k(x) 是已知的，而 N_k,r_k(x) 是未知的（当然，N_{k-1}-N_k 是已知的），则
$$\deg Q_k=\deg P_{k-1}-\deg P_k\le \deg P_k-N_k=\deg p_k$$
此时 $P_{k-1},P_{k}$ 的已知项足以让我们求得 $Q_k$。且
$$N_{k+1}-N_k=\deg Q_k=\deg P_{k-1}-\deg P_k$$
即，**迭代在消去高次的已知项的同时，在低次将相同数量的已知项变为了未知项**。
$$\deg Q_k=\deg P_{k-1}-\deg P_k>\deg P_k-N_k=\deg p_k$$
此时已知项太少，不足以求得 $Q_k$。
综上，在截断的情况下计算扩展欧几里得，
\deg P_k-N<\dfrac{\deg f}{2}\le\deg P_{k-1}-N\quad(1)
若 N=0 ，则该情形下再做一次辗转相除，问题规模将缩减至原来的一半。
折半欧几里得算法的执行过程
根据上述讨论，折半欧几里得算法需要实现一个函数 \operatorname{halfGcd}(f,g) ，f,g 的情况下计算 F,G 的扩展欧几里得，直至达到式 (1) ，Q_1,Q_2,\cdots,Q_k 所组成的，\textbf{M}(x)=\left<Q_1\right>\left<Q_2\right>\cdots\left<Q_k\right> 。
在实现了 \operatorname{halfGcd} 的情况下，
计算 \operatorname{halfGcd}(F,G) ，达到式 (1) 且 N=0 。
做一次辗转相除，问题规模缩减至原来的一半。
该过程的伪代码如下：
\begin{matrix}
& \textbf{Algorithm Polynomial } \operatorname{polyGcd}(F,G):\qquad\qquad\quad\quad\ \\
& \textbf{Input}: F,G\in \mathbb{F}[x],\deg F>\deg G\ge 0.\qquad\qquad\qquad\quad\ \\
& \textbf{Output}: \text{a matrix } \textbf{M}(x) \text{ such that} \displaystyle{\gcd(F,G)\brack 0}=\textbf{M}{F\brack G}&\\
&(\deg \textbf{M}＝\deg F-\deg\gcd(F,G)).\quad\\
1&\textbf{M}\leftarrow\operatorname{halfGcd}(F,G);\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
2&\displaystyle{A\brack B}=\textbf{M}{F\brack G};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \\
3&\textbf{if } B=0 \textbf{ then return M};\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
4&\displaystyle{Q\brack P}\leftarrow\begin{bmatrix}
\left\lfloor\dfrac{A}{B}\right\rfloor\\
A\bmod B
\end{bmatrix};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
5&\textbf M\leftarrow\left<Q\right>\textbf{M};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \ \\
6&\textbf{if } P=0 \textbf{ then return M};\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
7&\textbf{return}\left(\operatorname{polyGcd(B,P)}\textbf{M}\right);\qquad\qquad\qquad\qquad\qquad\qquad\quad\
\end{matrix}
若 \operatorname{halfGcd} 的时间复杂度为 \Theta(n\operatorname{polylog}(n)) ，则多项式欧几里得的时间复杂度由主定理得为
T(n)=T\left(\dfrac{n}{2}\right)+T_{\operatorname{halfGcd}}(n)+\Theta(n\log n)=\Theta(T_{\operatorname{halfGcd}}(n))
接下来考虑如何实现 \operatorname{halfGcd}(f,g) 。不难发现其也可以以折半的方式实现：
令 m=\left\lceil\dfrac{\deg f}{2}\right\rceil ，计算 \operatorname{halfGcd}\displaystyle\left(\left\lfloor\dfrac{f}{x^m}\right\rfloor,\left\lfloor\dfrac{g}{x^m}\right\rfloor\right) ，作用到 \displaystyle{f\brack g} 上再执行一次辗转相除，\displaystyle{s\brack p} ，则 \deg s<\dfrac{3}{4}\deg f 。这就是前半段迭代。
令 k=2m-\deg s ，则 \left\lfloor\dfrac{s}{x^k}\right\rfloor,\left\lfloor\dfrac{p}{x^k}\right\rfloor 就是 s(x),p(x) 当前状态下的已知项，
该过程的伪代码如下：
\begin{matrix}
& \textbf{Algorithm Polynomial } \operatorname{halfGcd}(f,g):\qquad\qquad\qquad\qquad\qquad\ \ \\
& \textbf{Input}: f,g\in \mathbb{F}[x],\deg f>\deg g\ge 0.\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
& \textbf{Output}: \text{a matrix } \textbf{M}(x) \text{ such that}\deg b<\dfrac{\deg f}{2}\le \deg a\le \deg f\\
&\text{ where }\displaystyle{a\brack b}=\textbf{M}{f\brack g}\ \ (\deg \textbf{M}＝\deg f-\deg a).\\
1&m\leftarrow\left\lceil\dfrac{\deg f}{2}\right\rceil;\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \\
2&\textbf{if}\deg g<m\textbf{ then return}
\begin{bmatrix}1&0\\0&1\end{bmatrix};\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \ \\
3&\textbf{M}\leftarrow\operatorname{halfGcd}\displaystyle\left(\left\lfloor\dfrac{f}{x^m}\right\rfloor,\left\lfloor\dfrac{g}{x^m}\right\rfloor\right);\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \\
4&\displaystyle{r\brack s}=\textbf{M}{f\brack g};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \\
5&\textbf{if }^{^{^{}}} \deg s<m \textbf{ then return M};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \ \ \\
6&\displaystyle{q\brack p}\leftarrow\begin{bmatrix}
\left\lfloor\dfrac{r}{s}\right\rfloor\\r\bmod s\end{bmatrix};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \\
7&\textbf M\leftarrow\left<q\right>\textbf{M};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \ \\
8&\textbf{if } \deg p<m \textbf{ then return M};\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \\
9&k\leftarrow 2m-\deg s;\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \ \ \\
10&\textbf{return}\left(\operatorname{halfGcd}\displaystyle\left(\left\lfloor\dfrac{s}{x^k}\right\rfloor,\left\lfloor\dfrac{p}{x^k}\right\rfloor\right)\textbf{M}\right);\qquad\qquad\qquad\qquad\qquad\quad\ \ \ \
\end{matrix}
因此 \operatorname{halfGcd} 的时间复杂度由主定理得为
T_{\operatorname{halfGcd}}(n)=2T_{\operatorname{halfGcd}}\left(\dfrac{n}{2}\right)+\Theta(n\log n)=\Theta(n\log^2 n)
综上，折半欧几里得算法的时间复杂度为 \Theta(n\log^2 n) 。