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.