Polynomial GCDs and resultants

A greatest common divisor of two univariate polynomials A and B can be computed over any field by the Euclidean algorithm, and over any unique factorisation domain by a pseudo-remainder sequence, which avoids divisions. These algorithms are implemented once in gr_poly for any coefficient ring. Which one is fast depends on the ring: the Euclidean algorithm and its asymptotically fast form, the half-gcd algorithm, are the methods of choice over finite fields; over \mathbb{Q} and \mathbb{Z} their intermediate coefficients grow, and the best algorithms reduce the problem modulo primes instead.

Throughout, A and B are nonzero with m = \deg A \ge n = \deg B, and \operatorname{lc} denotes the leading coefficient.

The Euclidean algorithm#

Over a field K, the Euclidean remainder sequence is R_0 = A, R_1 = B and R_{i+1} = R_{i-1} \bmod R_i, until R_{k+1} = 0; then R_k is a GCD of A and B, and \gcd(A, B) is its monic associate. Each step costs O(\deg R_i \cdot (\deg R_{i-1} - \deg R_i + 1)) operations, and the whole sequence O(mn) operations [GG2003, Ch. 3]. Keeping track of cofactors S_i, T_i with S_i A + T_i B = R_i (starting from S_0 = 1, T_0 = 0, S_1 = 0, T_1 = 1 and applying the quotients) gives the extended GCD at the same asymptotic cost.

The resultant follows from the remainder sequence. If R = A \bmod B has degree k, then

while \operatorname{res}(A, B) = 0 if R = 0 and n > 0, and \operatorname{res}(A, B) = \operatorname{lc}(B)^m if n = 0. The Euclidean resultant algorithm accumulates these factors along the sequence (Resultants).

In gr_poly these are _gr_poly_gcd_euclidean, _gr_poly_xgcd_euclidean and _gr_poly_resultant_euclidean. They need to invert the leading coefficient of each remainder and to decide when a remainder vanishes. Over a field with decidable equality this always works. Over \mathbb{Z}/n\mathbb{Z} with composite n pretending to be a field (gr_ctx_set_is_pretend_field), a non-invertible leading coefficient makes the algorithm fail with an error status, and the factor of n it has found can be read with gr_ctx_recover_zero_divisor, which is the basis of dynamic evaluation [DDDD1985]. Over inexact rings such as real balls the algorithm succeeds when the zero tests are decided, typically for exactly represented inputs, and returns GR_UNABLE when a remainder cannot be shown to be zero or nonzero: the GCD of polynomials known only approximately is not a well-posed problem.

The half-gcd algorithm#

The quotients of the first half of the Euclidean remainder sequence depend only on the top halves of A and B. The half-gcd algorithm computes the 2 \times 2 polynomial matrix that maps (A, B) to a pair of consecutive remainders of about half the degree, recursively, and repeats until the remaining polynomials are short, giving the GCD, the extended GCD and the resultant with O(M(n) \log n) operations [Moe1973; BGY1980; ThullYap1990]. The generic implementation _gr_poly_hgcd uses the Euclidean algorithm below an inner cutoff, and _gr_poly_gcd_hgcd, _gr_poly_xgcd_hgcd and _gr_poly_resultant_hgcd call it repeatedly above an outer cutoff. Its matrix products are done by the ring method _gr_poly_hgcd_mat_mul, with Strassen's seven products and, when the ring has an FFT representation, products of transformed entries (Greatest common divisors). Over infinite fields such as \mathbb{Q} the half-gcd algorithm is not used by default, because its intermediate coefficients are no smaller than those of the Euclidean algorithm while its overhead is larger, and over approximate fields because of numerical instability.

Pseudo-remainder sequences#

Over an integral domain D that is not a field, A \bmod B is not defined, but the pseudo-remainder

has coefficients in D. Repeated pseudo-division leads to a sequence whose remainders are associate (over the fraction field) to the Euclidean remainders, but whose coefficients grow exponentially in the number of steps. The subresultant pseudo-remainder sequence of Collins and of Brown and Traub [Col1967; Bro1971; BT1971] divides each pseudo-remainder by a factor known in advance to divide it exactly, so that its elements are, up to sign, the subresultants of A and B: determinants of submatrices of the Sylvester matrix, whose size is bounded by Hadamard's inequality and grows only linearly in the step number.

The generic _gr_poly_gcd_subresultant follows Cohen's formulation [Coh1996, Algorithm 3.3.1] over a unique factorisation domain:

1. Compute the contents a = \operatorname{cont}(A) and b = \operatorname{cont}(B) (GCDs of the coefficients, using the ring's GCD) and d = \gcd(a, b); replace A and B by their primitive parts, and set g = h = 1.
2. Let \delta = \deg A - \deg B and R = \operatorname{prem}(A, B). If R = 0, stop with B; if \deg R = 0, the GCD is d.
3. Set A \leftarrow B, B \leftarrow R / (g h^\delta) (an exact division), g \leftarrow \operatorname{lc}(A) and h \leftarrow h^{1 - \delta} g^\delta (also exact), and go to step 2.
4. The GCD is d times the primitive part of the final B.

The extended version _gr_poly_xgcd_subresultant and the resultant version _gr_poly_resultant_subresultant follow the same sequence. The resultant version needs only exact divisions and works over any integral domain with exact division.

As a last resort, _gr_poly_resultant_sylvester evaluates the determinant of the Sylvester matrix with gr_mat_det, which uses a division-free algorithm by default and therefore works over any commutative ring, including rings with zero divisors.

Generic algorithm selection#

gr_poly_gcd, gr_poly_xgcd and gr_poly_resultant call ring methods, whose defaults choose from the properties of the ring:

For the resultant, each failing algorithm falls back to the next one in the list subresultant, Sylvester determinant. The "small formulas" (_gr_poly_resultant_small) are division-free closed forms for a divisor of length at most 3 or a dividend of length at most 2. The non-underscore functions canonicalise the GCD: monic over a field, with a canonical unit as leading coefficient over other rings (positive over \mathbb{Z}); see Greatest common divisors over generic rings for the conventions and an example.

Specialised algorithms#

Several rings override the methods, mostly to tune the cutoffs between the Euclidean and half-gcd algorithms, which depend on the cost of the ring's multiplication:

gr_poly_gcd over fmpz calls _fmpz_poly_gcd through the ring method, but over fmpq it runs the generic Euclidean algorithm, which is much slower than fmpq_poly_gcd for large inputs. Resultants are better served: polynomial rings over \mathbb{Z}/p\mathbb{Z}, \mathbb{Z} and \mathbb{Q} override the resultant method with bivariate algorithms by evaluation and interpolation and by Chinese remaindering (Bivariate resultants by evaluation and interpolation, Resultants and discriminants over Z).