# How can I seed a parallel linear congruential pseudo-random number generator for maximal period?

Normally when I seed a sequential random number generator in C, I use the call

srand(time(NULL))


then use

rand() mod N


to obtain a random integer between 0 and N-1. However, when I do this in parallel, the calls to time(NULL) are so close to each other that they end up being exactly the same number.

I have tried using a linear congruential random number generator:

$x_{n+1} := a x_n + c\;\;\; (\bmod m)$ for some integers $a,c,$ and $m$.

I know that choosing $m=2^k$ for some large integer $k$ produces fast results because the modulus operator can be calculated by truncating digits. However, I find it difficult to establish seeds that produce random sequences with a large period in parallel. I know that a period length is maximal if

1. c and m are primes with respect to eachother
2. a-1 is divisible by all prime factors of m
3. if m is a multiple of 4, a-1 must also be a multiple of 4.

(source: wikipedia)

But how do I ensure that all random number streams have this maximal property? In terms of MPI, how do I incorporate the rank and size to produce maximal periods using the linear congruential method? Would it be easier to use a Lagged Fibonacci or Mersenne Twister to produce longer parallel random streams?

-
Just as a side note you almost certainly do not want to use a linear congruent PRNG for scientific computation. They can't properly sample spaces of dimension higher than 1. That is to say: they can't even draw a proper sample of points in a plane. – dmckee Feb 12 '12 at 17:09
@dmckee: Marsaglia's theorem is certainly relevant to some scientific computation, but it would be unfair to say it disqualifies LCG's for all scientific computation. Sometimes having a fast PRNG is equally important, and the number of spanning vectors (through the origin) is often more important than the number of hyperplanes that cover the sample space. – Jack Poulson Feb 12 '12 at 17:24
@dmckee: you may be right about LCG, but for lots of applications, 1D is enough. – Paul Feb 12 '12 at 17:28
@JackPerhaps I'm biased by the nature of my work, and I am aware of that which is why I specified the nature of the defect. Paul: There are a variety of well studied PRNGs, and lots of trade off (speed, period, cryptographic safety, and ability to draw high dimensionality tuples). I use Mersenne Twister both in stuff I code by hand and as the default PRNG in ROOT. It's not the fastest one going, but then drawing numbers is a generally a modest contribution to my running time. – dmckee Feb 12 '12 at 17:32
As another side note, if you need to use LC PRNGs (perhaps for speed reasons) you definitely don't want to use mod to grab the low order bits - as Jonathan Dursi suggested, they are much less random. Instead divide your (int) random number by maxint/range to get the range you need. It costs you a divide, but it's probably a cheaper option to improve the quality of your random number stream than switching to another PRNG. – Mark Booth Feb 13 '12 at 11:24

The trick is to interleave each process's LCG stream: for $p$ processes, we modify the LCG

$$x_{n+1} := a x_n + c\;\;\; (\bmod m),$$

to be

$$x_{n+p} := A_p x_n + C_p\;\;\; (\bmod m),$$

where $A_p$ and $C_p$ effectively step forward $p$ steps. We can quickly derive them by expanding the original LCG step:

$$x_{n+2} = a (a x_n + c) + c\;\;\; (\bmod m) = a^2 x_n + (a+1)c\;\;\; (\bmod m),$$

and the pattern is that $A_p \equiv a^p \mod m,$ and $C_p$ is the result of, starting with the number $0$, successively multiplying by $a$ and adding $1$ $p$ times, then multiplying by $c$, all $\bmod m$.

The last step is to ensure that each process's $p$-stride LCG does not overlap: simply initialize the process with rank $r$ with $x_r$ and a parallel LCG with individual periods of $N/p$ is ready, where $N$ is the original period and $p$ is the number of processes. If each process's modified LCG is used equally, the full period is recovered in parallel.

I implemented this about six months ago (perhaps naively), and you can find the code here.

-
That's an interesting approach. It is effectively taking N tuples from a single LC PRNG stream in a distributed way. It still suffers the correlation problems mentioned on wikipedia but requires no synchronisation overhead of a centralised PRNG source. I would be interested to see how the quality of these streams compare to correlation between streams created by multiple LC PRNGs with different constants. – Mark Booth Feb 13 '12 at 11:40
Nice idea; this seems to beg to be written for GPU. – Chinasaur Sep 2 '13 at 7:32
Is your implementation tuned for memory coherence? I imagine trying to give each thread a larger contiguous block and have them jump over each other in bigger steps will work more smoothly? On the GPU, on the other hand, the fully strided version is perfect. – Chinasaur Sep 7 '13 at 3:55

There's a very nice tutorial overview paper by Katzgrabber, Random Numbers In Scientific Computing: An Introduction, which I point people to who want to be a user of PRNGs for scientific computation. Linear congruential generators are fast, but that's about all they have going for them; they have short periods, and they can very easily go wrong; perfectly reasonable looking combinations of a, c, and m can end up with horrifically correlated outputs, even if you satisfy the usual requirements between a, c, and m.

Worse, in one common case where m is a power of two (so the mod operation is fast), the lower order bits have a much shorter period than the sequence as a whole, so if you're doing rand() % N, you have an even shorter period than you would expect.

As a general rule of thumb, lagged-fibbonacci, MT, and WELL generators have much better properties and they're still pretty fast.

In terms of seeding in parallel, Jack Poulson's method is nice because it gives a well-defined sequence of numbers divided evenly between the processors. If that doesn't matter, you can do anything reasonable to seed the different PRNGs; the same paper suggests something lots of people have come up with independantly, hashing the PID or MPI task number with the time. The particular formula suggested there is

long seedgen(void)  {
long s, seed, pid;

pid = getpid();
s = time ( &seconds );

seed = abs(((s*181)*((pid-83)*359))%104729);
return seed;
}


I have no particular opinions about that specific implementation, but the general approach is certainly reasonable.

-
I like the overview paper, thanks! I suppose I should defend LCG's by arguing that one can now simply pick good values for $(a,c,m)$ from Wikipedia or Knuth's Seminumerical Algorithms, and that the example in the Katzgrabber paper is a little unfair since it is an LCG without a shift ($c=0$). – Jack Poulson Feb 13 '12 at 3:14
There's that one, but even the one with a=106, c=1283, m=6075 (fig 2) is a whole pile of fail. But yes, there are known good triples available. – user389 Feb 13 '12 at 16:37
@JackPoulson: I have always found it to be very frustrating having to pick a,c, and m off the top of my head... when I do it this way, it always seems to lead to very small periods. Thanks for the citation! I'll look up Knuth's Seminumerical Algorithms. – Paul Mar 4 '12 at 15:12

A simple idea for spreading a typical sequential RNG over a decent number of threads is to have a single thread advance the seed as fast as possible and send only every thousandth or so seed out to memory. Then have each of your other threads pick up one of these spaced out reference seeds and process the 1000 values in that block, i.e. regenerate again the 1000 seeds in the block, generate their psuedorandom draws, and then do whatever other processing your task involves.

This works because for RNGs that don't compute all that much (LCG is certainly one, but many others should be in this category) the real bottleneck is sending the seeds out to memory (and also probably subsequent processing). If you run a LCG without sending anything out to memory, the whole thing ought to stay in CPU registers and be extremely fast. Even for a more complicated RNG, you should stay in L1 cache and be very fast.

I've used this very simple approach with an LCG that for legacy reasons we have to keep. We basically get linear speedup up to the 4-8 threads in a typical multicore workstation. But now I will try the method from Jack Poulson's answer and hopefully get even faster :).

OTOH, I believe this simple trick should work for other inherently sequential RNGs.

-

This answer is late in coming, but you should have a look at SPRNG. It's specifically designed for scalability in parallel and supports a handful of types of PRNGs.

-