Montgomery multiplication (REDC)

Montgomery's method [Mon1985] replaces the division in ab \bmod n by a division by a power of two, which is free. The price is that residues are kept in a scaled form, and that the modulus must be odd. FLINT implements it in nmod.h for full words (R = 2^{\mathtt{FLINT\_BITS}}) and half words (R = 2^{\mathtt{FLINT\_BITS}/2}). A textbook treatment is in [BZ2011].

Throughout, B = 2^{\mathtt{FLINT\_BITS}}, n is odd, and R is a power of two with n < R \le B.

Montgomery form#

The Montgomery form of x \in \mathbb{Z}/n\mathbb{Z} is \tilde x = xR \bmod n. Since R is invertible modulo n, this is a bijection, and

Addition and subtraction are therefore unchanged, and multiplication needs the map T \mapsto T R^{-1} \bmod n, called REDC. Conversion into Montgomery form costs one ordinary reduction of xR, and conversion back is one REDC of \tilde x.

Reduction#

Let n' = -n^{-1} \bmod R, which exists because n is odd. FLINT stores it in the field nred of an nmod_redc_ctx_t and computes it with n_binvert.

With R = B and T = ab for reduced a and b, step 1 is a low product, mn is a full product, and the division by R is taking the high word.

The subtractive variant#

n_ll_redc, used by nmod_redc_mul, works with n^{-1} \bmod R = -n' instead of n'. This avoids the carry from the addition in step 2, which matters when n is close to B and T + mn does not fit in two words.

Every intermediate value fits in a word, so this works for all odd moduli up to B - 1. It is what the plain nmod_redc_ functions use.

Lazy reduction#

Skipping step 3 of Algorithm 1 saves a comparison and a conditional subtraction per multiplication, at the cost of an output in [0, 2n) instead of [0, n). Such lazy residues can be fed back into the next multiplication if the modulus is small enough. This observation is due to Walter [Wal1999].

The fast functions use this with R = B, so they require n < 2^{\mathtt{FLINT\_BITS} - 2}; nmod_redc_can_use_fast tests exactly this. The condition is needed: for 2^{\mathtt{FLINT\_BITS} - 2} < n < 2^{\mathtt{FLINT\_BITS} - 1} and a, b close to 2n, the output of nmod_redc_fast_mul can exceed 2n, and the sum T + mn can overflow two words, so that the result is not even congruent to ab R^{-1}. The bound also keeps all intermediate values in range: T + mn < 4n^2 + Bn < B^2, and 2n < B for the lazy addition and subtraction modulo 2n.

A lazy residue is turned into an ordinary one by subtracting n once if needed (nmod_redc_fast_normalise). Comparing two lazy residues for equality requires normalising both first.

Half-word arithmetic#

With R = 2^{\mathtt{FLINT\_BITS}/2} and n < R, the product T = ab of two reduced operands fits in one word, and so does T + mn provided that T + (R - 1) n < B. REDC then needs only single-word multiplications: step 1 is a product reduced modulo R (a mask), step 2 a product, an addition and a shift.

• For a, b < n (the half functions) we need n^2 + (R - 1)n < B, which holds for n < R/2 = 2^{\mathtt{FLINT\_BITS}/2 - 1}.
• For lazy a, b < 2n (the half_fast functions) we need 4n \le R by Theorem 3, so n < 2^{\mathtt{FLINT\_BITS}/2 - 2}. Then also 4n^2 + Rn < B.

On 64-bit machines this limits the modulus to 31 or 30 bits. Whether the half-word version is faster depends on the processor; it saves the two-word arithmetic but not the number of multiplications.

Costs#

The number of products is the same as with Möller–Granlund reduction, but Montgomery reduction needs no normalisation shifts and fewer corrections. This shortens the chain of dependent operations, which matters most in exponentiation and other long chains of dependent products. The conversions cost about one ordinary modular multiplication each, so Montgomery arithmetic pays off only when many operations are done between them.

Use in FLINT#

• nmod_pow_ui switches to Montgomery arithmetic for odd n and exponents of at least 2^{11}, and nmod_2_pow_ui for exponents of at least 2^{20}. Both use the lazy variant when n < 2^{\mathtt{FLINT\_BITS} - 2}.
• n_is_prime performs its strong probable prime tests for n \ge 2^{32} in Montgomery form, and n_factor_pollard_brent_single iterates x \mapsto x^2 + c in Montgomery form.
• The multimodular algorithms for Bell numbers (arith_bell_number_multi_mod) and Euler numbers (in arb_euler_number_ui) use lazy Montgomery arithmetic modulo each prime.
• nmod_sqrt does all its arithmetic in Montgomery form, using the half_fast, fast or plain functions depending on the size of n.
• _nmod_vec_scalar_mul_nmod uses nmod_redc_mul for full-word odd moduli and vectors of length at least 8, where Shoup's method is not available.
• nmod_poly computes Euclidean gcds and resultants of medium length in a gr_ctx_init_nmod_redc_fast context when 2^{\mathtt{FLINT\_BITS}/2 - 1} \le n < 2^{\mathtt{FLINT\_BITS} - 2}.