Multivariate polynomials: representations and algorithms

FLINT represents multivariate polynomials sparsely, as sorted lists of terms with packed exponent vectors, in one module per coefficient ring. The algorithms are shared between the modules as far as possible: the monomial code is in mpoly.h, and most GCD and factoring algorithms exist in a version over \mathbb{Z} and versions over finite fields, with the rational and rational function modules built on top of the integer one.

Coefficients	Polynomials	Factorisation	Fractions
any `gr` ring	`gr_mpoly`
\mathbb{Z}	fmpz_mpoly.h	`fmpz_mpoly_factor`	fmpz_mpoly_q.h
\mathbb{Q}	fmpq_mpoly.h (content times `fmpz_mpoly`)	`fmpq_mpoly_factor`
\mathbb{Z}/n\mathbb{Z}, word-size n	`nmod_mpoly`	`nmod_mpoly_factor`
\mathbb{Z}/n\mathbb{Z}, any n	`fmpz_mod_mpoly`	`fmpz_mod_mpoly_factor`	`fmpz_mod_mpoly_q`
\mathbb{F}_q	`fq_nmod_mpoly`, `fq_zech_mpoly`	`fq_nmod_mpoly_factor`, `fq_zech_mpoly_factor`

The generic gr_mpoly implements arithmetic (heap multiplication and division, threaded versions) for any coefficient ring; the gr contexts for the specialised types (gr_ctx_init_fmpz_mpoly and others) call the specialised functions, which are faster and also provide GCDs and factorisation.

Representation#

AlgorithmMonomial orderings and exponent packingHow FLINT stores sparse multivariate polynomials: terms sorted by a monomial ordering, exponent vectors packed into machine words so that comparing, multiplying and dividing monomials are word operations, with guard bits for overflow detection; and the deflations, repackings and changes of variables that the algorithms apply to this representation.

A polynomial is stored with its terms sorted for one of the orderings ORD_LEX, ORD_DEGLEX, ORD_DEGREVLEX, and the exponent vectors packed into words with a guard bit per field, so that comparing, multiplying and testing divisibility of monomials are word operations. The packing width grows automatically. The ordering affects the speed of some operations (the dense array methods differ for lex and degree orderings) but not the results of GCDs and factorisations, which are computed in lex order internally and converted back.

Arithmetic#

AlgorithmSparse and dense multiplicationMultiplication of sparse multivariate polynomials with Johnson's heap algorithm and the improvements of Monagan and Pearce (packed monomials, chaining, small heaps, three-word accumulation, parallel splitting of the output), dense array accumulation, Kronecker substitution for dense inputs, and powering by the FPS recurrence.

Operation	Method	Used when	Functions
multiplication	Johnson's heap with chaining	sparse inputs, any ring	`fmpz_mpoly_mul_johnson`, `gr_mpoly_mul_heap`
	parallel heap	large inputs, several threads	`fmpz_mpoly_mul_heap_threaded`
	dense array accumulation	product dense in a box of moderate size	`fmpz_mpoly_mul_array`
	Kronecker substitution	dense inputs	`fmpz_mpoly_mul_dense`
powering	FPS recurrence	\mathbb{Z}, exponent \ge 3	`fmpz_mpoly_pow_fps`
	repeated multiplication	finite fields	`nmod_mpoly_pow_rmul`

fmpz_mpoly_mul and nmod_mpoly_mul choose automatically; the rules are described in Multiplication.

AlgorithmDivision, exact division and square rootsMultivariate division with remainder with respect to a monomial ordering, and its heap implementation by Monagan and Pearce, in which a heap of size at most the number of terms of the divisor supplies the terms of the current remainder; early termination of exact division, quasi-division over the integers, division by several polynomials, parallel exact division, and the analogous heap algorithm for square roots with its non-square certificates.

Operation	Method	Functions
exact division	heap division with early termination; threads for large dividends	`fmpz_mpoly_divides`
division with remainder	heap division (floor division of coefficients over \mathbb{Z})	`fmpz_mpoly_divrem`
pseudo-division	heap division with scaling	`fmpz_mpoly_quasidivrem`
division by several polynomials	heap	`fmpz_mpoly_divrem_ideal`
square root	heap, with non-square certificates	`fmpz_mpoly_sqrt`

GCDs and resultants#

AlgorithmGreatest common divisorsHow FLINT computes GCDs of multivariate polynomials over the integers: reductions that make most inputs easy (monomial content, deflation, missing variables, univariate cases), degree bounds from random images, and the four modular algorithms it chooses between (Brown's dense algorithm, Zippel's sparse algorithm, the Berlekamp–Massey based algorithm of Hu and Monagan, and EEZ-GCD Hensel lifting), with their correctness proofs and the cost model behind the choice.

fmpz_mpoly_gcd first reduces the problem (monomial contents, deflation, missing variables, univariate cases, divisibility tests) and then chooses between Brown's dense modular algorithm (two variables, dense inputs), the Berlekamp–Massey based sparse algorithm of Hu and Monagan (sparse inputs in three or more variables) and Zippel's sparse algorithm, using degree and term count estimates from random images. Hensel lifting (EEZ-GCD) and subresultants are available on request. Over finite fields, nmod_mpoly_gcd also uses Hensel lifting, and computes over extension fields when the field is too small:

AlgorithmGCDs over small finite fieldsHow the multivariate GCD algorithms over \mathbb{F}_p and \mathbb{F}_q cope with fields too small for random evaluation points, by computing images over extension fields \mathbb{F}_{p^k} or modulo irreducible polynomials of increasing degree; Brown's algorithm over \mathbb{F}_p with pairs of points and early termination; and how nmod_mpoly_gcd chooses among Brown's, Zippel's and Hensel-lifting algorithms.AlgorithmSparse interpolationRecovering a sparse polynomial from its values at a few points: Zippel's method for a known support (transposed Vandermonde systems solved in quadratic time, with extra unknowns for scaled images), and the method of Ben-Or and Tiwari for an unknown support (Berlekamp–Massey, roots of the recurrence polynomial, discrete logarithms modulo a smooth prime after Kronecker substitution) with the Berlekamp–Massey algorithm computed incrementally by half-gcd steps, as used by FLINT's sparse GCD and factoring algorithms.

Resultants and discriminants (fmpz_mpoly_resultant, fmpz_mpoly_discriminant) use the subresultant algorithm of Ducos [Duc2000] on univariate polynomials with multivariate coefficients, except that bivariate dense inputs go through a multimodular evaluation and interpolation method:

AlgorithmBivariate resultants by evaluation and interpolationThe resultant \mathrm{Res}_y(A, B) \in \mathbb{F}_p[x] of two bivariate polynomials, computed from univariate resultants at many values of x (a geometric progression or roots of unity) and interpolation, with the degree bound, the correction for specialisations that lower the degree, and FLINT's cutoffs.

Factorisation#

AlgorithmFactorisationHow FLINT factors multivariate polynomials over the integers: removal of contents and powers of variables, squarefree decomposition by Yun's algorithm, compression of the exponents by unimodular changes of variables, irreducibility certificates from Newton polytopes, direct formulas for low degree, dense bivariate factorisation with Hensel lifting and recombination, degree pruning with univariate images, and the three lifting strategies of Wang, Zippel and Zassenhaus.AlgorithmMultivariate Hensel liftingLifting a factorisation of A(x_0, \alpha_1, \ldots, \alpha_m) to a factorisation of A one variable at a time: uniqueness given the leading coefficients, the multivariate partial fraction (diophantine) equations solved by Taylor expansion, the quartic and quintic error updates of FLINT's fmpz_mpoly_hlift, bivariate lifting modulo prime powers, sparse lifting, and the leading coefficient methods of Wang and Kaltofen.

Stage	Method	Functions
content	GCDs of coefficients with respect to each variable	`fmpz_mpoly_factor_content`
squarefree	Yun's algorithm (with p-th roots in characteristic p)	`fmpz_mpoly_factor_squarefree`
simplification	unimodular compression of exponents; Newton polytope test	`mpoly_compression_set`
low degree	linear and quadratic in the main variable	`fmpz_mpoly_factor`
bivariate dense	Hensel lifting and recombination	`fmpz_bpoly_factor`
general	evaluation, degree pruning, then Wang (dense lifting), Zippel (sparse lifting) or Zassenhaus (lifting with recombination)	`fmpz_mpoly_factor_wang`, `fmpz_mpoly_factor_zippel`, `fmpz_mpoly_factor_zassenhaus`

Factoring over finite fields follows the same stages, with the changes described in:

AlgorithmFactoring over finite fieldsThe multivariate factorisation pipeline of nmod_mpoly_factor and fq_nmod_mpoly_factor, with what is specific to characteristic p: separable factorisation with p-th roots, deflation by powers of p, factoring over extension fields with recombination of Frobenius orbits when the field is too small for good evaluation points, quadratic equations in characteristic 2, bivariate factorisation with recombination by linear algebra, and Wang's leading coefficients along random curves.

Practical advice#

Use the specialised modules rather than gr_mpoly over \mathbb{Z}, \mathbb{Q} and finite fields: they accumulate coefficients without intermediate reductions, use dense methods when they pay, and provide GCDs and factorisation.
The ordering only matters for the leading terms (and hence for division with remainder and Gröbner bases) and for the speed of the array methods; lex order is the natural one for recursive computations.
Threads (flint_set_num_threads) are used by multiplication and exact division of large polynomials, and by Brown's GCD algorithm in nmod_mpoly_gcd.
GCDs and factorisations are randomised but never return wrong results: every algorithm verifies its output, by exact division or by a proof of its stopping criterion. Their running time can vary with the random choices.