MULTIDIGIT MODULAR MULTIPLICATION WITH

THE EXPLICIT CHINESE REMAINDER THEOREM

Daniel J. Bernstein

950518 (draft T)

Abstract. Fix coprime moduli m_{1}, …, m_{s}, of a few digits each. Let n be an integer of a few hundred digits. We show how arithmetic modulo n may be performed upon integers u represented as vectors (u mod m_{1}, …, u mod m_{s}). This method involves no multiprecision arithmetic, except in an easy precomputation; it is practical in software and extremely well suited for hardware. Our main tool is the Explicit Chinese Remainder Theorem, which says exactly how u differs from a particular linear combination of the remainders u mod m_{i}.

1. Introduction

Let m_{1}, …, m_{s} be positive integers. One may reconstruct a sufficiently small integer u from the sequence (u mod m_{1}, u mod m_{2}, …, u mod m_{s}), which we call the remainder representation of u; here u mod m_{i} = u − m_{i} ⌊u/m_{i}⌋ is the remainder when u is divided by m_{i}.

We may calculate uv mod m_{i} as (u mod m_{i})(v mod m_{i}) mod m_{i}. So we may find the remainder representation of uv from the remainder representations of u and v. This is an example of remainder arithmetic. Remainder arithmetic is particularly convenient in hardware since each m_{i} is handled independently.

In this paper we consider the problem of performing arithmetic modulo an integer n. To fix ideas we consider the problem of computing x^{k} mod n, where n is a few hundred digits long; the sufficiently imaginative reader will perhaps be able to find some use for such a computation.

Usual solution: Start from x mod n. Multiply it by itself to obtain (x mod n)^{2}. Reduce that modulo n to obtain x^{2} mod n. Continue with an appropriate sequence of multiplications and reductions so as to end up with x^{k} mod n. (See [3, section 4.6.3] and [1, section 1.2].)

All we have to do, then, is multiply and reduce integers modulo n. Conventional wisdom states that it is pointless to use the remainder representation here, since one cannot compute xy mod n without first, in effect, converting xy back to its usual representation as a binary expansion.

The conventional wisdom is wrong. In this paper we present a new method of reducing a remainder representation modulo n. We do not generally obtain

Typeset by AMS-TEX

1991 Mathematics Subject Classification. Primary 11Y16.

This paper was included in the author’s thesis at the University of California at Berkeley. The author was supported in part by a National Science Foundation Graduate Fellowship.

the smallest possible result, but we do obtain a small enough result that we may perform further operations safely. This method is highly parallel; we draw it to the attention of hardware designers. It is also practical in software.

The crux of our method is a fact that we dub the Explicit Chinese Remainder Theorem. Write P = m_{1}m_{2} ⋯ m_{s}. The usual Chinese Remainder Theorem states that if m_{i} and m_{j} are coprime for every i ≠ j then the remainder representation of u determines u mod P. In fact, u mod P is a particular linear combination of the remainders u mod m_{i}. The Explicit Chinese Remainder Theorem says, in a computationally useful form, exactly what multiple of P must be subtracted from that linear combination to obtain u.

The Explicit Chinese Remainder Theorem has been exploited by Montgomery, Silverman, and Couveignes; see [6], [5], and [2]. In fact, our method may be viewed as an optimization of a slight variant of [5, equation 3.4.4].

In section 2 we state and prove the Explicit Chinese Remainder Theorem. In section 3 we show how to compute one of the quantities in this theorem. In section 4 we explain how to convert an integer from the remainder representation to the usual representation. In section 5 we present our method to reduce a remainder representation modulo n. When n is fixed we may precompute various constants; we discuss the precomputation in section 6. In section 7 we discuss possible software and hardware implementation approaches. Everything here has a polynomial analogue, as outlined in section 8.

Remainder arithmetic is commonly called modular arithmetic. See [3, section 4.3.2] for an introduction. See also [3, exercise 4.5.2–15] for a method of computing the inverse of v modulo m_{i}, or more generally the ratio of u and v modulo m_{i}, when v and m_{i} are coprime.

2. The Explicit Chinese Remainder Theorem

Before stating this theorem we establish two pieces of notation. First, we write u ≡ v (mod P) if u mod P = v mod P; in other words, if u − v is divisible by P. Second, we write round z for the unique integer i with |i − z| < 1/2, when such an integer exists.

The Explicit Chinese Remainder Theorem. Fix coprime positive integers m_{1}, m_{2}, …, m_{s} with product P = m_{1}m_{2} ⋯ m_{s}. Let k_{1}, k_{2}, …, k_{s} be integers with k_{i}P/m_{i} ≡ 1 (mod m_{i}). Let u be an integer with |u| < P/2. If x_{i} ≡ k_{i}u (mod m_{i}) and z = ∑_{i} x_{i}/m_{i} then u = P z − P round z.

Observe that u is a linear combination of round z and the modified remainders x_{i}, with coefficients that do not depend on u. We will take advantage of this linearity in section 5.

Proof. We work modulo m_{j}:

P z = ∑_{i} x_{i} \frac{P}{m_{i}} ≡ x_{j} \frac{P}{m_{j}} ≡ k_{j}u\frac{P}{m_{j}} ≡ u (mod m_{j}).

The m_{j}’s are pairwise coprime so P z ≡ u (mod P). Write r = z − u/P; then P r = P z − u is a multiple of P, so r is an integer. But |z − r| = |u/P | < 1/2. Thus r = round z, and u = P z − P r = P z − P round z. □

3. Computing the nearest integer to a sum of rationals

We retain the notation of the Explicit Chinese Remainder Theorem. Write r = round z. Assume that |u| < P/4, so that |z − r| = |u/P | < 1/4.

Montgomery and Silverman [6] observed that to find r, the nearest integer to z, it suffices to have a low-precision approximation to z. They suggested computing z in floating-point arithmetic. We suggest computing z in fixed-point arithmetic:

Lemma 3.1. Let t_{1}, …, t_{s} be real numbers. Let r be an integer with |∑_{i} t_{i} − r| < 1/4. If 2^{a} ≥ 2s and q_{i} = ⌊2^{a}t_{i}⌋ then r = ⌊3/4 + 2^{−a} ∑_{i} q_{i}⌋.

For us t_{i} = x_{i}/m_{i} and z = ∑_{i} t_{i}. We may compute q_{i} by repeatedly doubling x_{i} modulo m_{i} and checking for overflows.

Proof. By construction 0 ≤ 2^{a}t_{i} − q_{i} < 1. Add: 0 ≤ 2^{a} ∑_{i} t_{i} − ∑_{i} q_{i} < s ≤ 2^{a−1}, so |∑_{i} t_{i} − 1/4 − 2^{−a} ∑_{i} q_{i}| ≤ 1/4, so |r − 1/4 − 2^{−a} ∑_{i} q_{i}| < 1/2. Thus 3/4 + 2^{−a} ∑_{i} q_{i} differs from r + 1/2 by less than 1/2. □

4. Converting from the remainder representation

To illustrate the practical use of the Explicit Chinese Remainder Theorem we consider the following problem. Again fix coprime positive integers m_{1}, m_{2}, …, m_{s} with product P = m_{1}m_{2} ⋯ m_{s}. Let u be an integer with |u| < P/4; say we know u’s remainder representation (u mod m_{1}, u mod m_{2}, …, u mod m_{s}). How do we reconstruct the usual representation of u as a binary expansion?

Precomputation: Find P/m_{i} mod m_{i} = ∏_{j≠i} m_{j} mod m_{i} by successive multiplications. The result is coprime to m_{i}; invert it modulo m_{i} to obtain k_{i}. Do this for each i. (See section 6 for another method.)

Next: For each i find an integer x_{i} ≡ k_{i}u (mod m_{i}). There are at least three plausible ways to do this: (1) set x_{i} = k_{i}u; (2) set x_{i} = k_{i}u mod m_{i}; (3) let x_{i} be the least remainder—either k_{i}u mod m_{i} or (k_{i}u mod m_{i}) − m_{i}, whichever has smaller absolute value. We will use method (2).

Next: Compute P z = ∑_{i} x_{i}(P/m_{i}). We could compute P on the fly, divide by m_{i} and multiply by x_{i}, and sum the results. We could precompute P, or a table of P/m_{i}, in a variety of ways. Several further ideas are discussed in [3, section 4.3.2]. A reasonable approach in practice is to add up z = x_{1}/m_{1} + ⋯ + x_{s}/m_{s} as a multiprecision rational number, either with binary splitting or term by term, and then multiply by P. The term-by-term version amounts to keeping track of the integers m_{1} ⋯ m_{i} and (x_{1}/m_{1} + ⋯ + x_{i}/m_{i})m_{1} ⋯ m_{i} as i increases. The binary splitting version is useful for large s in conjunction with fast multiplication. It is also useful in hardware.

Finally: Compute round z by the method of section 3. Multiply round z by P. Subtract P round z from P z to obtain u.

5. Reducing remainder representations modulo n

Again let u be an integer with |u| < P/4. Given the remainder representation of u we will, without multiprecision arithmetic, construct the remainder representation of an integer v congruent to u modulo n. Although |v| is not in general smaller than n, it is not too much bigger—it is less than n ∑_{i} m_{i}. Our method is expressed in the following corollary of the Explicit Chinese Remainder Theorem:

Lemma 5.1. Fix coprime positive integers m_{1}, m_{2}, …, m_{s} with product P. Let k_{1}, k_{2}, …, k_{s} be integers such that k_{i}P/m_{i} ≡ 1 (mod m_{i}). Let u be an integer with |u| < P/2. Write x_{i} = k_{i}u mod m_{i} and r = round ∑_{i} x_{i}/m_{i}. Then u ≡ v (mod n), where

v = ∑_{i} x_{i} (\frac{P}{m_{i}} mod n) − (P mod n)r.

Furthermore

v ≡ ∑_{i} x_{i} (\frac{P}{m_{i}} mod n mod m_{j}) − (P mod n mod m_{j})r (mod m_{j})

and |v| < n ∑_{i} m_{i}.

Here we see the advantage of focusing on u and v, which are linear combinations of x_{1}, …, x_{s}, r, rather than on u mod n = v mod n as in [6], [5], and [2].

Proof. u = ∑_{i} x_{i}P/m_{i} − P r ≡ v (mod n) by the Explicit Chinese Remainder Theorem.

The formula for v modulo m_{j} follows from the definition of v.

Finally, v is the difference of two nonnegative terms. The first is less than ∑_{i} m_{i}n, since x_{i} < m_{i}. The second is less than ns since r = round ∑_{i} x_{i}/m_{i} ≤ round ∑_{i} 1 = round s = s. □

It is easy to use these formulas in practice. Given n and the various m_{i} we precompute tables of k_{i}, P/m_{i} mod n mod m_{j}, and P mod n mod m_{j}, as explained in the next section. Then, given the remainder representation of u, we compute x_{i} directly as in section 4, and we compute r as in section 3. Finally, for each j, we compute v mod m_{j} as an appropriate dot product modulo m_{j}.

In [6] and [2], x_{i} was chosen as k_{i}u. Like [5] we set x_{i} = k_{i}u mod m_{i}; this greatly improves our v bound for very little effort. We could reduce |v| down to at most (n/4) ∑_{i} m_{i} by using least remainders for x_{i} and for various numbers modulo n. If it is worth that much hassle to save two bits, then perhaps we should go even further, manipulating all the x_{i} to make v small.

6. Precomputing everything

To use the method of section 5 we must precompute the s^{2} + s quantities P/m_{i} mod n mod m_{j} and P mod n mod m_{j}, as well as the inverse of P/m_{j} modulo m_{j}. We can do all this without too much trouble, if P ≥ 4nm_{j} for each j:

Lemma 6.1. Fix coprime positive integers m_{1}, m_{2}, …, m_{s} with product P. Let n be a positive integer. Assume that P/m_{j} ≥ 4n. For i ≠ j let c_{ij} be the inverse of m_{i} modulo m_{j}. Set k_{j} = ∏_{i≠j} c_{ij} mod m_{j}, n_{j} = n mod m_{j}, p_{j} = P mod n mod m_{j}, q_{j} = ⌊P/n⌋ mod m_{j}, and e_{ij} = q_{i}n_{j}+p_{j} mod m_{j}. For i ≠ j set d_{ij} = c_{ij}e_{ij} mod m_{j} and x_{ij} = k_{j}e_{ij} mod m_{j}. Set r_{j} = round ∑_{i≠j} x_{ji}/m_{i} and let d_{jj} be the ratio modulo m_{j} of ∑_{i≠j} c_{ij}x_{ji}−r_{j} and k_{j}. Then k_{j}P/m_{j} ≡ 1 (mod m_{j}) and q_{j}n_{j} +p_{j} ≡ 0 (mod m_{j}) for every j; and d_{ij} = P/m_{i} mod n mod m_{j} for every i and j. Also | ∑_{i≠j} x_{ji}/m_{i} − r_{j}| < 1/4.

Proof. First k_{j}P/m_{j} ≡ ∏_{i≠j} c_{ij}m_{i} ≡ 1 (mod m_{j}).

Abbreviate q = ⌊P/n⌋, so that P = qn + (P mod n). Then q ≡ q_{j} (mod m_{j}) so q_{j}n_{j} + p_{j} ≡ qn + (P mod n) = P ≡ 0 (mod m_{j}).

Next qn ≡ q_{i}n (mod m_{i}n) so P ≡ q_{i}n + (P mod n) (mod m_{i}n). But 0 ≤ q_{i}n + (P mod n) ≤ (m_{i} − 1)n + (n − 1) < m_{i}n. Hence P mod m_{i}n = q_{i}n + (P mod n). Reduce modulo m_{j}:

m_{i} (\frac{P}{m_{i}} mod n) = P mod nm_{i} = q_{i}n + (P mod n) ≡ q_{i}n_{j} + p_{j} ≡ e_{ij} (mod m_{j}).

To prove that d_{ij} = P/m_{i} mod n mod m_{j} we consider separately i ≠ j and i = j. For i ≠ j we have m_{i}d_{ij} ≡ m_{i}c_{ij}e_{ij} ≡ e_{ij} ≡ m_{i}(P/m_{i} mod n) (mod m_{j}). Since m_{i} is coprime to m_{j} we have d_{ij} ≡ P/m_{i} mod n (mod m_{j}).

To handle P/m_{j} mod n mod m_{j} we will apply the Explicit Chinese Remainder Theorem to all m_{i} except m_{j}. The product of this restricted set is P/m_{j}, and we have (k_{i}m_{j})(P/m_{j})/m_{i} ≡ 1 (mod m_{i}) for i ≠ j. Define u = P/m_{j} mod n; observe that |u| < n ≤ P/4m_{j}. Furthermore

x_{ji} ≡ k_{i}e_{ji} ≡ k_{i}m_{j}(P/m_{j} mod n) = (k_{i}m_{j})u (mod m_{i}).

Hence u = (P/m_{j})z − (P/m_{j}) round z with z = ∑_{i≠j} x_{ji}/m_{i}. By construction r_{j} = round z, so |z − r_{j}| = |um_{j}/P | < 1/4. Finally

k_{j}u = k_{j} \frac{P}{m_{j}} ∑_{i≠j} \frac{x_{ji}}{m_{i}} − k_{j} \frac{P}{m_{j}} r_{j} = k_{j} ∑_{i≠j} x_{ji} \frac{P}{m_{i}m_{j}} − k_{j} \frac{P}{m_{j}} r_{j}

≡ k_{j} ∑_{i≠j} x_{ji} \frac{P}{m_{i}m_{j}} m_{i}c_{ij} − k_{j} \frac{P}{m_{j}} r_{j} = k_{j} ∑_{i≠j} c_{ij}x_{ji} \frac{P}{m_{j}} − k_{j} \frac{P}{m_{j}} r_{j}

≡ ∑_{i≠j} c_{ij}x_{ji} − r_{j} ≡ k_{j}d_{jj} (mod m_{j}).

Hence d_{jj} ≡ u = P/m_{j} mod n (mod m_{j}) as desired. □

We may compute all the quantities in Lemma 6.1 exactly as they are defined. It is convenient to store x_{ij} on top of e_{ij}. To find r_{j} we use the method of section 3.

We use multiprecision arithmetic in computing ⌊P/n⌋ and P mod n, and in reducing n, ⌊P/n⌋ , P mod n modulo each m_{j}. We can save some time and effort by computing p_{j} as −q_{j}n_{j} mod m_{j}. Alternatively, when n_{j} is coprime to m_{j}, we can compute q_{j} as the ratio of −p_{j} and n_{j} modulo m_{j}, as suggested in [2]. (The situation of [2] was that P had far more digits than n, so it was impractical to compute ⌊P/n⌋. Fortunately, n was coprime to P, so n_{j} was coprime to m_{j}.)

7. Implementation ideas

We return to the practical problem of computing x^{k} mod n.

Select moduli m_{1}, m_{2}, …, m_{s} with product P ≥ 4(n ∑_{i} m_{i})^{2}. Precompute the quantities described in section 6. Now, given the remainder representations of integers t, t′ smaller than n ∑_{i} m_{i} in absolute value, we may (1) form the remainder representation of tt′, which is an integer of magnitude less than P/4; and (2) reduce tt′ modulo n by the method of section 5. We end up with the remainder representation of a new integer v ≡ tt′ (mod n). Since v is smaller than n ∑_{i} m_{i} in absolute value, it can participate in further multiplications.

So we start from the remainder representation of x mod n, which is smaller than n ∑_{i} m_{i}. We repeatedly multiply and reduce in the remainder representation, keeping all intermediate results below n ∑_{i} m_{i}, until we have reached v ≡ x^{k} (mod n). We convert back to the usual representation of v by the method of section 4. Finally we compute v mod n = x^{k} mod n.

The same strategy works for arbitrary sequences of ring operations modulo n. For example, we can compute multiples of points on an elliptic curve modulo n, compute convolutions modulo n as in [6] or [5], and so on. We simply have to make sure that P is large enough for the intermediate results.

Software. In practice we spend almost our entire time computing

v mod m_{j} = (∑_{i} x_{i} (\frac{P}{m_{i}} mod n mod m_{j}) − (P mod n mod m_{j})r) mod m_{j}

as per Lemma 5.1. This is a matrix-vector product: we multiply our precomputed matrix, with entries of the form P/m_{i} mod n mod m_{j} and −(P mod n) mod m_{j}, by the vector x_{1}, x_{2}, …, x_{s}, r. We could reduce the jth component of the output modulo m_{j} in each step, to keep the numbers small, or in one sweep at the end.

Existing computers perform many operations more efficiently than multiplication and addition modulo m_{j}. In light of the restricted range of inputs we could easily replace all arithmetic with memory lookups.

Extreme example: We could precompute the outputs x_{i}(P/m_{i} mod n) mod m_{j} for each possible x_{i}. Then, given x_{1}, x_{2}, …, x_{s}, we can simply look up and add up the corresponding outputs. If n is fixed for a long time this might be worthwhile. Less extreme would be, e.g., to break each x_{i} into a sum of terms of the form (2^{b})(1, 3, 5, or 7) and to precompute corresponding tables. This may be worthwhile on computers with slow multiplication.

Another helpful technique is to insert the matrix entries into carefully written machine code. On some computers it will be possible to keep most or all of the x_{i} in registers at once. On computers with very fast arithmetic it would be worthwhile to eliminate housekeeping by expanding all loops out of the machine code.

Optimal modulus sizes will depend heavily on such choices and on further details of the computer’s arithmetic facilities.

Note that we could store x_{i} = k_{i}u mod m_{i} rather than u mod m_{i}, and we could modify our precomputed matrix to produce k_{j}v mod m_{j} rather than v mod m_{j}.

Our moduli are not engraved in stone. Given n we could perhaps look for moduli permitting an unusually fast matrix-vector product: for example, where our matrix has many zeros or many repeated entries. Does a “fast matrix” exist for any n? If so, how could we find it? In some situations n itself is not fixed. How quickly can we find a good (n, m_{1}, m_{2}, …, m_{s})?

Hardware. Our reduction method is so dramatically parallel that it cries out for a hardware implementation. The author imagines several levels of chips, of increasing complexity and increasing effectiveness:

1. A straightforward integer matrix-vector multiplication chip.

2. A chip that multiplies a matrix by a vector and reduces the jth component of the output modulo m_{j}. It would be reasonable for the selection of m_{j} to be hardwired into the chip.

3. A chip that understands the entire procedure of section 5: given the remainder representation of u, it computes x_{i}, r, and the remainder representation of v.

4. Like level 3, but the chip stores several remainder representations in named “registers,” and can add, subtract, multiply, or reduce registers upon command.

5. Level 4 plus high-level operations such as modular exponentiation.

6. Level 5 plus conversion to and from the remainder representation.

7. Level 6 plus the precomputation procedure of section 6.

8. The polynomial analogue

There is an analogue of the Explicit Chinese Remainder Theorem for polynomials over a field:

The Explicit Chinese Remainder Theorem for Polynomials. Fix coprime polynomials m_{1}, m_{2}, …, m_{s} with product P = m_{1}m_{2} ⋯ m_{s}. Let k_{1}, k_{2}, …, k_{s} be polynomials such that k_{i}P/m_{i} ≡ 1 (mod m_{i}). Let u be a polynomial with deg u < deg P. Set x_{i} = k_{i}u mod m_{i}. Then u = ∑_{i} x_{i}P/m_{i}.

Note that “round z” for integers has been replaced by 0 for polynomials. Our integer algorithms, mutatis mutandis, work for polynomials.

Proof. Work modulo m_{j}: ∑_{i} x_{i}P/m_{i} ≡ x_{j}P/m_{j} ≡ k_{j}uP/m_{j} ≡ u (mod m_{j}). The m_{j}’s are pairwise coprime so ∑_{i} x_{i}P/m_{i} ≡ u (mod P). By construction deg x_{i} < deg m_{i}. Then both ∑_{i} x_{i}P/m_{i} and u have degree smaller than deg P, so they must be equal. □

When all m_{i} are monic and linear, the Explicit Chinese Remainder Theorem for Polynomials reduces to the familiar fact that a polynomial is a linear combination of its values at a sufficiently large set of points. See [3, pages 484–486].

References

1. Henri Cohen, A Course in Computational Algebraic Number Theory, Springer-Verlag, Berlin, 1993.

2. Jean-Marc Couveignes, Computing a square root for the number field sieve, in [4], 95–102.

3. Donald E. Knuth, The Art of Computer Programming, volume 2: Seminumerical Algorithms, 2nd edition, Addison-Wesley, Reading, Massachusetts, 1981.

4. Arjen K. Lenstra, Hendrik W. Lenstra, Jr. (editors), The development of the number field sieve, Lecture Notes in Mathematics 1554, Springer-Verlag, Berlin, 1993.

5. Peter L. Montgomery, An FFT Extension of the Elliptic Curve Method of Factorization, Dissertation, University of California at Los Angeles, 1992.

6. Peter L. Montgomery, Robert D. Silverman, An FFT Extension to the P − 1 Factoring Algorithm, Mathematics of Computation 54 (1990), 839–854.

5 Brewster Lane, Bellport, NY 11713