MODULAR EXPONENTIATION VIA THE EXPLICIT CHINESE REMAINDER THEOREM

DANIEL J. BERNSTEIN AND JONATHAN P. SORENSON

Abstract. Fix pairwise coprime positive integers p_{1}, p_{2}, …, p_{s}. We propose representing integers u modulo m, where m is any positive integer up to roughly \sqrt{p_{1}p_{2} ⋯ p_{s}}, as vectors (u mod p_{1}, u mod p_{2}, …, u mod p_{s}). We use this representation to obtain a new result on the parallel complexity of modular exponentiation: there is an algorithm for the Common CRCW PRAM that, given positive integers x, e, and m in binary, of total bit length n, computes x^{e} mod m in time O(n/lg lg n) using n^{O(1)} processors. For comparison, a parallelization of the standard binary algorithm takes superlinear time; Adleman and Kompella gave an O((lg n)^{3}) expected time algorithm using exp(O(\sqrt{n lg n})) processors; von zur Gathen gave an NC algorithm for the highly special case that m is polynomially smooth.

1. Introduction

In this paper we consider the problem of computing x^{e} mod m for large integers x, e, and m. This is the bottleneck in Rabin’s algorithm for testing primality, the Diffie-Hellman algorithm for exchanging cryptographic keys, and many other common algorithms. See, e.g., [18, Section 4.5.4].

The usual solution for computing x^{e} mod m is to compute small integers that are congruent modulo m to various powers of x. See, e.g., [18, Section 4.6.3], [11, Section 1.2], [16], [20], and [9]. For example, say e = 10. One can compute x_{1} = x mod m, then x_{2} = x_{1}^{2} mod m, then x_{4} = x_{2}^{2} mod m, then x_{5} = x_{1}x_{4} mod m, and finally x_{10} = x_{5}^{2} mod m. It is often convenient to allow x_{1}, x_{2}, x_{4}, x_{5} to be slightly larger than m.

The output x^{e} mod m and inputs x, e, m are conventionally written in binary. Standard practice is to also use the binary representation for x_{1}, x_{1}^{2}, x_{2}, x_{2}^{2}, etc. We instead use the residue representation: x_{1} is represented as (x_{1} mod p_{1}, x_{1} mod p_{2}, …, x_{1} mod p_{s}), x_{1}^{2} is represented as (x_{1}^{2} mod p_{1}, x_{1}^{2} mod p_{2}, …, x_{1}^{2} mod p_{s}), etc. Here p_{1}, p_{2}, …, p_{s} are small pairwise coprime positive integers, typically primes, such that P = p_{1}p_{2} ⋯ p_{s} is sufficiently large. Note that, because p_{j} is small, it is easy to compute (e.g.) x_{1}x_{4} mod p_{j} from x_{1} mod p_{j} and x_{4} mod p_{j}.

The usual Chinese remainder theorem says that (for example) x_{1}x_{4} mod P is determined by the residue representation of x_{1}x_{4}. In fact, any integer u is congruent modulo P to a particular linear combination of u mod p_{1}, u mod p_{2}, …, u mod p_{s}.

Date: 2005.11.09. Permanent ID of this document: 42983bacd3ca91e10217a79c66508b3b.

2000 Mathematics Subject Classification. Primary 11Y16; Secondary 68W10.

This paper combines and improves the preliminary papers [7] by Bernstein and [29] by Sorenson. Bernstein was supported by the National Science Foundation under grants DMS-9600083 and DMS–9970409. Sorenson was supported by the National Science Foundation under grant CCR–9626877. Sorenson completed part of this work while on sabbatical at Purdue University in Fall 1998.

Residue representations of x_{1} and x_{4}

easy

[figure] Residue representation of x_{1}x_{4}

CRT

[figure]

ECRT

[figure]

\begin{matrix}ECRT \\ mod m \end{matrix}

[figure] ECRT ; mod m ; mod p_{j}

[figure]

Binary representation of something congruent to x_{1}x_{4} modulo P

reduction ; modulo P

[figure] Binary representation of x_{1}x_{4}

reduction ; modulo m [figure] Binary representation of x_{5}

reduction ; modulo p_{j} [figure] Residue representation of x_{5}

The explicit Chinese remainder theorem says, in a computationally useful form, exactly what multiple of P must be subtracted from the linear combination to obtain u. See §2.

Once we know the binary representation of x_{1}x_{4}, we could reduce it modulo m to obtain the binary representation of x_{5}. We actually use the explicit CRT modulo m to obtain the binary representation of x_{5} directly from the residue representation of x_{1}x_{4}. See §3.

Once we know the binary representation of x_{5}, we reduce it modulo each p_{j} to obtain the residue representation of x_{5}. An alternative is to use the explicit CRT modulo m modulo p_{j} to obtain the residue representation of x_{5} directly from the residue representation of x_{1}x_{4}. See §4.

If P is sufficiently large then, as discussed in §5, one can perform several multiplications in the residue representation before reduction modulo m. This is particularly beneficial for parallel computation: operations in the residue representation are more easily parallelized than reductions modulo m. By optimizing parameters, we obtain the following result:

Theorem 1.1. There is an algorithm for the Common CRCW PRAM that, given the binary representations of positive integers x, e, and m, of total bit length n, computes the binary representation of x^{e} mod m in time O(n/lg lg n) using n^{O(1)} processors.

In §6 we define the Common CRCW PRAM, in §8 we present the algorithm, and in §7 we present a simpler algorithm taking time O(n).

For comparison: A naive parallelization of the standard binary algorithm takes superlinear time; see §7. Adleman and Kompella gave an O((lg n)^{3}) expected time algorithm using exp(O(\sqrt{n lg n})) processors; see [2]. Von zur Gathen gave an NC algorithm for the highly special case that m is polynomially smooth; see [35]. It is unknown whether general integer modular exponentiation is in NC. For more on parallel exponentiation algorithms, see Gordon’s survey [16].

Our techniques can easily be applied to exponentiation in finite rings more general than Z/m. For example, the deterministic polynomial-time primality test of Agrawal, Kayal, and Saxena in [4] can be carried out in sublinear time using a polynomial number of processors on the Common CRCW PRAM; the bottleneck is exponentiation in a ring of the form (Z/m)[x]/(x^{k} − 1).

2. The explicit Chinese remainder theorem

Algorithms for integer arithmetic using the residue representation were introduced in the 1950s by Svoboda, Valach, and Garner, according to [18, Section 4.3.2]. This section discusses conversion from the residue representation to the binary representation.

For each real number α such that α − 1/2 ∉ Z, define round α as the unique integer r with |r − α| < 1/2.

Theorem 2.1. Let p_{1}, p_{2}, …, p_{s} be pairwise coprime positive integers. Write P = p_{1}p_{2} ⋯ p_{s}. Let q_{1}, q_{2}, …, q_{s} be integers with q_{i}P/p_{i} ≡ 1 (mod p_{i}). Let u be an integer with |u| < P/2. Let u_{1}, u_{2}, …, u_{s} be integers with u_{i} ≡ u (mod p_{i}). Let t_{1}, t_{2}, …, t_{s} be integers with t_{i} ≡ u_{i}q_{i} (mod p_{i}). Then u = P α − P round α where α = ∑_{i} t_{i}/p_{i}.

Proof. P α = ∑_{i} t_{i}(P/p_{i}) ≡ ∑_{i} u_{i}q_{i}(P/p_{i}) ≡ u_{j}q_{j}(P/p_{j}) ≡ u (mod p_{j}) for each j, so P α ≡ u (mod P). Write r = α − u/P. Then r is an integer, and |r − α| = |u/P | < 1/2, so r = round α, i.e., u/P = α − round α. □

The usual Chinese remainder theorem says that u ≡ ∑_{i} t_{i}P/p_{i} (mod P); this is the integer version of Lagrange’s interpolation formula. This is a popular way to compute u from the residues u_{i}: first compute t_{i} = u_{i}q_{i} mod p_{i} or simply t_{i} = u_{i}q_{i}, then compute P α = ∑_{i} t_{i}P/p_{i}, then reduce P α modulo P to the right range.

The explicit Chinese remainder theorem, Theorem 2.1, suggests another way to divide P α by P. Use t_{i} and p_{i} directly to compute a low-precision approximation to α = ∑_{i} t_{i}/p_{i} with sufficient accuracy to determine round α; see, for example, Theorem 2.2 below. Then subtract P round α from P α to obtain u.

As far as we know, the first use of the explicit Chinese remainder theorem was by Montgomery and Silverman in [23, Section 4]. It has also appeared in [22], [13, Section 2.1], [7], [29], and [3, Section 5].

Theorem 2.2. Let β_{1}, β_{2}, …, β_{s} be real numbers. Let r and a be integers. If |r − ∑_{i} β_{i}| < 1/4 and 2^{a} ≥ 2s then r = ⌊3/4 + 2^{−a} ∑_{i} ⌊2^{a}β_{i}⌋⌋.

Proof. r < 1/4 + ∑_{i} β_{i} = 1/4 + 2^{−a} ∑_{i} 2^{a}β_{i} ≤ 1/4 + 2^{−a}(s + ∑_{i} ⌊2^{a}β_{i}⌋) ≤ 1/4 + 2^{−a}(2^{a}/2 + ∑_{i} ⌊2^{a}β_{i}⌋) = 3/4 + 2^{−a} ∑_{i} ⌊2^{a}β_{i}⌋ ≤ 3/4 + ∑_{i} β_{i} < r + 1. □

In the situation of Theorem 2.1, assume further that |u| < P/4. We can then use Theorem 2.2 to quickly compute r = round α. We choose a with 2^{a} ≥ 2s, then compute the fixed-point approximations 2^{−a} ⌊2^{a}t_{i}/p_{i}⌋ to t_{i}/p_{i}, then compute r = ⌊3/4 + 2^{−a} ∑_{i} ⌊2^{a}t_{i}/p_{i}⌋⌋.

3. The explicit CRT mod m

Theorem 3.1. Let p_{1}, p_{2}, …, p_{s} be pairwise coprime positive integers. Write P = p_{1}p_{2} ⋯ p_{s}. Let q_{1}, q_{2}, …, q_{s} be integers with q_{i}P/p_{i} ≡ 1 (mod p_{i}). Let u be an integer with |u| < P/2. Write t_{i} = uq_{i} mod p_{i}. Let m be a positive integer. Write v = ∑_{i} t_{i}(P/p_{i} mod m) − (P mod m) round ∑_{i} t_{i}/p_{i}. Then u ≡ v (mod m), and |v| ≤ m ∑_{i} p_{i}.

The hypotheses on p, P, q, u are the same as in Theorem 2.1. The hypothesis on t is more restrictive: Theorem 2.1 allowed any integer t_{i} ≡ uq_{i} (mod p_{i}), whereas Theorem 3.1 insists that 0 ≤ t_{i} < p_{i}.

One can allow a slightly wider range of t_{i}, often making t_{i} easier to compute, at the expense of a larger bound on |v|. In the other direction, one can reduce the bound on |v| by taking t_{i} between −p_{i}/2 and p_{i}/2, and one can reduce the bound further by similarly adjusting P/p_{i} mod m, P mod m.

Proof. By Theorem 2.1, u = ∑_{i} t_{i}P/p_{i} − P round ∑_{i} t_{i}/p_{i} ≡ v. Furthermore, v ≤ ∑_{i} t_{i}(P/p_{i} mod m) ≤ ∑_{i} p_{i}m since t_{i} ≤ p_{i}, and −v ≤ (P mod m) round ∑_{i} t_{i}/p_{i} ≤ m ∑_{i} p_{i} since t_{i} ≤ p_{i}^{2}. □

Theorem 3.1 suggests the following representation of integers modulo m. Select moduli p_{1}, …, p_{s} whose product P = p_{1} ⋯ p_{s} exceeds 4(m ∑_{i} p_{i})^{2}. Use the vector (x mod p_{1}, …, x mod p_{s}), where x is any integer between −m ∑_{i} p_{i} and m ∑_{i} p_{i}, to represent x mod m. Note that each element of Z/m has many representations.

The following procedure, given two such vectors (x mod p_{1}, …, x mod p_{s}) and (y mod p_{1}, …, y mod p_{s}), computes another such vector (v mod p_{1}, …, v mod p_{s}) with v ≡ xy (mod m):

(1) Precompute q_{i}, P/p_{i} mod m, and P mod m.

(2) Compute t_{i} = ((x mod p_{i})(y mod p_{i})q_{i}) mod p_{i}, so that t_{i} = (xy)q_{i} mod p_{i}.

(3) Compute round ∑_{i} t_{i}/p_{i} by Theorem 2.2.

(4) Compute v = ∑_{i} t_{i}(P/p_{i} mod m) − (P mod m) round ∑_{i} t_{i}/p_{i}. Then v is between −m ∑_{i} p_{i} and m ∑_{i} p_{i}, and v ≡ xy (mod m), by Theorem 3.1.

(5) Compute the residues (v mod p_{1}, …, v mod p_{s}).

The output can then be used in subsequent multiplications. One can carry out more componentwise operations before applying Theorem 3.1 if P is chosen larger; see §5 for further discussion.

Montgomery and Silverman in [23, Section 4] suggested computing u mod m by first computing v. (Another way to compute u mod m, well suited for FFT-based arithmetic, is to perform the computation of [8, Section 13] modulo m.) The idea of using v for subsequent operations, and not bothering to compute u mod m, was introduced in [7].

4. The explicit CRT mod m mod p_{j}

Theorem 4.1. Under the assumptions of Theorem 3.1,

v ≡ ∑_{i} t_{i}(P/p_{i} mod m mod p_{j}) − (P mod m mod p_{j}) round ∑_{i} t_{i}/p_{i} (mod p_{j}).

Proof. Reduce the definition of v modulo p_{j}. □

In §3 we discussed precomputing P/p_{i} mod m, precomputing P mod m, computing v, and reducing v modulo p_{j}. Theorem 4.1 suggests a different approach: precompute P/p_{i} mod m mod p_{j}; precompute P mod m mod p_{j}; compute v mod p_{j} as (∑_{i} t_{i}(P/p_{i} mod m mod p_{j}) − (P mod m mod p_{j}) round ∑_{i} t_{i}/p_{i}) mod p_{j}. This idea was introduced in [7].

Theorem 4.1 is particularly well suited to parallel computation, as briefly observed in [7]. By using additional ideas described in the remaining sections of this paper, one can exponentiate on a polynomial-size parallel computer in sublinear time, as shown in [29]. This paper includes all the results from [7] and [29].

5. Higher powers

We begin by reviewing the left-to-right base-2 algorithm for computing x^{e} mod m:

Let l denote the number of bits in e; Write e = ∑_{k=0}^{l−1} e_{k}2^{k}, where e_{k} ∈ {0, 1}; y ← 1; For(k ← l − 1; k ≥ 0; k ← k − 1) do:

y ← y^{2}x^{e_{k}} mod m;

Output(y);

More generally, let b be a power of 2. The base-b algorithm is as follows:

Let l denote the number of base-b digits in e; Write e = ∑_{k=0}^{l−1} e_{k}b^{k}, where 0 ≤ e_{k} < b; y ← 1; For(k ← l − 1; k ≥ 0; k ← k − 1) do:

y ← y^{b}x^{e_{k}} mod m;

Output(y);

Sane people take every opportunity to reduce intermediate results modulo m: for example, they compute y^{b} mod m by squaring y, reducing modulo m, squaring again, reducing again, etc. We call this quite reasonable practice into question. It may be better to first compute y^{b}x^{e_{k}}, then reduce the result modulo m. The benefit of performing fewer reductions may outweigh the cost of performing arithmetic on larger numbers. (Similar comments apply to other exponentiation algorithms.)

When we use the explicit Chinese remainder theorem to perform multiplications modulo m, the cost of handling larger numbers is in using more primes p_{1}, p_{2}, . . .. In the context of parallel computation, increasing the number of primes mainly affects the number of processors, not the run time. In fact, as explained in the next three sections, we can save an asymptotically non-constant factor in the run time in various models of parallel computation, at a reasonable expense in the number of processors.

Perhaps the same idea can also save a small constant factor in the serial case.

6. Definition of the CREW PRAM and the Common CRCW PRAM

We use two models of parallel computation in this paper. This section defines the models.

In both models, computers are “parallel random-access machines” (PRAMs) in which many processors can access variable locations in a shared memory. The models differ in how they handle memory conflicts: the “Common CRCW PRAM” allows many processors to simultaneously write a common value to a single memory location, while the “CREW PRAM” does not. Our exponentiation algorithm in §8 relies on the extra power of the Common CRCW PRAM; the simplified algorithm in §7 works in either model.

We caution the reader that phrases such as “CREW PRAM” are not sufficient to pin down a model of computation in the literature. For example, [25, Sections 2.6 and 15.2] defines unrealistic PRAMs in which the processor word size is permitted to grow linearly as a computation progresses; at the other extreme, [14] presents computations that fit into logarithmic-size index registers and constant-size data registers. The complexity of a computation can depend heavily on details of the model. See [15] for a comparison of several models.

Computers. A parallel computer is parametrized by positive integers p, s, w, r, i such that p ≤ 2^{w} and s ≤ 2^{w}. The computer has p processors, numbered from 0 through p − 1, operating on s words of shared memory, labelled mem[0] through mem[s−1]. Each word is a w-bit string, often interpreted as the binary representation of an integer between 0 and 2^{w} − 1. Each processor has r registers, labelled reg[0] through reg[r − 1]; each register contains one word. The computer also has space for i instructions; each processor has an instruction pointer.

Algorithms. A parallel algorithm is a sequence of instructions. (The algorithm fails on a computer of size p, s, w, r, i if it has more than i instructions.) Here are the possible instructions:

(1) Clear j: Set reg[j] ← 0. (The algorithm fails if j ≥ r.)

(2) Increment j, k: Set reg[j] ← (reg[k] + 1) mod 2^{w}.

(3) Add j, k, ℓ: Set reg[j] ← (reg[k] + reg[ℓ]) mod 2^{w}.

(4) Subtract j, k, ℓ: Set reg[j] ← (reg[k] − reg[ℓ]) mod 2^{w}.

(5) Shift left j, k, ℓ: Set reg[j] ← (reg[k]2^{reg[ℓ]}) mod 2^{w}.

(6) Shift right j, k, ℓ: Set reg[j] ← ⌊reg[k]/2^{reg[ℓ]}⌋.

(7) Word size j: Set reg[j] ← w.

(8) Identify j: Set reg[j] ← this processor’s number.

(9) Read j, k: Set reg[j] ← mem[reg[k]]. (The algorithm fails if reg[k] ≥ s.)

(10) Write j, k: Set mem[reg[k]] ← reg[j].

(11) Jump j, k, ℓ: If reg[k] ≥ reg[ℓ], go to instruction j, rather than proceeding to the next instruction as usual.

A n-bit input to the algorithm is placed into memory at time 0, with n in the first word of memory, and the n bits packed into the next ⌈n/w⌉ words of memory. (The algorithm fails if 2^{w} ≤ n, or if s < 1+⌈n/w⌉.) All other words of memory, registers, etc. are set to 0. Each processor performs one instruction at time 1; each processor performs one instruction at time 2; and so on until all the processors have halted. The output is then in memory, encoded in the same way as the input.

Memory conflicts. During one time step, a single memory location might be accessed by more than one processor.

The concurrent-read exclusive-write PRAM, or CREW PRAM, allows any number of processors to simultaneously read the same memory location, but it does not allow two processors to simultaneously write to a single memory location.

The concurrent-read common-concurrent-write PRAM, or Common CRCW PRAM, allows any number of processors to simultaneously read the same memory location, and allows any number of processors to simultaneously write the same memory location, if they all write the same value. (If two processors attempt to write different values, the algorithm fails.)

A memory location cannot be read and written simultaneously.

Asymptotics. When we say that a parallel algorithm handles an n-bit input using (for example) time O(n) on O(n^{3}) processors with a word size of O(lg n) using O(n^{4}) memory locations, we mean that it works on any sufficiently large computer: there are functions p(n) ∈ O(n^{3}), s(n) ∈ O(n^{4}), w(n) ∈ O(lg n), and t(n) ∈ O(n) such that, for every n, every p ≥ p(n), every s ≥ s(n), every w ≥ w(n), every r larger than the highest register number used in the algorithm, every i larger than the length of the algorithm, and every input string of length n, the algorithm runs without failure on that input in time at most t(n) on a parallel computer of size p, s, w, r, i.

When we do not mention the word size (and number of memory locations), we always mean that the required word size is logarithmic in the time-processor product (and, consequently, the required number of memory locations is polynomial in the time-processor product).

7. Modular exponentiation on the CREW PRAM

In this section we present three successively better parallel algorithms, using a polynomial number of processors, for modular exponentiation on the CREW PRAM.

Time O(n(lg n)^{2}). Given the binary representations of n-bit integers x and y, one can compute xy in time O(lg n) using n^{O(1)} processors. One can compute x/y and x mod y by Newton’s method in time O((lg n)^{2}) using n^{O(1)} processors. See [32, Theorem 12.1].

The base-2 exponentiation algorithm shown in §5, using these subroutines, computes x^{e} mod m in time O(n(lg n)^{2}) if x, e, m have n bits. There are O(n) iterations in the algorithm, each iteration involving a constant number of multiplications and divisions of O(n)-bit integers. The number of processors is polynomial in n.

Time O(n lg n). A division algorithm of Beame, Cook, and Hoover takes time only O(lg n) using n^{O(1)} processors, after a precomputation taking time O((lg n)^{2}). See [5].

The main subroutine in the Beame-Cook-Hoover algorithm computes powers x^{v}, where x has n bits and v ∈ {0, 1, …, n}, in time O(lg n). The idea is to use the CRT to recover x^{v} from the remainders x^{v} mod p for enough small primes p. Beame, Cook, and Hoover precompute the primes p and the powers u^{v} mod p for all u ∈ {0, 1, …, p − 1} and all v ∈ {0, 1, …, n}; then they can compute x^{v} mod p as (x mod p)^{v} mod p. This is one of the ideas that we use in §8.

Time O(n). To save another lg n factor, we choose an integer b ≥ 2 so that lg b ∈ Θ(lg n), and we use the base-b exponentiation algorithm.

The number of iterations drops to O(log_{b} e) = O(n/lg n). Each iteration involves computing y^{b}x^{v} mod m for some v ∈ {0, 1, …, b − 1}; we compute y^{b} and x^{v} in time O(lg n) by the Beame-Cook-Hoover subroutine, then multiply, then use the Beame-Cook-Hoover algorithm again to divide by m.

8. Modular exponentiation on the Common CRCW PRAM

In this section, we present our sublinear-time algorithm, which uses a polynomial number of processors, for modular exponentiation on the Common CRCW PRAM. We begin by informally explaining how we save an additional factor of lg lg n from the running time of our linear-time algorithm from §7. Next we review some FFTbased subroutines that we use in the algorithm. We then present our algorithm along with a proof of its complexity, followed by some discussion.

How we save a lg lg n factor. The algorithm in this section, like the simplified linear-time algorithm of §7, performs several multiplications before each reduction modulo m. It achieves better parallelization than the algorithm of §7 as follows:

(1) It avoids the binary representation in the main loop. It uses the explicit Chinese remainder theorem modulo m modulo p_{j}, as described in §4, to work consistently in the residue representation.

(2) It uses the Cole-Vishkin parallel-addition algorithm to add O(n) integers, each having O(lg n) bits, in time O((lg n)/lg lg n) using n^{1+o(1)} processors. See [12] and [34]. This is where we are making use of the increased power of the CRCW PRAM over the CREW PRAM. Computing such sums is the bottleneck in using the explicit Chinese remainder theorem modulo m modulo p_{j}.

Perhaps the same lg lg n speedup can be obtained with other redundant representations of integers modulo m; we leave this exploration to the reader.

FFT-based subroutines. We take advantage of several standard FFT-based tools here:

(1) Given n bits representing two integers x, y (in binary), one can use the Schönhage-Strassen multiplication algorithm to compute the product xy in time (lg n)^{O(1)} using n^{1+o(1)} processors. See [28].

(2) Given n bits representing two integers x, y with y ≠ 0, one can compute the quotient ⌊x/y⌋ in time (lg n)^{O(1)} using n^{1+o(1)} processors. See [32, Theorem 12.1].

(3) Given n bits representing s integers, one can compute the product of all the integers in time (lg n)^{O(1)} using n^{1+o(1)} processors, by multiplying pairs in parallel.

(4) Given n bits representing integers u, p_{1}, p_{2}, …, p_{s}, one can use the Borodin-Moenck remainder-tree algorithm to compute u mod p_{1}, …, u mod p_{s} in time (lg n)^{O(1)} using n^{1+o(1)} processors. See [21], [10, Sections 4–6], and [8, Section 18].

Without the FFT-based tools, step 4 below requires roughly n^{3} processors to be carried out in polylogarithmic time. This is (after some serialization) adequate for Theorem 8.1, but it becomes a bottleneck when e is much shorter than m.

The algorithm. Fix ϵ > 0. Some of the O constants below depend on ϵ; we indicate when this happens by adding ϵ as a subscript.

Theorem 8.1. There is an algorithm for the Common CRCW PRAM that, given the binary representations of positive integers x, e, and m, of total bit length n, computes the binary representation of x^{e} mod m in time O_{ϵ}(n/lg lg n) using O(n^{2+ϵ}) processors.

Outside the inner loop (step 6 below), the algorithm takes time only (lg n)^{O(1)} using O(n^{2+ϵ}) processors.

Step 1: build multiplication tables. Fix a positive rational number δ ≤ ϵ/3. Compute an integer a ≥ 1 within 1 of δ lg n. Note that an a-bit integer fits into O_{ϵ}(1) words, since a ∈ O_{ϵ}(lg n).

Compute xy and (for y ≠ 0) ⌊x/y⌋ for each pair x, y of a-bit integers in parallel. Store the results in a table of O(2^{2a}) words. This takes time O(a) using 2^{2a} processors; i.e., time O_{ϵ}(lg n) using O(n^{2δ}) processors.

Note that, if x and y are O(lg n)-bit integers, then x and y are also O_{ϵ}(a)-bit integers, so one processor can compute x + y, x − y, xy, and ⌊x/y⌋ in time O_{ϵ}(1) with the help of this table.

Step 2: find primes. Define b = 2^{a}. Note that b ∈ Θ(n^{δ}).

Find the smallest integer k ≥ 6 such that 2^{k} ≥ 4 + 4b(n + k). Use a parallel sieve, as described in [30], to find the primes p_{1}, p_{2}, …, p_{s} between 2 and 2^{k}. This takes time O((1 + δ) lg n) using O(n^{1+δ}) processors; note that s ≤ p_{s} ∈ O(n^{1+δ}).

Define P = p_{1}p_{2} ⋯ p_{s}. Note that lg P ∈ O(n^{1+δ}) by, e.g., [27, Theorem 9]. Note also that P ≥ 4(m^{2} ∑ p_{i})^{b}. (Indeed, 2^{k} ≥ 41, so log P = log p_{1} + log p_{2} + ⋯ + log p_{s} ≥ 2^{k}(1 − 1/ log 2^{k}) by [27, Theorem 4]; also p_{1} + p_{2} + ⋯ + p_{s} ≤ 2^{2k}. Thus lg P ≥ log P ≥ 2^{k−1} ≥ 2 + b(2n + 2k) ≥ lg 4 + b(lg m^{2} + lg ∑ p_{i}).)

Multiply p_{1}, …, p_{s} to compute P. This takes time (lg lg P)^{O(1)} using (lg P)^{1+o(1)} processors; i.e., time (lg n)^{O(1)} using n^{1+δ+o(1)} processors.

Step 3: build power tables. For each i ∈ {1, 2, …, s} in parallel, for each u ∈ {0, 1, …, p_{i} − 1} in parallel, for each v ∈ {0, 1, …, b} in parallel (or serial), compute u^{v} mod p_{i}. Store the results in a table. This takes time O_{ϵ}(lg b) using O(sp_{s}(b + 1)) processors; i.e., time O_{ϵ}(lg n) using O(n^{2+2δ} lg n) processors.

(An alternative approach, with smaller tables, is to find a generator g_{i} for the multiplicative group (Z/p_{i})^{∗}, build a table of powers of g_{i}, and build a table of discrete logarithms base g_{i}. See [30].)

Step 4: compute ECRT coefficients. For each i in parallel, compute P_{i} = P/p_{i}; P_{i} mod p_{i}; P_{i} mod m; and P_{i} mod m mod p_{1}, P_{i} mod m mod p_{2}, …, P_{i} mod m mod p_{s}. This takes time (lg lg P)^{O(1)} using s(lg P)^{1+o(1)} processors; i.e., time (lg n)^{O(1)} using n^{2+2δ+o(1)} processors.

Next, for each i in parallel, compute q_{i} = (P_{i} mod p_{i})^{p_{i}−2} mod p_{i}, so that q_{i} is the inverse of P_{i} modulo p_{i}. This takes time O(lg p_{s}) using s processors; i.e., time O_{ϵ}(lg n) using O(n^{1+δ}) processors.

Step 5: convert to the residue representation. Set x ← x mod m. This takes time (lg n)^{O(1)} using n^{1+o(1)} processors.

Set x_{1} ← x mod p_{1}, x_{2} ← x mod p_{2}, …, x_{s} ← x mod p_{s}. This takes time (lg n)^{O(1)} using n^{1+δ+o(1)} processors. Also set y_{1} ← 1, y_{2} ← 1, …, y_{s} ← 1.

Let l be the number of base-b digits in e. Write e = ∑_{k=0}^{l−1} e_{k}b^{k} with 0 ≤ e_{k} < b. Note that l ∈ O_{ϵ}(n/lg n).

Step 6: inner loop. Repeat the following substeps for k ← l − 1, k ← l − 2, …, k ← 0. Each substep will take time O_{ϵ}((lg n)/lg lg n), so the total time here is O_{ϵ}(n/lg lg n).

Invariant: The integer y represented by (y_{1}, y_{2}, …, y_{s}) is between −m ∑_{i} p_{i} and m ∑_{i} p_{i}, and is congruent to x^{e_{l−1}b^{l−k−2}+···+e_{k+2}b+e_{k+1}} modulo m. The plan is to compute the residue representation of u = y^{b}x^{e_{k}}, and then change y to the integer v identified in Theorem 3.1. Note that |u| < (m^{2} ∑_{i} p_{i})^{b} ≤ P/4; hence v ≡ u ≡ x^{e_{l−1}b^{l−k−1}+···+e_{k+2}b^{2}+e_{k+1}b+e_{k}} (mod m).

For each i in parallel, compute t_{i} ← (y_{i}^{b} mod p_{i})(x_{i}^{e_{k}} mod p_{i})q_{i} mod p_{i}. This takes time O(1) using s processors, thanks to the precomputed power table.

Compute r ← round ∑_{i} t_{i}/p_{i} with the help of Theorem 2.2. The fixed-point divisions take time O(1) using s processors. The addition takes time O((lg s)/lg lg s) using s^{1+o(1)} processors using the Cole-Vishkin addition algorithm.

For each j in parallel, compute y_{j} ← ∑_{i} t_{i}(P_{i} mod m mod p_{j}). This takes time O((lg s)/lg lg s) using s^{2+o(1)} processors, again using the Cole-Vishkin algorithm; i.e., time O((lg n)/lg lg n) using n^{2+2δ+o(1)} processors.

For each j in parallel, compute y_{j} ← (y_{j} − r(P mod m mod p_{j})) mod p_{j}. This takes time O(1) using s processors.

Step 7: convert to the binary representation. Use Theorem 2.1 to compute the binary representation of the integer y whose residue representation is (y_{1}, y_{2}, …, y_{s}). This takes time (lg n)^{O(1)} using n^{2+2δ+o(1)} processors.

Now x^{e} = x^{e_{l−1}b^{l−1}+···+e_{1}b+e_{0}} ≡ y (mod m). The output of the algorithm is y mod m.

Discussion. The algorithm can be converted into a polynomial-size unboundedfan-in circuit of depth O_{ϵ}(n/lg lg n), using the techniques explained in [31]. The last algorithm discussed in §7 can be converted into a polynomial-size bounded-fan-in circuit of depth O(n) by the same techniques.

The only previous sublinear-time algorithms were algorithms that use many more processors or that drastically restrict m. See [16, Section 6] for a survey.

We do not know an exponentiation algorithm that takes time O(n/lg lg n) using n^{2+o(1)} processors. Allowing ϵ to vary with n would hurt the time bound: for example, choosing ϵ as 1/lg lg lg lg n would produce an algorithm that takes time O(n lg lg lg lg n/lg lg n) using O(n^{2+1/lg lg lg lg n}) processors.

In a few applications, the exponent e is only a small part of the input, perhaps \sqrt{n} or \sqrt[3]{n} bits. The running time depends primarily on the length of e: the algorithm takes time O_{ϵ}((lg e)/lg lg n) + (lg n)^{O(1)} using O(n^{2+ϵ}) processors.

References

[1] —, Proceedings of the 20th annual ACM symposium on theory of computing, Association for Computing Machinery, New York, 1988.

[2] Leonard M. Adleman, Kireeti Kompella, Using smoothness to achieve parallelism, in [1] (1988), 528–538.

[3] Amod Agashe, Kristin Lauter, Ramarathnam Venkatesan, Constructing elliptic curves with a known number of points over a prime field, in [33] (2004), 1–17. MR 2005m:11112. URL:

    http://research.microsoft.com/~klauter/.

[4] Manindra Agrawal, Neeraj Kayal, Nitin Saxena, PRIMES is in P, Annals of Mathematics

    160 (2004), 781–793. URL: http://ProjectEuclid.org/Dienst/UI/1.0/Summarize/euclid.
    annm/1111770735.

[5] Paul W. Beame, Stephen A. Cook, H. James Hoover, Log depth circuits for division and related problems, SIAM Journal on Computing 15 (1986), 994–1003. ISSN 0097–5397. MR 88b:68059.

[6] Daniel J. Bernstein, Detecting perfect powers in essentially linear time, and other studies in computational number theory, Ph.D. thesis, University of California at Berkeley, 1995.

[7] Daniel J. Bernstein, Multidigit modular multiplication with the explicit Chinese remainder

theorem, in [6] (1995). URL: http://cr.yp.to/papers.html#mmecrt.

[8] Daniel J. Bernstein, Fast multiplication and its applications, to appear. URL: http://cr.yp.

to/papers.html#multapps. ID 8758803e61822d485d54251b27b1a20d.

[9] Daniel J. Bernstein, Pippenger’s exponentiation algorithm, to be incorporated into author’s

High-speed cryptography book. URL: http://cr.yp.to/papers.html#pippenger.

[10] Allan Borodin, Robert T. Moenck, Fast modular transforms, Journal of Computer and System Sciences 8 (1974), 366–386; older version, not a subset, in [21]. ISSN 0022–0000. MR 51:7365.

URL: http://cr.yp.to/bib/entries.html#1974/borodin.

[11] Henri Cohen, A course in computational algebraic number theory, Graduate Texts in Mathematics, 138, Springer-Verlag, Berlin, 1993. ISBN 3–5440–55640–0. MR 94i:11105.

[12] Richard Cole, Uzi Vishkin, Faster optimal parallel prefix sums and list ranking, Information and Computation 81 (1989), 334–352. ISSN 0890–5401. MR 90k:68056.

[13] Jean-Marc Couveignes, Computing a square root for the number field sieve, in [19] (1993), 95–102. MR 13212222.

[14] Shawna Meyer Eikenberry, Jonathan P. Sorenson, Efficient algorithms for computing the Jacobi symbol, Journal of Symbolic Computation 26 (1998), 509–523. ISSN 0747–7171. MR 99h:11146.

[15] Faith E. Fich, The complexity of computation on the parallel random access machine, in [26] (1993), 843–899.

[16] Daniel M. Gordon, A survey of fast exponentiation methods, Journal of Algorithms 27 (1998), 129–146. ISSN 0196–6774. MR 99g:94014. URL: http://www.ccrwest.org/gordon/dan.html.

[17] Richard M. Karp (chairman), 13th annual symposium on switching and automata theory, IEEE Computer Society, Northridge, 1972.

[18] Donald E. Knuth, The art of computer programming, volume 2: seminumerical algorithms, 3rd edition, Addison-Wesley, Reading, 1997. ISBN 0–201–89684–2.

[19] Arjen K. Lenstra, Hendrik W. Lenstra, Jr. (editors), The development of the number field sieve, Lecture Notes in Mathematics, 1554, Springer-Verlag, Berlin, 1993. ISBN 3–540–57013– 6. MR 96m:11116.

[20] Alfred J. Menezes, Paul C. van Oorschot, Scott A. Vanstone, Handbook of applied cryptography, CRC Press, Boca Raton, Florida, 1996. ISBN 0–8493–8523–7. MR 99g:94015. URL:

http://cacr.math.uwaterloo.ca/hac.

[21] Robert T. Moenck, Allan Borodin, Fast modular transforms via division, in [17] (1972), 90– 96; newer version, not a superset, in [10]. URL: http://cr.yp.to/bib/entries.html#1972/

moenck.

[22] Peter L. Montgomery, An FFT extension of the ellpitic curve method of factorization, Ph.D. thesis, University of California at Los Angeles, 1992. URL: http://cr.yp.to/bib/entries.

html#1992/montgomery.

[23] Peter L. Montgomery, Robert D. Silverman, An FFT extension to the P − 1 factoring algorithm, Mathematics of Computation 54 (1990), 839–854. ISSN 0025–5718. MR 90j:11142.

[24] Andrew M. Odlyzko, Gary Walsh, Hugh Williams (editors), Conference on the mathematics of public key cryptography: the Fields Institute for Research in the Mathematical Sciences, Toronto, Ontario, June 12–17, 1999, book of preprints distributed at the conference, 1999.

[25] Christos M. Papadimitriou, Computational complexity, Addison-Wesley, Reading, Massachusetts, 1994. ISBN 0201530821. MR 95f:68082.

[26] John H. Reif (editor), Synthesis of parallel algorithms, Morgan Kaufman, San Mateo, California, 1993. ISBN 1–55860–135–X. MR 94c:68086.

[27] J. Barkley Rosser, Lowell Schoenfeld, Approximate formulas for some functions of prime numbers, Illinois Journal of Mathematics 6 (1962), 64–94. ISSN 0019–2082. MR 25:1139.

URL: http://cr.yp.to/bib/entries.html#1962/rosser.

[28] Arnold Schönhage, Volker Strassen, Schnelle Multiplikation großer Zahlen, Computing 7 (1971), 281–292. ISSN 0010–485X. MR 45:1431. URL: http://cr.yp.to/bib/entries.

html#1971/schoenhage-mult.

[29] Jonathan P. Sorenson, A sublinear-time parallel algorithm for integer modular exponentiation, in [24] (1999).

[30] Jonathan P. Sorenson, Ian Parberry, Two fast parallel prime number sieves, Information and Computation 144 (1994), 115–130. ISSN 0890–5401. MR 95h:11097.

[31] Larry Stockmeyer, Uzi Vishkin, Simulation of parallel random access machines by circuits, SIAM Journal on Computing 13 (1984), 409–422. ISSN 0097–5397. MR 85g:68018.

[32] Stephen R. Tate, Newton iteration and integer division, in [26] (1993), 539–572.

[33] Alf van der Poorten, Andreas Stein (editors), High primes and misdemeanours: lectures in honour of the 60th birthday of Hugh Cowie Williams, American Mathematical Society, Providence, 2004. ISBN 0–8218–3353–7. MR 2005b:11003.

[34] Uzi Vishkin, Advanced parallel prefix-sums, list ranking and connectivity, in [26] (1993), 215– 257.

[35] Joachim von zur Gathen, Computing powers in parallel, SIAM Journal on Computing (1987), 930–945. MR 89j:68060.

Department of Mathematics, Statistics, and Computer Science (M/C 249), The University of Illinois at Chicago, Chicago, IL 60607–7045

E-mail address: djb@cr.yp.to

Department of Computer Science and Software Engineering, Butler University, 4600 Sunset Avenue, Indianapolis, IN 46208

   E-mail address: sorenson@butler.edu