Greatest common divisors

$\mathbb{Z}[x_0, \ldots, x_{n-1}]$ is a unique factorisation domain, so two polynomials $A$, $B$ have a greatest common divisor $G$, unique up to sign; FLINT returns the one with positive leading coefficient, and also the cofactors $\bar A = A/G$ and $\bar B = B/G$ if requested. Euclid's algorithm with pseudo-remainders works in principle, but the coefficients of the intermediate remainders grow and their supports fill up even when $G$ and the cofactors are small and sparse. All competitive algorithms are therefore modular: they compute images of $G$ modulo primes and at evaluation points, where the problem is univariate or bivariate over a finite field, and reconstruct $G$ from these images [Bro1971; GCL1992, Ch. 7; GG2003, §6]. The difficulty is to know the shape of $G$ in advance or to detect the images that do not reflect it.

The algorithms below are those used for fmpz_mpoly (and through it fmpq_mpoly); the function-level description, with an example, is Greatest common divisors. The algorithms over finite fields are the same inside, with the additional problem that a small field may not have enough evaluation points, which is treated in GCDs over small finite fields. The reconstruction of sparse polynomials from values is described in Sparse interpolation and the lifting of factorisations on Multivariate Hensel lifting. There is no generic multivariate GCD in gr_mpoly; the gr contexts of these modules call the specialised functions.

Contents and normalisation#

Write $R = \mathbb{Z}[x_1, \ldots, x_{n-1}]$, so that $A \in R[x_0]$. The content $\operatorname{cont}_{x_0}(A) \in R$ is the GCD of the coefficients of $A$ as a polynomial in $x_0$, and $A$ is primitive if this content is 1. By Gauss's lemma products of primitive polynomials are primitive, so

and the algorithms below assume primitive inputs, computing contents by recursive GCDs in fewer variables (fmpz_mpolyl_content). For an integer content this is the usual $\gcd(\operatorname{cont}(A), \operatorname{cont}(B))$ in $\mathbb{Z}$.

The images of $\tilde G$ modulo a prime or at an evaluation point are determined up to the normalisation of the leading coefficient, which is known: this is how images computed independently are made consistent before they are combined. At the end, $G$ is the primitive part of $\tilde G$.

The same reasoning applies with several main variables, which gives a common frame for the algorithms below [Sch-mpoly]. View $A$ and $B$ in $R'[x_{m+1}, \ldots, x_{n-1}]$ with $R' = D[x_0, \ldots, x_m]$, make them primitive over $D[x_{m+1}, \ldots]$, compute GCDs in $m + 1$ variables of images at evaluation points of the other variables, impose the image of $\gamma = \gcd(\operatorname{lc}(A), \operatorname{lc}(B))$ (leading coefficients with respect to $x_0, \ldots, x_m$) as their leading coefficient, and interpolate. If the primitive part of the interpolated polynomial divides $A$ and $B$, it is the GCD. Brown's algorithm takes $m = n - 2$, with dense interpolation in the last variable and recursion down to the Euclidean algorithm; Zippel's algorithm takes $m = 0$, and the algorithm of Hu and Monagan $m = 1$, both with sparse interpolation.

Easy cases#

Before any modular computation, fmpz_mpoly_gcd removes as much structure as possible; the list of steps is described in Easy cases. They rest on three facts.

• Monomials: the GCD of the monomial contents (term-wise minima of the exponents) can be split off, and $\gcd(c x^e, B)$ is $\gcd(c, \operatorname{cont}(B)) \cdot \min(x^e, \operatorname{mcont}(B))$.
• Deflation: by Lemma 4 of Monomial orderings and exponent packing, exponents can be divided by common strides without changing the GCD, which reduces degrees, and inputs can then be repacked into one word.
• Missing variables: if $x_v$ occurs in $A$ but not in $B$, then $\gcd(A, B) = \gcd(\operatorname{cont}_{x_v} A, B)$, because every divisor of $B$ is free of $x_v$ and divides $A$ if and only if it divides all coefficients of $A$ with respect to $x_v$.

With one common variable the GCD is univariate (fmpz_poly_gcd, see GCDs of polynomials over Z), and when one input is a monomial times a univariate polynomial of degree at most 100 in some variable, the GCD is computed with the coefficients of the other input with respect to the remaining variables.

Degree bounds#

mpoly_gcd_info_set_estimates uses this with a random 63-bit prime and random evaluation points, for each variable in turn: the degrees of the univariate GCDs are upper bounds for the degrees of $G$ (exactly the degrees with high probability), and the numbers of their nonzero terms estimate the term counts of $G$. If all bounds are zero, $G$ is the GCD of the contents; if they equal the degrees of $A$ (or $B$), the code tests whether $A \mid B$ (or $B \mid A$) by exact division. The estimates then feed the cost model of Choosing an algorithm.

Brown's algorithm#

Brown's algorithm [Bro1971] computes images of $\tilde G$, $\bar A$ and $\bar B$ modulo primes $p$ and reconstructs them by the Chinese remainder theorem. Modulo $p$ it recurses on the variables: the last variable is evaluated at points of $\mathbb{F}_p$, the GCDs of the images in one variable fewer are computed recursively, normalised to the image of the leading coefficient $\gamma$, and interpolated (nmod_mpolyn_gcd_brown_smprime, described in GCDs over small finite fields). Over $\mathbb{Z}$:

The test of step 6 replaces trial division, which is the expensive final step of most modular GCD algorithms. Brown's algorithm costs about as much as a few dense multiplications of the inputs: it is the method of choice for dense inputs and for two variables, and hopeless for sparse inputs of high degree in many variables, where the dense box has $\prod_i (d_i + 1)$ entries.

Zippel's algorithm#

Zippel's algorithm [Zip1979] exploits sparsity through the shape assumption: the image of $G$ modulo a random prime has the same monomials as $G$. The first prime gives the monomials; for the others, only the coefficients are unknown, and they are found from few evaluations by solving linear systems (sparse interpolation, Known support: transposed Vandermonde systems).

The result is correct because of the final divisions and the degree bound: a primitive common divisor of $A$ and $B$ whose degree in $x_0$ is the degree of the GCD of images, which by Lemma 2 is at least $\deg_{x_0} G$, is $G$ itself. The leading coefficient of $G$ in $x_0$ is in general a polynomial in $x_1, \ldots$, not a monomial, so the univariate images computed during sparse interpolation are only known up to scalar multiples; FLINT treats these scalars as additional unknowns of the linear systems, as proposed by de Kleine, Monagan and Wittkopf [dKMW2005].

The bound of step 3, in bits, is the bit size of the smaller of the bounds of fmpz_mpoly_factor_bound_si for the divisors of $A$ and of $B$, which use the degrees in all variables, plus the bit size of $\gamma$.

Sparse interpolation with Berlekamp–Massey#

Zippel's algorithm needs one recursive sparse interpolation per variable. Hu and Monagan [HuMon2016] instead treat two main variables $x_0, x_1$ densely and all others at once: the coefficients of $G$ as a polynomial in $x_0, x_1$ are sparse polynomials in $x_2, \ldots$, which are recovered by the algorithm of Ben-Or and Tiwari [BenOrTiw1988] from their values at the powers of one point (Unknown support: Ben-Or–Tiwari).

The number of evaluations is about twice the number of terms of the largest coefficient of $\tilde G$ (or of the cofactor used), against roughly the number of terms times the number of variables for Zippel's algorithm. Using $\Gamma$ to normalise the images solves the leading coefficient problem at the price of possibly introducing a factor $\Gamma / \operatorname{lc}(G)$, removed at the end.

Hensel lifting#

The EEZ-GCD algorithm [MosYun1973; Wan1980] views the GCD problem as a factorisation problem: if $g = \gcd(A(x_0, \alpha), B(x_0, \alpha))$ for an evaluation point $\alpha \in \mathbb{Z}^{n-1}$ and $g$ is coprime to the cofactor $\bar a = A(x_0, \alpha)/g$, then $A = G \bar A$ is the unique lifting of $A(x_0, \alpha) = g \cdot \bar a$ with the correct leading coefficients (Lemma 1 of Multivariate Hensel lifting).

When the GCD has few terms but a dense structure in few variables, Hensel lifting is efficient, and it needs only one good evaluation point, which matters over small finite fields (GCDs over small finite fields). Its weakness is the cost of dense lifting for sparse polynomials in many variables. fmpz_mpoly_gcd does not use it; it is available as fmpz_mpoly_gcd_hensel. A last method, fmpz_mpoly_gcd_subresultant, computes the last nonzero subresultant in a well chosen main variable with Ducos' algorithm [Duc2000] and is useful only for very small inputs.

Choosing an algorithm#

After the easy cases, the inputs are converted to a lex context in the variables that occur, ordered differently for each algorithm (the variables of largest degree become the dense variables of Brown's algorithm and the main variables of the others), and the estimates of Degree bounds give costs:

• for Brown's algorithm, roughly the size of the dense boxes of $A$, $B$, $G$ and the cofactors, adjusted by the densities $\ell_A/\prod_i (\deg_{x_i} A + 1)$ and the expected number of points before stabilisation (mpoly_gcd_info_measure_brown);
• for the Berlekamp–Massey algorithm, the number of evaluations (the maximal number of terms of a coefficient of $G$ in the two main variables) times the cost of a bivariate GCD and of evaluating $A$ and $B$ (mpoly_gcd_info_measure_bma).

With at most two variables Brown's algorithm is used. With three or more, the Berlekamp–Massey algorithm is tried first if its estimated time is less than Brown's and either the inputs are sparse ($n(\delta_A + \delta_B) < 1$ for the densities $\delta$) or its estimate is below 1% of Brown's; otherwise Brown's algorithm comes first. If both fail, Zippel's algorithm is tried. Each algorithm can fail (for instance when it runs out of good points, or when the exponents do not fit), never returning a wrong result, and the next one is used.

Use in FLINT#

• Driver: fmpz_mpoly_gcd, fmpz_mpoly_gcd_cofactors and the forced variants through _fmpz_mpoly_gcd_algo, with the bookkeeping of mpoly_gcd_info_t (mpoly_gcd_info_limits, mpoly_gcd_info_stride, mpoly_gcd_info_set_perm, mpoly_gcd_info_measure_brown, mpoly_gcd_info_measure_bma, mpoly_gcd_info_measure_zippel, mpoly_gcd_info_measure_hensel).
• Brown: fmpz_mpolyl_gcd_brown (Algorithm 1, Theorem 3), with nmod_mpolyn_gcd_brown_smprime for the images and the interpolation helpers fmpz_mpoly_interp_reduce_p_mpolyn, fmpz_mpoly_interp_crt_p_mpolyn; a threaded version over several primes is fmpz_mpolyl_gcd_brown_threaded_pool.
• Zippel: fmpz_mpolyl_gcd_zippel (Algorithm 2), with nmod_mpolyl_gcdp_zippel_smprime and nmod_mpolyl_gcds_zippel.
• Berlekamp–Massey: fmpz_mpolyl_gcd_zippel2 (Algorithm 3).
• Hensel: fmpz_mpolyl_gcd_hensel (Algorithm 4), with fmpz_mpoly_hlift.
• Subresultants: fmpz_mpoly_gcd_subresultant.
• fmpq_mpoly_gcd and friends apply these to the integer parts; the finite field versions are described in GCDs over small finite fields.