A new class of scalable parallel pseudorandom number generators based on PohligHellman exponentiation ciphers
Abstract
Parallel supercomputerbased Monte Carlo applications depend on pseudorandom number generators that produce independent pseudorandom streams across many separate processes. We propose a new scalable class of parallel pseudorandom number generators based on Pohlig–Hellman exponentiation ciphers. The method generates uniformly distributed floating point pseudorandom streams by encrypting simple sequences of integer messages into ciphertexts by exponentiation modulo prime numbers. The advantages of the method are: the method is trivially parallelizable by parameterization with each pseudorandom number generator derived from an independent prime modulus, the method is fully scalable on massively parallel computing clusters due to the large number of primes available for each implementation, the seeding and initialization of the independent streams is simple, the method requires only a few integer multiply–mod operations per pseudorandom number, the state of each instance is defined by only a few integer values, the period of each instance is different, and the method passes a battery of intrastream and interstream correlation tests using up to pseudorandom numbers per test. The 32bit implementation we propose has millions of possible instances, all with periods greater than . A 64bit implementation depends on 128bit arithmetic, but would have more than possible instances and periods greater than .
pacs:
02.70.c, 05.10.a, 05.10.Gg, 05.10.Ln, 05.40.a, 07.05.Tp, 95.75.WxI Introduction
We propose a new class of scalable parallel pseudorandom number generators for use in largescale Monte Carlo and other computer applications that require large numbers of independent streams of pseudorandom numbers. The method we propose is based on Pohlig–Hellman exponentiation ciphers.(1); (2) The method creates a pseudorandom stream by encrypting a simple sequence of short integer plaintext messages into ciphertexts using the transformation
(1) 
where each generator instance is based on an independent prime modulus , and an exponent that is coprime to . Here and throughout, means is the remainder of upon division by , with . The method is based upon elementary number theory.(3); (4); (5) For every prime number , the set of integers forms a finite field that is closed under addition and multiplication modulo . Also, the set of nonzero elements forms a group that is closed under multiplication, and for every element in there exists a unique inverse such that . The pseudorandom number generator algorithm described here cycles through a sequence of messages with selected from using some simple pattern that uniformly samples the set. The exponentiation step equation (1) then gives a sequence of ciphertexts that is uniformly distributed on . Uniformly distributed double precision floating point values on the open real interval are formed with a floating point division: .
Most pseudorandom number generators generate the next pseudorandom integer from either the previous pseudorandom integer in the sequence, or by using two or more pseudorandom integers from earlier in the sequence. This method is different in that the pseudorandom sequence arises from the encryption of a sequence of integer messages. In this way, it is similar to cryptographically secure pseudorandom number generators(2); (6) and floating point pseudorandom number generators(7) based on block ciphers. The quality of the pseudorandom sequence from our method is based on modular exponentiation being a good oneway cryptographic function, meaning that it is computationally difficult to ascertain the message from the ciphertext without knowing both the modulus and exponent.(2); (6); (5); (8); (9) The algorithm we propose here produces excellent, longperiod pseudorandom sequences even for 32bit prime moduli. The period of the 32bit implementation we propose is greater than , and has millions of possible instances. Sixtyfour bit implementations have periods in excess of , with more than possible instances.
Two qualitatively different schemes have been used to create scalable systems of pseudorandom number generators for use on massively parallel supercomputers: stream splitting and parameter splitting.(10) Stream splitting is based on a single pseudorandom number generator with an extremely long period, with parallelization accomplished by subdividing the full period into independent nonoverlapping subsequences. This method requires careful seeding to ensure that no two subsequences will overlap for any feasible set of processes. Parameter splitting uses a single algorithm, but produces independent pseudorandom sequences by assigning different parameters to each process. The SPRNG web site(11); (12) gives examples of parallel generators of both classes. For example, linear congruential generators of the form , with a power of two, can be parallelized by parameter splitting.(13); (11); (12); (14) The parameters are , a welltested multiplier that ensures a maximal period and , an odd number less than . Each process uses the same values of and , but the value of is chosen to be different for each process.(15) Another widely used class of parallel generators are based on the lagged Fibonacci method.(16); (17); (11); (12); (18) They are usually of the form , where is a power of two, is one of the operations addition, exclusive or, subtraction, or multiplication, and are integer parameters chosen based on primitive polynomials modulo 2.(2); (5); (19) Lagged Fibonacci generators have periods of the order of , with typically of the order of several hundred to several thousand. The state of the generator is defined by a table of integers or bits, which represent the most recent pseudorandom values in the sequence. Parallel implementation of these algorithms can be accomplished by either stream splitting or parameter splitting. (11); (12)
The parallel pseudorandom number generator class we propose here is based on parameter splitting. In our algorithm, every independent process is assigned a different prime modulus , which produces an independent pseudorandom sequence. The sequence of pseudorandom numbers produced using equation (1) are effectively uncorrelated within a single stream and between streams, and the period of every stream is different. The number of independent streams is limited only by the number of prime numbers in the range defined by the implementation. For example, in a 32bit implementation the prime moduli can be chosen from the set of 98,182,656 primes in the range .(20)
To demonstrate the relation to encryption and explain some of the useful properties of the method, the message can be decrypted from the ciphertext using the decryption exponent , i.e. for some .(1); (2); (3); (4); (5) The decryption exponent exists and is unique if and are coprime, i.e if . The decryption exponent can be determined quickly using the extended Euclidean algorithm.(2); (3); (4); (5) The decryption is based on Fermat’s little theorem:(3); (4); (5) for any prime n and for all in , . Therefore
(2) 
For any allowed exponent , equation (1) maps onto a permutation of the same set. Therefore, any message sequence that uniformly samples will produce ciphertexts that uniformly sample . The Pohlig–Hellman algorithm is closely related to the more widely used RSA public key encryption method(2); (3); (4); (5); (6); (21) which uses moduli that are products of two large primes. Using a composite modulus creates a serous weakness for this application since every message that is a multiple of or gives a ciphertext that is also a multiple of or , respectively.
In the Pohlig–Hellman cryptographic scheme, the key consisting of the prime and exponent must be known only to the sender and receiver. Otherwise, an eavesdropper can easily determine the decryption exponent. Cryptography theory suggests using safe primes for the moduli, i.e. primes for which is also prime.(1); (9); (2); (22) This helps ensure eavesdroppers will find it computationally difficult to take the discrete logarithm to recover the message from the ciphertext.(1); (3); (4); (5) While our empirical tests do not give noticeably better statistics for safe primes than for primes in general, the number of safe primes does not seriously limit the scaling with 32bit moduli since there are 3,060,794 safe primes in the range .(23)
The algorithm presented here is similar to a block cipher operating in the counter mode which outputs a sequence of pseudorandom numbers that are ciphertexts resulting from a simple sequence of plaintexts.(2); (6) The Pohlig–Hellman method can be used to generate cryptographically secure pseudorandom sequences, but cryptographic applications require primes with hundreds of digits.(1); (2); (3); (4); (5); (6); (8); (9) Our application, instead, uses Pohlig–Hellman encryption of 32bit messages to produce pseudorandom sequences suitable for use in Monte Carlo simulations and other applications. The exponentiations can be accomplished using the method of repeated squaring,(2); (3); (4); (5) so the number of multiply–mod operations needed to calculate is less than . As we will show below, the algorithm can be implemented using small exponents that require very few multiply–mod operations per pseudorandom number.(24) If the modulus is chosen to be a prime number less than , each multiply–mod operation can be performed in hardware on a 64bit processor with one 64bit multiply, and one 64bit mod. Our algorithm leads to a new scalable class of fast and effective pseudorandom generators based on parameter splitting. There are several advantages to this pseudorandom number generation method.

The algorithm is based on elementary number theory and cryptography.

The method is trivially parallelizable by parameterization, with each instance derived from an independent prime modulus.

Pseudorandom sequences that result from different prime moduli are independent, and have different periods.

The method is fast since it requires only a few integer multiply–mod operations per pseudorandom number.

The method is fully scalable on massively parallel computing clusters due to the millions of available 32bit primes.

As we will show below, the period of the generator can be greatly extended by using message skips derived from a prime number linear congruential pseudorandom number generator.

The seeding and initialization of the independent streams is simple, and it is possible to initialize each process without needing information about the states of any of the other processes. The state of each generator is defined by only seven integers: , , , and , and integer pseudorandom skip parameters , and defined below. All of these values can be stored in local memory, and require no global storage. The values , , and are fixed in each instance. Only the three values , , and change during calls to the generator.

The algorithm allows one to quickly jump far forward or backward in the pseudorandom sequence.

The algorithm is simple enough to allow the generator to be implemented as an inline function for efficiency.

The algorithm can be implemented in parallel on vector computers and GPUs.

The method passes a battery of intrastream and interstream correlation tests using up to pseudorandom numbers per test.
Ii Message skipping algorithms
The selection of messages to be encrypted for a given modulus can be accomplished in many different ways, but they can all be expressed in terms of an integer skip sequence with the skips chosen from :
(3a)  
(3b) 
Note that if the skip values are chosen uniformly and randomly (not pseudorandomly) from , then the sequence of messages forms a uniform random sequence on the same set. Message is, in effect, a onetime pad encryption of message .(2) Since the set of ciphertexts is onetoone with the set of messages, the sequence of ciphertexts also forms a uniformly distributed random sequence.
In encryption, it is essential to avoid cribs, i.e. messages that result in easily decoded ciphertexts. For example, the messages are cribs for all allowed exponents since , respectively. Calculation efficiency makes it desirable to use small exponents to reduce the number of multiply–mod operations needed to generate each pseudorandom number. RSAbased cryptographic applications often use small exponents such as .(6) Messages with and result in trivially decodable ciphertexts, so exponents result in additional cribs. Cribs can be avoided by padding messages to ensure no messages are close to or . For cryptographic applications, the padding needs to have a random component.(2); (6) For our application, eliminating cribs from the message stream biases the ciphertext distribution due to the elimination of ciphertexts formed from messages with and , resulting in a slight but systematic undersampling of small and large values of . Even though the effect of this is small, we choose not to implement an algorithm that does not give a uniform distribution of ciphertexts over the full period. For our purpose, it is not necessary to eliminate cribs, since they would appear in a random message sequence, but rather to prevent correlated sequences of cribs. Our goal is to select a simple skip pattern that ensures a uniform sampling of the set of all messages, is computationally fast, has a long period, allows the use of small encryption exponents, and avoids correlated cribs. We will accomplish this by using a pseudorandom skip sequence that, over the full period of the generator, uniformly samples the messages in , but we will first examine the properties of pseudorandom ciphertext sequences derived from simpler skip patterns.

The simplest skip sequence that uniformly samples is the unit skip, i.e. for all . Pohlig–Hellman encryption of this sequence is, in effect, a block cipher operating in counter mode.(2); (6) The message sequence is , with the period of the generator being . For 32bit moduli, one can easily test the full period of the generator. In spite of the cribs near and , the sequence with passes most statistical tests for randomness even for small exponents. Naturally, the fullperiod sequence produces a perfectly uniform onedimensional histogram since after steps every value will appear once, and only once, in the sequence. Except for the onedimensional frequency test, and other tests that are affected by the uniformity of the sampling of numbers in the range , and in spite of the symmetry noted below, the exponentiation cipher with unit skips passes all of our other statistical correlations tests for exponents and . By comparison, all linear congruential generators fail all dimensional correlations tests once a substantial fraction of the period has been exhausted, since all dimensional histograms become uniform as the period of the generator is approached.(19) The exponentiation cipher does not suffer from this; see figures 1 and 2 which display the twodimensional correlation patterns for a prime number linear congruential generator, and a unitskip Pohlig–Hellman generator using the same prime modulus. As with most pseudorandom number generators, there is a symmetry in the pseudorandom sequence. The Pohlig–Hellman cipher has the symmetry
i.e. the ciphertexts derived from messages in are strongly correlated with the messages in , but in reverse order.

The next simplest skip algorithm is the constant skip: with , so . This message pattern is especially important for understanding the quality of our proposed pseudorandom skip algorithm discussed next. For most values of , the constant skip removes sequential and closely spaced cribs. However, other than breaking up closely spaced cribs, the pattern produced by constant skips is not substantially better than the unit skip pattern since
(4) i.e. this is the same as the unitskip sequence except for a single constant factor . Like the unit skip case, there is a symmetry in the pseudorandom sequence since . As with the unit skip, the constant skip passes a battery of fullperiod statistical tests for and , except ones affected by the uniformity of the onedimensional distribution.

Our recommended skip pattern is a pseudorandom skip produced by a prime number linear congruential pseudorandom number generator:(19); (13); (25)
(5) with prime modulus . The multiplier is chosen to be a primitive root(3); (4); (5) that delivers a full period, welltested pseudorandom sequence.(19); (13); (25) This gives pseudorandom skips in the range , with the period of the skip sequence being . This results in the following pseudorandom ciphertext sequence:
(6a) (6b) (6c) Using a pseudorandom skip sequence serves several important purposes:

A pseudorandom skip effectively eliminates the problem of correlated cribs, allowing the use of small exponents.

Using a pseudorandom skip extends the period of the generator to .

The method provides a uniform sampling of ciphertexts over the full period of the generator. Each message in , and hence each ciphertext in , will appear exactly times in the sequence.

The state of each generator is defined by only seven integers: .

Implementations using 32bit primes are fast. With small exponents, each pseudorandom number requires only a few 64bit multiply–mod operations that can be implemented in hardware on 64bit processors.

The implementations suggested below pass a battery of statistical tests with up to pseudorandom numbers per test.
First, let’s prove that the period of the generator is . Since is a primitive root mod , after pseudorandom skips will have cycled through every value in , so and . This means the sequence of messages separated by multiples of steps in the sequence are derived from a constant skip . Since and are coprime to , the constant skip is also coprime to . Therefore, the state of the generator after steps, with and , is given by
(7a) (7b) (7c) This demonstrates the mathematical form of the pseudorandom sequence across the full period of the generator. Since and , and each subsequence of messages of length is different from every other subsequence, the period the generator is . Over the full period, every ciphertext in will appear exactly times. One can use equation (7) to skip steps forward directly to any point in the pseudorandom sequence, with skips being fast if is close to a multiple of . Backward skips of are accomplished by replacing with , subtracting rather than adding in equation (7b), and reordering equations (7a) and (7b).(27)
The simple form of the ciphertext sequence allows one to determine the fullperiod dimensional correlations pattern for rectangular regions with volume in steps, i.e. without needing to exhaust the generator. The dimensional density of points over the full period is so the average number of points in the volume above is . One can select points , with the succeeding points given by , , etc. where the skips take on all values in . A magnified region near the origin of the fullperiod twodimensional pattern for safe prime , , ,(13) and is shown in figure 3, with the fullperiod threedimensional pattern shown in figure 4.
We tested the quality of 32bit pseudorandom sequences based on pseudorandom skip sequences using a battery of independent statistical tests. We first tested uniform double precision floating point pseudorandom sequences over the period of the skip generator for more than tenthousand different safe primes between and , and used a single prime number linear congruential pseudorandom number generator recommended by L’Ecuyer:(13) and . Even using exponents as small as , the intrastream pseudorandom sequences pass all of our statistical tests across the period of the skip generator.
We then used equation (7) to test for intrastream correlations among ciphertexts separated by steps in the sequence across the constant skip period using thousands of safe primes between and . This pseudorandom sequence is given by , with and . These constant skip sequences pass all of our the statistical tests for exponents as small as for up to pseudorandom numbers, and all but the onedimensional tests across the period . (As noted earlier, the onedimensional constant skip distribution becomes uniform over the period. These are the only tests that the constant skip sequences appear to fail.) Based on this, we recommend using exponents or , which require only four or five multiply–mod operations per exponentiation, respectively.
Since ciphertexts separated by steps demonstrate good intrastream statistics because of the pseudorandom skip, and ciphertexts separated by multiples of steps demonstrate good statistics with constant skips , one has good reason to believe that the entire pseudorandom sequence should pass a battery of statistical tests until the length of the pseudorandom sequence approaches the full period. To test this, we performed intrastream statistical tests of the pseudorandom 32bit skip algorithm over as large a fraction of the period and for as many different moduli as possible. We tested sequences of pseudorandom numbers using hundreds of different primes, sequences of pseudorandom numbers using dozens of different primes, and sequences of pseudorandom numbers using a few different primes. This latter test corresponds to several thousand periods of the skip sequence. The method passed every test we applied.
We also tested to ensure that the algorithm displayed lack of correlation between streams. Suppose one has processes each with a different prime modulus with , and pseudorandom sequences . Our interstream correlations tests draw the pseudorandom numbers from the streams in the order . We tested groupings of , , , and different streams, using exponents . We also included tests using safe primes, i.e. all of the safe primes in . To test that seeding coincidences do not create spurious correlations, we performed many of the interstream tests by starting every sequence with exactly the same initial values and , but avoiding cases in which one of the early messages in the sequence is a crib. The interstream correlations passed every test for up to pseudorandom numbers per test.

Iii Tests
We first applied wellestablished pseudorandom number test suites DIEHARD,(28) NIST(29), and TestU01,(30) to ensure the generator passes a wide variety of tests. We then applied the following ten additional tests that produce histograms to which one can apply a chisquare test.(19) For each test, we calculated and the associated value, i.e. the onesided probability of having that value above or below the median. We applied these additional tests to sequences of up to pseudorandom numbers per test.

Onedimensional frequency test:(19) We distributed the sequence of pseudorandom numbers into a onedimensional histogram with bins, and compared the histogram to a uniform Poisson distribution.

Twodimensional serial test:(19) We distributed pairs of pseudorandom numbers into a twodimensional histogram with bins, and compared the histogram to a uniform Poisson distribution. This tests for sequential pair correlations in the sequence.

Threedimensional serial test:(19) We distributed triplets of pseudorandom numbers into a threedimensional histogram with bins, and compared the histogram to a uniform Poisson distribution. This tests for sequential triplet correlations in the sequence.

Fourdimensional serial test:(19) We distributed groups of four pseudorandom numbers into a fourdimensional histogram with bins, and compared the histogram to a uniform Poisson distribution. This tests for sequential fourpoint correlations in the sequence.

Five dimensional serial test:(19) We distributed groups of five pseudorandom numbers into a fivedimensional histogram with bins, and compared the histogram to a uniform Poisson distribution. This tests for sequential fivepoint correlations in the sequence.

Poker test:(19) We used groups of five pseudorandom numbers and counted the number of pairs, threeofakind etc. formed from five cards and ten denominations, and compared the resulting histogram to a Poisson distribution derived from the exact probabilities. This tests for a variety of fivepoint correlations in the sequence.

Collision test:(19) We distributed pseudorandom numbers into bins and compared the histogram of the number of collisions against a Poisson distribution derived from the exact probabilities. This tests for longrange correlations in the sequence.

Runs test:(19) We compared the histogram of the length of runs of 0’s () and 1’s () to the Poisson distribution derived from the exact probabilities. This tests for shortrange correlations of the leading bits.

Fourier test:(11) We used a fast Fourier transform(31) to calculate the Fourier coefficients of sequences of double precision floating point pseudorandom numbers,
(8) where . The real and imaginary parts of the Fourier coefficients with should each be gaussian distributed about zero with variance . We distributed the real and imaginary parts of the Fourier coefficients into a histogram, and compared the result to a Poisson distribution derived from the exact gaussian distribution. This tests for longrange pair correlations in the sequence.

Twodimensional Ising model energy distribution test:(32); (33) We performed a Wolff algorithm(34) Monte Carlo simulation at the critical point of the twodimensional Ising model on a square lattice, and compared the energy histogram to a Poisson distribution derived from the exact probabilities(32); (33) Since the Wolff algorithm is based on stochastically growing fractal critical clusters that can span the system, this tests for longrange correlations in the pseudorandom sequence, and has proven to be effective at identifying weak generators.(32); (35) Assigning an independent generator to each of the 32768 bonds in the lattice tests provides an additional test for exposing interstream correlations.
For several of these tests, we sampled the loworder bits to confirm they do not harbor any hidden correlations. We define passing our tests by there being no events among the tests. Since there were tens of thousands of independent tests, we also counted the number of events in our samples to confirm that the number was consistent with the expected number. Finally, we applied chisquare tests to histograms of the fullperiod distribution of sequences in dimensional rectangular volumes of size with and , for . The generator passed every test.
Iv Implementation
The algorithm for generating uniformly distributed double precision floating point pseudorandom numbers on the open interval is:
(9a)  
(9b)  
(9c)  
(9d)  
(9e) 
In a 32bit implementation, operations (9a)(9c) of the the algorithm can be implemented using 64bit unsigned integer arithmetic, which can be executed in hardware on 64bit processors. This algorithm generates the sequence given by equations (7). For efficiency, one can precalculate the double precision floating point value and perform a floating point multiplication instead of a floating point divide in step (9d).
For a multiprocessor environment, each of the processes can be assigned an independent prime modulus using wellestablished primality tests. The RabinMiller test,(5); (36); (37) which is the same as Algorithm P in Knuth,(19) provides a simple probabilistic test for primality. Every odd prime with odd satisfies one of the following conditions for every base in : either , or for some some in the range . A composite number satisfying these criteria is called a strong pseudoprime to base . For any odd composite, the number of bases for which is a strong pseudoprime to the base is less than , so if the test is applied repeatedly with randomly chosen bases in , the probability that a composite will pass every test is less than .(5); (19); (36); (37) Better yet, the RabinMiller test can quickly and deterministically identify all primes below . There are no composite numbers below that are strong pseudoprimes to all of the twelve smallest prime bases ().(38); (39); (40); (41); (42) Therefore, any number less than that passes the RabinMiller test for all twelve of these bases is prime. Likewise, any number less than that passes the RabinMiller test for all of the five smallest prime bases () is prime. For efficiency, one first checks to see if any small primes divide before applying the RabinMiller test.
One can initialize the generators using the system time variable and the process identifier to construct unique values of , and different starting values and for each of the processes. One might use equation (7) to skip each generator forward relative to some base state such as and to ensure the initial skips and messages are all unsynchronized. Even if two generators later become synchronized with the same values of and , the values for the ciphertexts for different moduli will be different due to the encryption step, and the message synchrony is removed after a few skips, even if the skips were to remain synchronized.
The implementation for 32bit moduli is simple and fast since the multiply–mod operations can be executed in hardware using standard 64bit unsigned arithmetic on 64bit processors. The quality of the pseudorandom sequences does not appear to be dependent on the values of the prime moduli, but we recommend using safe primes selected from unless the number of processes exceeds 3,060,794, the number of safe primes in that range. If the implementation requires more than 3,060,794 instances, there are 49,091,941 primes in with coprime to , and 92,045,560 primes in that range with coprime to . The number of welltested 31bit prime number pseudorandom skip generators is not scalable, but one can use the ones available in the literature(25); (13) to substantially extend the number of possible instances.(43)
The full period of each generator is . This 32bit implementation passes all of our intrastream and interstream correlations tests for and , for up to pseudorandom numbers per prime modulus. The exponent requires only five multiply–mod operations per pseudorandom number, one for the skip and four for the exponentiation, and the exponent requires only six. Consequently, the code is compact and simple enough to be implemented as an inline function.
If far longer periods or far more instances are needed, one can implement the algorithm using 64bit primes. There are approximately safe primes in . L’Ecuyer(13) provides suitable pseudorandom skip parameters near so the period of each process is . Using 64bit prime moduli results in a speed penalty with current processors since the multiply–mod operations need to be implemented in using 128bit unsigned integers. Since a single process is unlikely to exhaust a 63bit pseudorandom skip, one might consider using smaller encryption exponents for efficiency. Preliminary statistical tests indicate the method works well for 64bit safe prime moduli even when using , the smallest allowed exponent.
V Conclusion
We propose a new class of parallel pseudorandom number generators based on Pohlig–Hellman exponentiation ciphers. The method creates pseudorandom streams by encrypting simple sequences of integer messages. The method is fully scalable based on parametrization since each process can be assigned a unique prime modulus. By using pseudorandom skips among messages, one can use small exponents and the period is greatly extended. For 32bit implementations, only a few 64bit multiply–mod operations are needed per pseudorandom number. There are millions of possible independent instances, and the period of each instance is is greater than . We have tested thousands of different pseudorandom streams for intrastream and interstream correlations using up to pseudorandom numbers per test, and all pass a battery of statistical tests. A 64bit implementation would have more than possible instances and periods greater than . Sample C++ code for a 32bit multiprocessor MPI implementation of the Pohlig–Hellman pseudorandom number generator with pseudorandom skip can be found at http://works.bepress.com/paul_beale/.
Vi Acknowledgements
We thank Matt Glaser, Rudy Horne, Nick Mousouris, Ethan Neil, Robert Blackwell, John Black, and David Grant for helpful discussions. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research.
References
 S. Pohlig and M. Hellman, IEEE Trans. Inform. Theory (24): 106 (1978).
 B. Schneier, Applied Cryptography (Wiley, New York, 1994).
 T. Koshy, Elementary Number Theory and Applications (Academic, San Diego, 2002).
 J.H. Silverman, A Friendly Introduction to Number Theory (Pearson, New York, 2006).
 N. Koblitz, A Course in Number Theory and Cryptography (SpringerVerlag, New York, 1987).
 N. Ferguson, B. Schneier, T. Kohno, Cryptography Engineering: Design Principles and Applications, (Wiley, Indianapolis, 2010).
 W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, (Cambridge, New York, 1992), 2nd ed.
 M. Blum and S. Micali, SIAM J. Comput. 13, 850, (1984).
 S. Patel and G.S. Sundaram, CRYPTOÕ98, LNCS 1462, 304, (1998).
 H. Bauke and S. Mertens, Phys. Rev. E 75, 066701 (2007).
 M. Mascagni and A. Srinivasan, ACM Trans. Math. Software 26, 436 (2000).
 The Scalable Parallel Random Number Generators Library (SPRNG), http://www.sprng.org.
 P. L’Ecuyer, Math. Comp. 68, 249 (1999).
 M. Mascagni, Parallel Comput. 24, 923 (1998).
 The period of every such generator is , and the low order bits are highly correlated within each stream, as well as between streams.
 M. Mascagni, M. L. Robinson, D. V. Pryor and S. A. Cuccaro, Lec. Notes Statistics 106, 263 (1995).
 M. Mascagni, S. A. Cuccaro, D. V. Pryor and M. L. Robinson, J. Comp. Phys. 119, 211 (1995).
 M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput. Sim. 8(1), 3 (1998).
 D. Knuth, The Art of Computer Programming, vol. 2 (AddisonWesley, Reading, Massachusetts, 1999).
 The OnLine Encyclopedia of Integer Sequences, http://oeis.org/A036378.
 R. Rivest, A. Shamir, and L. Adelman, Commun. ACM 21, 120 (1978).
 This makes a Sophie Germain prime.
 The OnLine Encyclopedia of Integer Sequences, http://oeis.org/A211395; http://oeis.org/A211397.
 A cryptographically secure implementation of the Pohlig–Hellman cipher would ordinarily be implemented with randomly chosen exponents of the same order as , since that would force an eavesdropper to resort to solving the discrete logarithm problem, even if she knew .(3); (4); (5)
 P. L’Ecuyer, Commun. ACM 31, 741 (1988).
 C. Bays and S.D. Durham, ACM Trans. Math. Software 2, 59 (1976).
 The mathematical form of the message sequence can be used to demonstrate why the method requires . If , one can evaluate the sum in equation (6b) in closed form to give . Therefore, using Fermat’s little theorem, and , so period is reduced to .
 G. Marsaglia, DIEHARD: a battery of tests of randomness (1996); see http://stat.fsu.edu/pub/diehard/.
 A. Rukhin, J. Soto, J. Nechvatal, M. Smid, E. Barker, S. Leigh, M. Levenson, M. Vangel, D. Banks, A. Heckert, J. Dray, and S. Vo, NIST special publication 80022, National Institute of Standards and Technology (NIST), Gaithersburg, Maryland, USA, 2001; see http://csrc.nist.gov/rng/.
 P. L’Ecuyer and R. Simard, ACM Trans. Math. Software 33, 22 (2007); see http://www.iro.umontreal.ca/~simardr/testu01/tu01.html.
 J.W. Cooley and J.W. Tukey, Math. Comput. 19, 297 (1965).
 P.D. Beale, Phys. Rev. Lett. 76, 78 (1996).
 R.K. Pathria and P.D. Beale, Statistical Mechanics (Academic, Boston, 2011), 3rd ed.
 U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
 A.M. Ferrenberg, D.P. Landau and Y.J. Wong, Phys. Rev. Lett. 69, 3382 (1992).
 G.L. Miller, J. Comp. Sys. Sci. 13, 300 (1976).
 M.O. Rabin, J. Number Theory 12, 128 (1980).
 C. Pomerance, J. L. Selfridge and S. S. Wagstaff, Jr., Math. Comp. 35, 1003 (1980).
 The OnLine Encyclopedia of Integer Sequences http://oeis.org/A014233.
 Zhenxiang Zhang, Math. Comp. 70, 863 (2001).
 Y. Jiang and Y. Deng, http://arxiv.org/abs/1207.0063v1.
 See http://mathworld.wolfram.com/RabinMillerStrongPseudoprimeTest.html.
 One might consider using poweroftwo linear congruential generators for the skip generator since they form a scalable class,(11); (12); (13) but we have not investigated the effects the strong correlations in the loworder bits of the skips might have on correlations of the ciphertexts.