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.