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.


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.
Generating the next number in any stream is (almost) as fast as generating the next number with $X$.


Put otherwise: can we get multiple independent streams "for free"?

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 $&lt; 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.
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.