Sparse and dense multiplication

Let B and C be polynomials with \ell_B and \ell_C terms in the distributed representation of Monomial orderings and exponent packing. Their product has at most \ell_B \ell_C terms, and in the sparse case about that many; in the dense case far fewer, because many products b_i c_j share a monomial. The classical algorithm forms all products and sorts or merges them. Done naively (adding b_i C to an accumulator one row at a time) it costs O(\ell_B^2 \ell_C) monomial comparisons in the worst case. FLINT's algorithms avoid this in different ways: a heap (sparse inputs), a dense array (moderately dense inputs) or a reduction to univariate multiplication (dense inputs).

The heap method#

Johnson's algorithm [Joh1974] merges the \ell_B sorted sequences b_i C = (b_i c_0, b_i c_1, \ldots) with a priority queue holding the next unused product of each sequence. Monagan and Pearce [MonPea2007] made it practical with three refinements: packed monomials (so that the heap keys are words compared with one instruction), chaining (products with equal monomials are linked in one heap node, so the heap only holds distinct monomials) and a rule that keeps the heap small. We describe the algorithm as implemented, with rows indexed by B and columns by C.

In the dense case, chaining makes the cost per product close to constant: many products are popped together and the heap stays much smaller than \ell_B. The packed monomial of each candidate product is a word sum (N = 1) or an N-word sum, and the comparisons use the packed comparison.

Coefficient accumulation#

Over \mathbb{Z}, if all coefficients of B and C are small fmpz (absolute value below 2^{62} on a 64-bit machine), each product b_i c_j is a signed two-word integer, and a sum of at most \ell_B of them has absolute value below \ell_B 2^{124} < 2^{190}; it is accumulated in three words with smul_ppmm and three-word additions and converted to an fmpz once per output term. Otherwise fmpz_addmul is used. Over \mathbb{Z}/n\mathbb{Z} (nmod_mpoly), products are accumulated unreduced in up to three words and reduced once per output term, the delayed reduction of Choosing a method for modular arithmetic. This accounts for much of the speed of the specialised versions compared with the generic one.

Parallel heap multiplication#

Monagan and Pearce [MonPea2009] parallelise by splitting the output into ranges of monomials, each computed by a heap restricted to the products in the range. FLINT (fmpz_mpoly_mul_heap_threaded) cuts the monomials at 4t boundaries for t threads: the i-th boundary is a monomial e such that about ((i+1)/4t)^2 \, \ell_B \ell_C of the products b_k c_l are at most e (found by mpoly_search_monomials, a binary search along the rows for each candidate). The quadratic spacing makes the ranges for the largest monomials, which are processed first, contain more products, so that no thread ends long before the others. For each range and each row of B the first and last columns of C in the range are known, and the restricted heap algorithm is run on these segments. The pieces are concatenated in order at the end.

Dense methods#

Array accumulation#

If the product is dense in a box of moderate size, it is faster to add the products into a dense array indexed by the exponents than to sort them. FLINT's array method (fmpz_mpoly_mul_array) slices the product:

• for ORD_LEX, by the exponent of x_0: with d_i = \deg_{x_i} B + \deg_{x_i} C + 1, the product of the slices of B and C of degrees s and t in x_0 contributes to slice s + t, and each slice is a dense array of \prod_{i \ge 1} d_i \le 300000 entries, with d_0 \le 300;
• for the degree orderings, by total degree: within a slice the exponent of x_{n-1} is determined by the others, and the array is indexed by the exponents of x_0, \ldots, x_{n-2}, each below a bound d on the total degree of the product, so it has d^{n-1} entries.

Each output slice is computed by adding all contributing products of input slices into a zeroed array and then reading off the nonzero entries in decreasing order.

FLINT evaluates this bound in bits for each slice and accumulates in one, two or three words if it fits (with one bit for the sign), and in fmpz otherwise. Slices are distributed over threads in fmpz_mpoly_mul_array_threaded. The method declines (returns 0) when the array would be too large or the exponents do not fit one word.

Kronecker substitution#

For dense inputs of moderate degree the fastest method is to substitute a power of a single variable for each variable and multiply univariate polynomials [Kro1882].

fmpz_mpoly_mul_dense converts B and C to dense arrays in this order (no padding is needed for x_0), multiplies them with fmpz_poly_mul, which uses the univariate algorithms of Polynomial arithmetic over generic rings and its specialisations (Kronecker substitution into integers, multimodular and Schönhage–Strassen methods), and converts back. The cost is that of a univariate product of length D = \prod_i d_i, so the method only pays when D is not much larger than the number of terms of the product. For nmod_mpoly the same is done with nmod_poly_mul.

Squaring and powering#

Squaring with the heap method can use symmetry: B^2 = \sum_i b_i^2 x^{2\beta_i} + 2 \sum_{i < j} b_i b_j x^{\beta_i + \beta_j} needs only about half the products (gr_mpoly_sqr_commutative_heap). For higher powers, repeated multiplication computes B^2, B^3, \ldots with a cost dominated by the last product, and binary powering does not help much for sparse polynomials because the last squaring of a dense intermediate result is as expensive as all the rest. Monagan and Pearce [MonPea2012] instead compute B^k directly by the recurrence for the coefficients of a power of a power series (FPS, after the classical method of Euler; see [Knu1997, §4.7]).

Write B = b_0 x^{\beta_0} + B' and compare the coefficients of a monomial x^{\mu} on both sides of B D(P) - (k-1) P D(B) = 0:

where P = \sum_j p_j x^{\pi_j}. The term with i = 0 is (\phi(\mu - \beta_0) - (k-1)\phi(\beta_0)) b_0 p_{\mu - \beta_0} = (\phi(\mu) - k\phi(\beta_0)) b_0 p_{\mu - \beta_0}. So if \phi(\mu) \ne k \phi(\beta_0), the coefficient of x^{\mu - \beta_0} in P is

in terms of coefficients of P at larger monomials. Generating the products b_i p_j with a heap in decreasing order of \beta_i + \pi_j, the algorithm computes, for each monomial \mu, the coefficient a_\mu = \sum b_i p_j of A = BP and the coefficient p_{\mu - \beta_0} of P at the same time; the latter is needed only when x^{\beta_0} divides x^\mu, which is tested with the packed divisibility test.

Since \ell_P is the number of terms of B^{k-1}, this is about the cost of the final multiplication in repeated multiplication, while avoiding the intermediate products. fmpz_mpoly_pow_ui uses FPS for k \ge 3 (after checking that the output length k(\ell_B - 1) + 1 cannot overflow); the divisions by (k\phi(\beta_0) - \phi(\mu)) b_0 make it specific to characteristic zero. Over \mathbb{Z}/n\mathbb{Z}, nmod_mpoly_pow_ui multiplies repeatedly (nmod_mpoly_pow_rmul), and for prime n = p uses B^{p} = \sum b_i x^{p\beta_i} (the Frobenius map, since b^p = b in \mathbb{F}_p) to handle the base-p digits of k.

Generic and specialised implementations#

The generic implementation is gr_mpoly_mul, for any coefficient ring given by a gr context: it uses gr_mpoly_mul_monomial for a monomial factor, the heap method gr_mpoly_mul_heap (and gr_mpoly_mul_heap_threaded when the product has more than 10^4 term products and threads are available), and the symmetric squaring gr_mpoly_sqr_commutative_heap for B^2 in commutative rings. It has no array or Kronecker methods, and accumulates coefficients with the ring's own operations. The specialised modules add:

The choice of method in fmpz_mpoly_mul, with its thresholds, is described in Multiplication; nmod_mpoly_mul follows the same pattern.

• Heap: Algorithm 1 in _fmpz_mpoly_mul_johnson (one-word exponents: _fmpz_mpoly_mul_johnson1), with the accumulation of Coefficient accumulation; parallel splitting in fmpz_mpoly_mul_heap_threaded with mpoly_search_monomials.
• Array: fmpz_mpoly_mul_array, _fmpz_mpoly_mul_array_chunked_LEX and _fmpz_mpoly_mul_array_chunked_DEG, the accumulators _fmpz_mpoly_addmul_array1_slong1, _fmpz_mpoly_addmul_array1_slong2, _fmpz_mpoly_addmul_array1_slong and _fmpz_mpoly_addmul_array1_fmpz (Lemma 3).
• Kronecker: fmpz_mpoly_mul_dense (Lemma 4).
• Powering: fmpz_mpoly_pow_fps, fmpz_mpoly_pow_ui, fmpz_mpoly_pow_fmpz (Theorem 6).