Modular multiplication with a multiword modulus

FLINT computes $ab \bmod n$ for residues $0 \le a, b < n$ when $n$ has $\ell \ge 2$ limbs in three places: in mpn_mod (mpn_mod_mul), in mpn_extras (flint_mpn_mulmod_preinvn and its relatives, which mpn_mod, fmpz and the primality tests use) and in fmpz_mod for two-limb moduli. All methods follow Barrett's idea [Bar1987] of replacing the division by $n$ with multiplications by a precomputed approximation of $1/n$; for a single limb, the analogous methods are described in Division by an invariant word, Barrett and Lemire reduction and Multiplication by a fixed operand (Shoup). A textbook treatment of multiword Barrett and Montgomery reduction is [BZ2011, §2.4].

Throughout, $B = 2^{\mathtt{FLINT\_BITS}}$, $n$ has exactly $\ell$ limbs ($B^{\ell - 1} \le n < B^{\ell}$), and $\nu$ is the number of leading zero bits of its top limb.

Normalisation#

The reduction works with the normalised modulus $D = n 2^{\nu}$, which satisfies $B^{\ell}/2 \le D < B^{\ell}$. Reducing modulo $D$ is enough, since for every integer $T$

To reduce a product $ab$, one reduces $T = (2^{\nu} a) b$ modulo $D$ and shifts the remainder right by $\nu$ bits. mpn_mod stores $D$ in its context (MPN_MOD_CTX_MODULUS_NORMED) and shifts $a$ before the multiplication; flint_mpn_mulmod_preinvn instead expects operands that are already shifted.

The precomputed inverse#

flint_mpn_preinvn computes the $\ell$-limb number

Since $B^{\ell}/2 < D + 1 \le B^{\ell}$, the quotient $\mu = \lfloor B^{2\ell}/(D+1) \rfloor$ lies in $[B^{\ell}, 2B^{\ell})$, so $0 \le v < B^{\ell}$ and $\mu = B^{\ell} + v$ is an approximation of $B^{2\ell}/D$ from below whose leading 1 need not be stored. (For $D = B^{\ell} - 1$, $v = 0$.) This is the $\ell$-limb analogue of the one-limb inverse of Möller and Granlund [GraMol2010]. It costs one division of a $2\ell$-limb number, which mpn_mod performs once per context.

Reduction of a product#

The input bound holds for $T = (2^{\nu} a) b$ with $a, b < n$, since then $T < D n < D B^{\ell}$. It implies $T_1 < D$.

The bound 5 is crude; it only matters that the number of corrections is bounded by a small constant. The cost is one $\ell \times \ell$ product (to form $T$), one high product $T_1 v$ and one low product $q D$ of $\ell + 1$ limbs, which is about three products when all three are computed in full. flint_mpn_mul_or_mulhigh_n and flint_mpn_mul_or_mullow_n compute the high and low products with the short product algorithms of mpn_extras when these are faster (High, low and middle products of integers). The approximate high product is a lower bound that is short by at most one unit (flint_mpn_mulhigh_n), so $q$ is still at most $Q$ and the loop in step 4 runs at most once more.

For $\ell = 2$, flint_mpn_mulmod_preinvn_2 (used by mpn_mod_mul and the two-limb vector loops) is the same algorithm written with inline two-limb arithmetic; the low product $qD$ is computed as three limbs.

Variants#

• flint_mpn_fmmamod_preinvn computes $a_1 b_1 + a_2 b_2 \bmod n$. The shifted sum $T = 2^{\nu}(a_1 b_1 + a_2 b_2)$ is less than $2 D n < 2 D B^{\ell}$; one conditional subtraction of $D B^{\ell}$ (from the high half) brings it below $D B^{\ell}$, and the algorithm above finishes. This costs two products and one reduction, against two products and two reductions for separate modular products.
• flint_mpn_mod_preinvn and flint_mpn_divrem_preinvn reduce numbers of any length $m \ge \ell$ by applying the algorithm to $\ell$ limbs of the dividend at a time, from the top.
• mpn_mod_set_mpn, which reduces unreduced sums in mpn_mod, uses flint_mpn_mod_preinvn for $\ell \ge 3$ and GMP's mpn_divrem_2 for $\ell = 2$.

Two-limb Barrett reduction in fmpz_mod#

For $B < n < B^2$, fmpz_mod_ctx_init precomputes the three-limb number $m = \lfloor B^4 / n \rfloor$, and fmpz_mod_mul (through _fmpz_mod_mul2) reduces $x = ab < n^2$ by classical Barrett reduction:

The code computes the five-limb product $x_1 m$ with ten word products, skipping only the low half of the lowest one (which cannot carry into the limbs that are used), and the low three limbs of $qn$. The remainder fits in three limbs, so the corrections are three-limb comparisons. The modulus $n = B$ has its own code (_fmpz_mod_mul2s, a wrapping product), and one-limb moduli use nmod_mul.

A fixed multiplier#

When the same $a$ multiplies many $b$, as in scaling a vector or in the division steps of Euclid's algorithm, precomputation for $a$ pays off. flint_mpn_mulmod_want_precond chooses the method from $\ell$ and the number of products: none for fewer than 4 products, Shoup's method for $\ell \le 10$ (and for $\ell \le 12$ with at most 12 products), the matrix method for larger $\ell$ up to 64 limbs (and beyond with enough products), so that mpn_mod uses Shoup's method up to 10 limbs and the matrix method for 11 to 16 limbs. For $\ell = 2$ with the top bit of $n$ set, no precomputation is used through this function, but the vector functions of mpn_mod have their own inlined two-limb Shoup code.

Shoup's method#

Shoup's method [Sho2021; Har2014] precomputes the $\ell$-limb quotient

(flint_mpn_mulmod_precond_shoup_precompute, one division). A product is then computed as

where $r$ is obtained from the low halves of $ab$ and $qn$ alone.

The cost per product is one high product and two low products of $\ell$ limbs, against one full, one high and one low product without precomputation. If $n < B^{\ell}/2$ ($\nu > 0$), then $r < 2n < B^{\ell}$ fits in $\ell$ limbs and one conditional subtraction finishes.

The implementation computes $q$ with the approximate high product flint_mpn_mulhigh_n, a lower bound which may be $q - 1$. This can only happen when $f < c B^{\ell - 1}$ for a small constant $c$ (the error of the approximation, a few units of the limb below the high half). Then the computed remainder is $r + n < n (1 + c/B + b/B^{\ell}) < n (3/2 + c/B) < 2n$ when $b < n < B^{\ell}/2$, so for $\nu > 0$ one subtraction still suffices. If the top bit of $n$ is set, $r$ may exceed $B^{\ell}$; the code then computes modulo $B^{\ell + 1}$ and subtracts $n$ up to twice (for $\ell = 2$ it uses an exact high product and one subtraction).

The matrix method#

The matrix method precomputes the $\ell$ residues

(flint_mpn_mulmod_precond_matrix_precompute, $\ell - 1$ one-limb shifts and reductions). For $b = \sum_i b_i B^i$,

and the right-hand side is a matrix-vector product of $\ell^2$ word products, less than $\ell B D < B^{\ell + 2}$. It is reduced modulo $D$ with two steps of one-limb quotient estimation (flint_mpn_mod_preinv1) and shifted right by $\nu$ bits. Compared with the full product followed by the reduction of Algorithm 1, the high and low products of the reduction are replaced by $O(\ell)$ work, at the price of $\ell^2$ limbs of precomputed data. The fused version flint_mpn_fmmamod_precond_matrix adds two such matrix-vector products before reducing.

Inversion and powering#

Inverses modulo a multiword $n$ are computed with GMP's extended Euclidean algorithm (mpn_gcdext in mpn_mod_inv, fmpz_gcdinv in fmpz_mod_inv), which is subquadratic for large operands. Powers with large exponents use GMP's mpz_powm (in mpn_mod_pow_fmpz, and through fmpz_powm in fmpz_mod_pow_fmpz for moduli below 70000 bits). It uses Montgomery's representation [Mon1985], in which no quotient estimate is needed, with sliding windows; mpn_mod uses it for exponents with more than 6 bits. For moduli of 70000 bits and more (and when fft_small is available), fmpz_mod_pow_fmpz switches to generic sliding-window exponentiation with fmpz_mod_mul, whose products use the FFT.

Use in FLINT#

• mpn_mod_mul, mpn_mod_sqr, _mpn_mod_vec_mul: Algorithm 1, for $\ell = 2$ through flint_mpn_mulmod_preinvn_2.
• mpn_mod_fmma and the length-2 dot products: flint_mpn_fmmamod_preinvn.
• _mpn_mod_vec_mul_scalar, _mpn_mod_vec_addmul_scalar, _mpn_mod_vec_submul_scalar: Shoup's method or the matrix method for vectors of at least 4 (or 6) elements.
• _mpn_mod_poly_divrem_q1_preinv1_fmma_precond and _mpn_mod_poly_divrem_q1_preinv1_karatsuba_precond: the matrix method for the two quotient coefficients of a Euclidean step.
• fmpz_mod_mul for $B < n < B^2$: Proposition 2.
• flint_mpn_mulmod_preinvn, flint_mpn_powmod_preinvn and flint_mpn_sqrtmod_preinv in mpn_extras, used by mpn_mod_sqrt and the integer square root modulo a prime.
• Sums of many products are not reduced one product at a time; see Delayed reduction for multiword moduli.