FLINT3.7.0-dev
/
Algorithm

Half-GCD

Cost
O(M(n) \log n) field operations for the half-gcd, GCD, extended GCD and resultant of polynomials of degree n, against O(n^2) for the Euclidean algorithm.

The Euclidean algorithm computes \gcd(a, b) for polynomials of degree at most d with O(d^2) operations in the coefficient field. The half-gcd algorithm reduces this to O(M(d) \log d), where M(d) is the cost of a multiplication. It rests on the observation that the first quotients of the Euclidean algorithm depend only on the leading coefficients of the inputs, so that about half of the quotients can be computed recursively from the top halves of a and b. The idea is due to Knuth [Knu1970] and Schönhage [Sch1971] for integers and continued fractions, and to Moenck [Moe1973] for polynomials; a correct and complete version for polynomials was given by Brent, Gustavson and Yun [BGY1980], and Thull and Yap [ThullYap1990] gave a unified treatment for polynomials and integers. A textbook account is [GG2003, Ch. 11]. Thomé [Tho2002] generalised the recursion to matrix polynomials.

Throughout, the coefficients lie in a field K; for nmod_poly this means that the modulus must be prime (more precisely, every leading coefficient that occurs must be invertible).

Generic version#

The half-gcd algorithm is implemented once, in gr_poly, for any field in which the leading coefficients that occur can be inverted; a ring \mathbb{Z}/N\mathbb{Z} pretending to be a field (gr_ctx_set_is_pretend_field) runs it as well, and a non-invertible leading coefficient then reveals a factor of N. The generic GCD, extended GCD and resultant use it over finite fields from length 200, with inner cutoff 100, and the Euclidean algorithm over other fields, where the half-gcd algorithm suffers from coefficient growth (exact infinite fields such as \mathbb{Q}) or numerical instability (floating-point and ball arithmetic). The type-specific modules call the same functions with their own cutoffs (Specialised algorithms); those of nmod_poly are listed in Use in FLINT. The 2 \times 2 matrix products are the ring method _gr_poly_hgcd_mat_mul, so a ring can replace the products described below. Polynomial GCDs and resultants compares the half-gcd algorithm with the Euclidean and subresultant algorithms over all rings.

The remainder sequence and its matrices#

For a, b \in K[x] with \deg a > \deg b, the Euclidean algorithm produces remainders r_0 = a, r_1 = b and

r_{i-1} = q_i r_i + r_{i+1}, \qquad \deg r_{i+1} < \deg r_i,

until some r_{k+1} = 0; then r_k is a GCD. With the quotient matrices Q_i = \begin{pmatrix} q_i & 1 \\ 1 & 0 \end{pmatrix} we have (r_{i-1}, r_i)^T = Q_i (r_i, r_{i+1})^T, so

\begin{pmatrix} a \\ b \end{pmatrix} = M_j \begin{pmatrix} r_j \\ r_{j+1} \end{pmatrix}, \qquad M_j = Q_1 Q_2 \cdots Q_j, \qquad \det M_j = (-1)^j.

Since \deg r_{i-1} = \deg q_i + \deg r_i, every q_i has degree at least 1, and the entries of M_j have degree at most \sum_{i \le j} \deg q_i = \deg a - \deg r_j.

Lemma 1 (Uniqueness). Suppose that (a, b)^T = Q_1 \cdots Q_j (A, B)^T for polynomials q_1, \ldots, q_j of degree at least 1, with \deg A > \deg B. Then q_1, \ldots, q_j are the first j Euclidean quotients of (a, b), and A = r_j, B = r_{j+1}.

Proof. Let (u_{i-1}, u_i)^T = Q_i \cdots Q_j (A, B)^T, so that u_j = A, u_{j+1} = B and u_{i-1} = q_i u_i + u_{i+1}. By downward induction, \deg u_i > \deg u_{i+1} for all i: this holds for i = j, and then \deg u_{i-1} = \deg q_i + \deg u_i > \deg u_i. Hence each u_{i+1} is the remainder and q_i the quotient of the division of u_{i-1} by u_i, which are unique. By induction from i = 1, the u_i are the Euclidean remainders.

The half-gcd#

Definition 2. Let \deg a = d > \deg b. The half-gcd of (a, b) is the triple (M_j, r_j, r_{j+1}) for the unique j with

\deg r_j \ge \lceil d/2 \rceil > \deg r_{j+1},

that is 2 \deg r_j \ge d > 2 \deg r_{j+1}.

The index j exists and is unique because the degrees of the remainders decrease strictly from \deg r_0 = d (and \deg 0 = -\infty). The entries of M_j have degree at most d - \deg r_j \le \lfloor d/2 \rfloor.

The key fact is that the half-gcd of the top halves of a and b is an initial part of the half-gcd of a and b.

Lemma 3 (Lifting). Let \deg a = d > \deg b \ge m, write a = x^m a_1 + a_0 and b = x^m b_1 + b_0 with \deg a_0, \deg b_0 < m, and let d' = d - m and (M, A^*, B^*) be the half-gcd of (a_1, b_1). Then (A, B)^T = M^{-1} (a, b)^T satisfies

(A, B) = x^m (A^*, B^*) + M^{-1} (a_0, b_0), \qquad \deg A = m + \deg A^* \ge m + \lceil d'/2 \rceil > \deg B,

and A, B are consecutive remainders of (a, b), with (a, b)^T = M (A, B)^T.

Proof. The first identity is linearity. M^{-1} = \pm \operatorname{adj}(M) has entries of degree at most d' - \deg A^* \le \lfloor d'/2 \rfloor \le \deg A^*, so M^{-1}(a_0, b_0) has degree less than m + \deg A^*. Hence \deg A = m + \deg A^*. Also \deg B^* < \lceil d'/2 \rceil \le \deg A^*, so \deg B < m + \deg A^* = \deg A. Since M is a product of quotient matrices with quotients of degree at least 1 (those of (a_1, b_1), where \deg a_1 = d' > \deg b_1), Lemma 1 applies.

Algorithm 1: Half-gcd

Input: a, b with \deg a = d > \deg b.

Output: the half-gcd (M, A, B) of (a, b).

  1. m \gets \lceil d/2 \rceil. If \deg b < m, return (I, a, b).
  2. If d is below a cutoff, run the Euclidean algorithm, multiplying up the quotient matrices, until the remainder degree drops below m, and return.
  3. (R, \cdot, \cdot) \gets half-gcd of (a \operatorname{div} x^m, b \operatorname{div} x^m); (A, B) \gets R^{-1}(a, b), computed as in Lemma 3 from the outputs of the recursive call and the low parts a_0, b_0.
  4. If \deg B < m, return (R, A, B).
  5. Divide: A = qB + D. Let k = 2m - \deg B.
  6. (S, \cdot, \cdot) \gets half-gcd of (B \operatorname{div} x^k, D \operatorname{div} x^k); (A', B') \gets S^{-1}(B, D), again via Lemma 3.
  7. Return (R \, Q \, S, A', B') with Q = \begin{pmatrix} q & 1 \\ 1 & 0 \end{pmatrix}.

Theorem 4. Algorithm 1 returns the half-gcd of (a, b).

Proof. By induction on d. Step 2 is the definition. In step 3, with d' = d - m = \lfloor d/2 \rfloor, Lemma 3 shows that A, B are consecutive remainders with (a, b)^T = R (A, B)^T and \deg A \ge m + \lceil d'/2 \rceil \ge m. If \deg B < m, they straddle m = \lceil d/2 \rceil as required.

Otherwise m \le \deg B < m + \lceil d'/2 \rceil, so 0 < k \le m, and B, D are the next pair of remainders. The polynomial B \operatorname{div} x^k has degree d'' = \deg B - k = 2 \deg B - 2m, and \lceil d''/2 \rceil = \deg B - m. Applying Lemma 3 to (B, D) with k in place of m (if \deg D < k the recursive call returns the identity and (A', B') = (B, D)), A' and B' are consecutive remainders of (B, D), hence of (a, b), with

\deg A' \ge k + \lceil d''/2 \rceil = m > \deg B'.

The matrix is RQS because (a, b)^T = R(A, B)^T = RQ(B, D)^T = RQS(A', B')^T.

Proposition 5 (Complexity). Algorithm 1 uses O(M(d) \log d) operations in K.

Proof. Both recursive calls receive inputs of degree at most about d/2 (d' = \lfloor d/2 \rfloor, and d'' < 2 \lceil d'/2 \rceil). The other steps are a constant number of products of polynomials of degree O(d) (applying R^{-1} and S^{-1}, whose entries have degree at most d/4 or so, to vectors of degree O(d), one division, and one 2 \times 2 matrix product), which cost O(M(d)). Hence H(d) \le 2H(d/2) + O(M(d)), which gives O(M(d) \log d) since M(d)/d is nondecreasing.

GCD, extended GCD and resultant#

GCD. Starting from (a, b), one half-gcd step produces consecutive remainders (A, B) with \deg B < d/2. After one division the degree has halved, and repeating gives the GCD. The total cost is dominated by the first call, so it is also O(M(d) \log d). Only the remainders are needed, not the matrix, so the top-level call can skip the final matrix product.

Extended GCD. The cofactors s, t with sa + tb = \gcd(a, b) are read off from the product of all the matrices: if (a, b)^T = M (r_k, 0)^T, then M^{-1} = \pm \operatorname{adj}(M) gives r_k as a combination of a and b.

Resultant. For a = qb + r with \deg a = \alpha, \deg b = \beta, \deg r = \rho and leading coefficient c of b,

Lemma 6. \operatorname{Res}(a, b) = (-1)^{\alpha\beta} c^{\alpha - \rho} \operatorname{Res}(b, r) if r \ne 0, \operatorname{Res}(a, b) = c^\alpha if \beta = 0, and \operatorname{Res}(a, b) = 0 if r = 0 and \beta > 0.

Proof. \operatorname{Res}(b, a) = c^\alpha \prod_{b(\xi) = 0} a(\xi) over the roots \xi of b in an algebraic closure, with multiplicity. Since a(\xi) = r(\xi), this equals c^{\alpha - \rho} \operatorname{Res}(b, r), and \operatorname{Res}(a, b) = (-1)^{\alpha\beta} \operatorname{Res}(b, a). The other cases follow from the same product formula.

So the resultant is the product, over the steps of the Euclidean algorithm, of factors that depend only on the degrees and leading coefficients of the remainders [GG2003, Ch. 6 and §11.2]. The half-gcd algorithm performs the same divisions as the Euclidean algorithm, except that in the recursive calls the remainders are truncated: they have the correct leading coefficients, and their degrees are offset by the number of dropped coefficients. FLINT accumulates the factors of Lemma 6 during the half-gcd, tracking that offset. The last remainder produced inside a recursive call is not yet a true remainder (its low coefficients are corrected only by the lifting step of Lemma 3), so its factor is deferred until its true degree is known. Fast computation of resultants and subresultants along these lines is due to Reischert [Rei1997]; see also [GG2003, §11.2].

Use in FLINT#

The algorithm is implemented once, in gr_poly (_gr_poly_hgcd, _gr_poly_gcd_hgcd, _gr_poly_xgcd_hgcd, _gr_poly_resultant_hgcd), and nmod_poly calls it with the gr context for \mathbb{Z}/n\mathbb{Z} and cutoffs tabulated by the bit length \beta of the modulus:

CutoffFunctionRange of valuesMeaning
innernmod_poly_hgcd_iter_recursive_cutoff60–244step 2 of Algorithm 1: Euclid inside the recursion
outernmod_poly_hgcd_outer_cutoff310–1585in the GCD and resultant loops, finish with Euclid below this length
GCDnmod_poly_gcd_hgcd_cutoff470–2199nmod_poly_gcd and nmod_poly_resultant use the half-gcd from this length of the shorter input
extended GCDnmod_poly_xgcd_hgcd_cutoff72–388nmod_poly_xgcd uses the half-gcd from this length

The extended GCD switches earlier than the GCD, since the Euclidean algorithm must then also update the cofactors, while the half-gcd computes the transformation matrices in any case.

Below the GCD cutoff, nmod_poly_gcd and nmod_poly_resultant use the Euclidean algorithm, in Montgomery representation (nmod_poly_gcd_euclidean_redc_fast) when \ell \ge 32 and \mathtt{FLINT\_BITS}/2 \le \beta \le \mathtt{FLINT\_BITS} - 2; see Montgomery multiplication (REDC). nmod_poly_resultant first tries special methods for very short inputs.

The 2 \times 2 polynomial matrix products of step 7 use eight multiplications when an entry is shorter than 20, and Strassen's seven-multiplication scheme otherwise. With fft_small, when all entries have length at least 128, each of the eight entries is transformed once and the four products are accumulated in the transformed domain, so that 8 forward and 4 inverse transforms replace the 24 transforms of eight separate FFT products; from length 1024 the transforms are distributed over the thread pool.

The functions:

The functions are described, with examples, in Greatest common divisors and Resultants and discriminants.

slong nmod_poly_hgcd(nmod_poly_t m11, nmod_poly_t m12, nmod_poly_t m21, nmod_poly_t m22, nmod_poly_t A, nmod_poly_t B, const nmod_poly_t a, const nmod_poly_t b)
void nmod_poly_gcd_hgcd(nmod_poly_t G, const nmod_poly_t A, const nmod_poly_t B)
void nmod_poly_xgcd_hgcd(nmod_poly_t G, nmod_poly_t S, nmod_poly_t T, const nmod_poly_t A, const nmod_poly_t B)
ulong nmod_poly_resultant_hgcd(const nmod_poly_t A, const nmod_poly_t B)

Referenced from#

References#

  • [Knu1970] D. E. Knuth, "The analysis of algorithms", Actes du Congrès International des Mathématiciens (Nice, 1970), tome 3, Gauthier-Villars (1971), 269–274
  • [Sch1971] A. Schönhage, "Schnelle Berechnung von Kettenbruchentwicklungen", Acta Informatica 1 (1971), 139–144
  • [Moe1973] R. T. Moenck, "Fast computation of GCDs", Proceedings of the 5th Annual ACM Symposium on Theory of Computing (STOC 1973), 142–151
  • [BGY1980] R. P. Brent, F. G. Gustavson and D. Y. Y. Yun, "Fast solution of Toeplitz systems of equations and computation of Padé approximants", Journal of Algorithms 1:3 (1980), 259–295
  • [ThullYap1990] Thull, K. and Yap, C. : A Unified Approach to HGCD Algorithms for Polynomials and Integers, (1990)
  • [GG2003] J. von zur Gathen and J. Gerhard, Modern Computer Algebra, second edition, Cambridge University Press (2003)
  • [Tho2002] E. Thomé, "Subquadratic computation of vector generating polynomials and improvement of the block Wiedemann algorithm", Journal of Symbolic Computation 33:5 (2002), 757–775
  • [Rei1997] D. Reischert, "Asymptotically fast computation of subresultants", ISSAC 1997, 233–240
Cite this page

The FLINT team, “Half-GCD,” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/algorithms/polynomials/half-gcd.html (generated 2026-10-20).