Multimodular reduction and Chinese remaindering

A multimodular algorithm computes a result x (an integer, or the entries of a matrix or the coefficients of a polynomial) by reducing its inputs modulo many word-size primes m_1, \ldots, m_k, solving the problem modulo each prime, and recovering x from its residues by the Chinese remainder theorem, which determines x modulo P = m_1 \cdots m_k. The two conversions, multi-mod (from x to its residues) and multi-CRT (back), each cost O(k^2) limb operations if done naively; the subproduct tree makes them quasi-linear [BorMoe1974; Ber2008].

The generic and fixed-size Chinese remaindering used inside the small-prime FFT is described with crt_helpers; multimodular algorithms for matrices and polynomials are discussed in Multimodular and p-adic linear algebra over the integers and Multiplication with several primes and CRT.

The subproduct tree#

Group the moduli into the leaves of a balanced binary tree, and label each node with the product v of the moduli below it, so that the root carries P. With N the bit size of P, each level of the tree has total size about N bits, there are O(\log k) levels, and the tree is computed bottom-up in time O(M(N) \log k).

Since v_b \mid v, (x \bmod v) \bmod v_b = x \bmod v_b, so the algorithm is correct; each level costs divisions of total size O(N), so the whole costs O(M(N) \log k) with fast division.

Chinese remaindering with fractional cofactors#

The classical solution x \equiv \sum_i r_i \, c_i \, (P/m_i) \pmod P with c_i = (P/m_i)^{-1} \bmod m_i needs the cofactors P/m_i, of N bits each. The tree form combines the children b and c of a node a with v_a = v_b v_c by A = A_b v_c + A_c v_b, which never forms anything larger than v_a. For this to work the values must carry the right scaling.

The values R themselves are computed from the root down by a remainder tree: R_b = \bigl((R_a \bmod v_b) \, v_c\bigr) \bmod v_b, which follows from R_b \equiv R_a v_c \pmod{v_b} and v_b \mid v_a. At the bottom of the tree only the inverses R^{-1} modulo the leaves are needed, and a leaf product fits in a word, so no multiprecision inverse is ever computed. The precomputation therefore costs O(M(N) \log k): a product tree and a remainder tree.

So the values grow by about one bit per level, and the combination is a pair of unbalanced products and one addition per node, again O(M(N) \log k) in total.

Refinements in FLINT#

• Batched leaves. Moduli of at most 32 bits are multiplied in pairs, and 20-bit moduli in triples, so that the tree works with full-word leaves; single-word code splits and combines the individual moduli at the bottom.
• Basecases. Below a chunk of about 2048 bits the CRT is the linear sum of step 1 (by mpn_addmul_1), and below a node of about 2048 bits the reduction computes x \bmod m_j for all leaves at once as a dot product of the limbs of x with precomputed powers B^i \bmod m_j, with a two-limb accumulator and a single reduction per leaf.
• Fixed-length code. When P is a single small chunk, the CRT uses the straight-line kernels of crt_helpers.h.
• Small values. Results are often much smaller than P (for example the small entries of a matrix whose bound is dominated by a few large ones). The reconstruction first reconstructs x modulo the product Q of the first few moduli. If the result x' has at least a limb of slack below Q and x' \bmod m_i = x_i for all the other moduli, then x' \equiv x \pmod P and 0 \le x' < Q \le P, so x' = x and the tree is not traversed. A random residue passes the slack test with probability about 2^{-64}, so this costs almost nothing when it fails.
• Division. Reductions down the tree use precomputed inverses of the node products (flint_mpn_mod_preinvn) above 16 limbs.
• Vectors. For many values, the tables are traversed once per block of values (flint_mpn_multi_mod_vec, flint_mpn_multi_crt_vec).
• Memory. The structure stores O(N \log k) bits. The one-shot reconstruction flint_mpn_multi_crt_once instead traverses the tree depth first, recomputing subproducts as needed, with O(N) memory, which makes reconstructions of billions of bits possible. Large nodes are processed in parallel.

Use in FLINT#

The functions are flint_mpn_crt_init, flint_mpn_multi_mod, flint_mpn_multi_crt and their variants in mpn_extras (Multimodular reduction and Chinese remaindering). They are used through fmpz_comb_t (fmpz_multi_mod_ui, fmpz_multi_CRT_ui) by the multimodular algorithms for fmpz_mat and mpn_mod matrices (multiplication), fmpz_poly and fmpq_poly (GCDs, interpolation, Taylor shifts), resultants of bivariate polynomials, and by the multimodular computation of Bernoulli, Euler and Bell numbers and Stirling numbers in arith.