# Parallel pseudorandom number generators

This question is primarily related to a practical software-engineering problem, but I would be curious to hear if theoreticians could provide more insight in it.

Put simply, I have a Monte Carlo simulation that uses a pseudorandom number generator, and I would like to parallelise it so that there are 1000 computers running the same simulation in parallel. Therefore I need 1000 independent streams of pseudorandom numbers.

Can we have 1000 parallel streams with the following properties? Here $X$ should be a very well-known and widely-studied PRNG with all kinds of nice theoretical and empirical properties.

1. The streams are provably as good as what I would get if I simply used $X$ and split the stream generated by $X$ into 1000 streams.

2. Generating the next number in any stream is (almost) as fast as generating the next number with $X$.

Of course if we simply used $X$, always discarding 999 numbers and picking 1, then we certainly would have property 1, but we would lose in the running time by factor 1000.

A simple idea would be to use 1000 copies of $X$, with seeds 1, 2, ..., 1000. This certainly would be fast, but it is not obvious if the streams have good statistical properties.

After some Googling, I have found, for example, the following:

• The SPRNG library seems to be designed for exactly this purpose, and it supports multiple PRNGs.

• Mersenne twister seems to be a popular PRNG nowadays, and I found some references to a variant that is able to produce multiple streams in parallel.

But all this is so far from my own research areas, that I couldn't figure out what is really the state-of-the-art, and which constructions work well not only in theory but also in practice.

Some clarifications: I do not need any kind of cryptographic properties; this is for scientific computation. I will need billions of random numbers, so we can forget any generator with a period of $< 2^{32}$.

Edit: I cannot use a true RNG; I need a deterministic PRNG. Firstly, it helps a lot with debugging and makes everything repeatable. Secondly, it allows me to do, e.g., median-finding very efficiently by exploiting the fact that I can use the multi-pass model (see this question).

Edit 2: There is a closely related question @ StackOverflow: Pseudo-random number generator for cluster environment.

• why wouldn't you use the PRNG with $1000$ independently sampled seeds? i don't understand how this does not satisfy 1 and 2, since you require no coordination between the different machines Jun 7 '11 at 14:20
• I'm not an expert, but recently (searching information about a TCS question) I found this hardware: idquantique.com/true-random-number-generator/… ... a PCI board that can generate a 16Mbits/sec stream of (quantum) random bits. ... you can buy a bunch of them and implement a few random number generator servers ... not a great theoretical approach but the bits are guaranteed to be "good" :-) :-) Jun 7 '11 at 14:26
• @Vor: I would like to keep everything repeatable and deterministic. Given a fixed seed, I want to get exactly the same result if I re-run the experiment. And I want to be able to run the same experiment on a single machine and again get the same results. (For one, it helps a lot when debugging parallel algorithms...) Jun 7 '11 at 14:32
• @Jukka: ok! ... and I suppose that storing billions of unzippable wild bits along with the experiment results is not so feasible :-) ... a PRNG expert is needed! Jun 7 '11 at 14:44
• Thanks for the answers so far! Let's see if we get more participation with a bounty... Jun 11 '11 at 22:06

You can use an evolution of the Mersenne Twister algorithm developed by Saito and Matsumoto:

SIMD-oriented Fast Mersenne Twister (SFMT)

SFMT is a Linear Feedbacked Shift Register (LFSR) generator that generates a 128-bit pseudorandom integer at one step. SFMT is designed with recent parallelism of modern CPUs, such as multi-stage pipelining and SIMD (e.g. 128-bit integer) instructions. It supports 32-bit and 64-bit integers, as well as double precision floating point as output. SFMT is much faster than MT, in most platforms. Not only the speed, but also the dimensions of equidistributions at v-bit precision are improved. In addition, recovery from 0-excess initial state is much faster. See Master's Thesis of Mutsuo Saito for detail.

The period varies from $2^{607}-1$ to $2^{216091}-1$.

Using one same pesudorandom number generator for generating multiple independent streams by changing the initial values may cause a problem (with negligibly small probability). To avoid the problem, using different parameters for each generation is preferred. This technique is called dynamic creation of the MT parameters.

In the SFMT source code you can find some examples of parameter sets (of variable periods) and an awk script to convert a CSV file to a compilable parameter set. There is also a tool called "Dynamic Creation of Mersenne Twister generators".

The authors recently developed another modified version of the Mersenne Twister - Mersenne Twister for Graphic Processors - designed to run in GPUs and take advantage of their native parallel execution threads. The key feature is speed: $5 \times 10^7$ random integers every 4.6ms on a GeForce GTX 260.

The periods of generated sequence are $2^{11213}-1$ , $2^{23209}-1$ and $2^{44497}-1$ for 32-bit version, and $2^{23209}-1$, $2^{44497}-1$, $2^{110503}-1$ for 64-bit version. It It support 128 parameter sets for each period, in other words, it can generate 128 independent pseudorandom number sequences for each period. We have developed Dynamic Creator for MTGP, which generates more parameter sets

Indeed they provide a MTGPDC tool to create up to $2^{32}$ parameter sets (i.e. independent streams).

The algorithm passes the main randomness tests like Diehard and NIST.

A preliminary paper is also availbale on arXiv: A Variant of Mersenne Twister Suitable for Graphic Processors

• A related but older tool is Matsumoto and Nishimura (1998): Dynamic Creation of Pseudorandom Number Generators. But I haven't been able to figure out which of these tools are just a proof of concept and which are widely-used industry-strength software packages. Jun 13 '11 at 20:52
• @Jukka: perhaps you can ask it directly to the authors of the MTGP algorithm. From their site: "... Any feedback is welcome (send an email to Mutsuo Saito, saito "at sign" math.sci.hiroshima-u.ac.jp and m-mat "at sign" math.sci.hiroshima-u.ac.jp) ...". Perhaps they may not be 100% impartial, but they surely know well the strong and weak points of MTGP, and can tell you if it can be suitable for your purpouses. Jun 13 '11 at 22:13
• It seems that Mersenne Twister + Dynamic Creation is the recommended way to do it in Mathematica. Jun 16 '11 at 22:34
• @Jukka: MT+DC package can be found on Matsumoto's site, too (math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html); and I think that MTGP is only a variant suitable for GPUs. So MT+DC seems a better (and tested/stable) choice (unless you absolutely need $5 \times 10^{7}$ random integers every 4.6ms on each stream :-)))) Jun 16 '11 at 23:36
• @Vor: If you edit your answer and replace MTGP with dcmt, I can accept it. Jun 17 '11 at 12:13

There seem to be many ways to tackle this problem, but one simple way would be to use the Blum Blum Shub PRNG. This PRNG is defined by the recurrance relation $x_{i+1} = x_i^2 \mbox{ mod }N$, where $N$ is a semiprime. To get a random bit out of this you can simply take the bit parity of $x_i$. What's nice about this is that since $x_{i+k} = x_i^{2^k}\mbox{ mod }N = x_i^{2^k \mbox{ mod } \lambda(N)}\mbox{mod }N$ you can directly calculate any step in time constant in $k$ (i.e. $O(\log(N)^3)$ or faster depending on which multiplication algorithm you use for the modular exponential). Thus is you have $M$ machines, then for the machine indexed by $y$ you can use the generator $x_{i+1,y} = x_i^{2^M \mbox{mod }\lambda(N)}\mbox{ mod }N$, where $x_{0,y} = x_0^{2^y \mbox{ mod }\lambda(N)}\mbox{ mod }N$, where $x_0$ is your seed. Conveniently this generates exactly the same stream of numbers as if you used a single stream and distributed it's output to each of the machines in turn.

This isn't the fastest of PRNGs, though, so it will only be useful if the overhead of whatever you are doing in the simulation is significantly more than the cost of the PRNG. However it is worth pointing out that it will be much faster for certain combinations of $M$ and $N$ than others, particularly if the binary representation of $2^M \mbox{ mod }\lambda(N)$ contains few 1s or is small.

• I think it would be faster to let each machine generate a contiguous portion of the sequence, spacing them so far apart that they will not intersect. Anyway, using the Blum Blum Shub for non-cryptographic applications seems to me a bit an overkill. Jun 8 '11 at 0:55
• @Antonio: Yes, that would be slightly faster, particularly if you know ahead of time exactly how many trials you need. If you don't know, then I think you'll get the same scaling either way. Wierdly Blum Blum Shub was exactly the PRNG we were thaught in computational physics years ago. If you aren't using it for cryptographic purposes you can use a much smaller modulus, so it's not really that slow, and for many tasks the it will be fast compared to whatever function of the random variable you need to compute. Jun 8 '11 at 1:51

How about a preprocessing phase? Given a random seed $s$ (of size $n$), run $X$ to get a pseudorandom stream of size $1000n$. Denote this stream by $s_1,s_2,\ldots,s_{1000}$, where for $1\le i \le 1000$, $s_i$ is a contiguous portion of the stream of size $n$.

This preprocessing phase can be done with a very low overhead, given the fact that $X$ is an efficient PRNG (today, we have very fast PRNG's).

Now, give $s_i$ as the seed to the $i$th machine, which uses $X$ to generate its own pseudorandom stream.

Given the nice properties of $X$, unless $s$ is known, for any $1 \le i < j \le 1000$, the seeds $s_i$ and $s_j$ are computationally independent. Moreover, you only have to generate and save one small seed (i.e. $s$); therefore, this approach does not need a great deal of true randomness or storage.

• Isn't this essentially the same approach as what @Antonio suggested: use a PRNG to generate seeds for itself. I have a bit uneasy feeling about this... To give a trivial example of what might go wrong, consider a PRNG where output = internal state and the seed simply sets the internal state. Jun 8 '11 at 15:10
• @Jukka: My approach is similar to Antonio's, yet mine is more general. The PRNG in your example (where output = internal state) does not seem to be cryptographically secure. A PRNG is cryptographically secure if its output is computationally indistinguishable from the uniform distribution. See this for more info. PS: The Blum-Blum-Shub PRNG satisfies this condition. Jun 8 '11 at 19:05

You could use a pseudorandom function $f$ such as AES or ChaCha with a single random key, encrypting a counter. Assign each of the $M = 1000$ parallel processes a unique starting value in $\{ 0, 1, \ldots, M - 1 \}$, and then compute the $j$th random block of bits for process $i$ as $f(i + jM)$, i.e. increment the counter in each process by $M$ for every subsequent block of random bits.

This will give you a cryptographic RNG on every process, but it does not necessarily come with a performance cost. AES is fast if you have hardware that supports it, and ChaCha is fast regardless. Of course, you'll want to measure this in your specific setting to be sure.

Both desired properties 1 and 2 are directly satisfied by this. It's moreover convenient that the behavior of the entire system of parallel tasks is controlled by a single "seed" (the key for $f$).

• If I do not care about cryptographic strength, how does ChaCha(counter) compare with, e.g., Mersenne Twister? Is it faster or slower? Does it have have at least as good statistical properties? I tried to google, but failed to find any articles that compare these two in a non-cryptographic context. Jun 13 '11 at 21:30

There is now a jump function for SFMT (a fast Mersenne Twister implementation).

This allows me to initialise 1000 MTs so that there is no cycle overlap. And SFMT should be faster than MTGP. Almost perfect for my purposes.

You can just use 1000 instances of the Mersenne Twister initialized with different seeds.

You can sample the seeds from another Mersenne Twister, or, to be surer of their independence, from the OS cryptographic pseudorandom number generator (/dev/urandom in Linux).

The Mersenne Twister always operates on the same cyclic sequence, the seed controls where you start generating it. With indepenently sampled seeds, each generator will start at different, typically very far points, with a very small probability of intersection.

• So MT has some nice special properties that guarantee that seeding MT with another MT makes sense? Jun 7 '11 at 15:10
• does MT have any provable pseudorandomness properties? Jun 7 '11 at 15:41
• @Jukka: not any I'm aware of. That's why I suggested to use another type of PRNG for seeding if you are particularly afraid of some strange unknown kind of correlations. Jun 8 '11 at 0:45
• @Sasho: the Wikipedia page mentions k-distribution and the large period. Jun 8 '11 at 0:47
• these indirect measures perplex me; is it ever the case that all you want from a PRNG is a large period and $k$-distribution? i doubt that; those are just heuristic sanity checks; contrast with $k$-wise independence which actually is a pseudorandom property that guarantees accuracy in many settings. also even if you combine two PRNG's, you at least should still show that at least the heuristic "randomness" properties hold Jun 8 '11 at 1:15