Program: FIBMULT, a pseudo-random number generator using a multiplicative lagged Fibonacci algorithm, implemented in C for sequential and distributed memory parallel computers. Generates uniformly distributed pseudo-random floating point numbers in the interval [0.0, 1.0).
Version: 1.0, February 1995.
Author:
Paul Coddington,
Northeast Parallel Architectures Center
at Syracuse University,
paulc@npac.syr.edu/A>
Code: C code for the generator, plus
test programs which use the generator, and a sample input and
output file, are available at
ftp://nhse.npac.syr.edu/pub/random/generators/fibmult/C/
The Algorithm
The program implements a
lagged Fibonnacci pseudo-random number generator
using multiplication, with lags p and q,
i.e. F(p,q,*) in the notation of
[Marsaglia].
A sequence of odd random integers X(n) is obtained from the formula
X(n) = X(n-p) * X(n-q) mod Mwhere M is a large integer (taken to be 231 here).
Unlike lagged Fibonnacci generators which use subtraction, addition, or XOR, this generator seems to give "good" random numbers for any lag, not necessarily large (see [Coddington]). However larger lags give larger periods and empirically seem to give "better" random numbers (see [Marsaglia], [Knuth], [Coddington]), and thus are recommended.
The period of this generator is (2p-1)*2b-3, where p is the largest lag and b is the number of bits in the integers X(n), i.e. b = log2(M). The lag needs to be chosen large enough to give a period much longer than the number of random numbers to be generated.
Only certain lags (p,q) will give the maximal period of the generator (see [Marsaglia], [Knuth]). Below are some of those lags, with the corresponding periods for 32-bit integers. As a comparison, note that linear congruential generators have a period of about 10^{9} (which is too short for most purposes), or 10^{18} for 64-bit integers (which is fine for most purposes). A Teraflop-year is about 10^{21} floating point operations.
p q period (32 bit) 9689 4187 10^{2925} 4423 2098 10^{1340} 2281 1029 10^{695} 1279 418 10^{393} 607 273 10^{191} 521 168 10^{166} 250 103 10^{84} 127 63 10^{46} 97 33 10^{38} 55 24 10^{25} 43 22 10^{21} 31 13 10^{18} 24 10 10^{16} 17 5 10^{14} 7 3 10^{11} 5 2 10^{10}I recommend using at least (55,24), and ideally much larger lags such as (2281,1029). Larger lags should give higher quality random numbers, with the only possible drawbacks being increased memory use and initialization time. However the memory use is only of order p words, so p=2281 will add less than 10 Kbytes to the memory required for the user program.
Note that very little is known theoretically about this generator, however it performs very well in all empirical statistical tests (see [Marsaglia], [Coddington]). This does not necessarily mean it will perform well for your application. It is a good idea to check your results by also running your program using a good linear congruential generator with at least 48-bit arithmetic, such as RANF or DRAND48 which are standard on Unix systems, or a combined linear congruential generator such as the Lecuyer generator (see [James], [Coddington]).
Note that since different generators are being run on different processors, the results will be dependent on the number of processors. Thus a program that is run on a certain number of processors and checkpointed (i.e. the current state of the random number generator is stored) cannot be restarted using a different number of processors.
The program accesses the "lag table", i.e. the array storing the previous p values in the sequence, using a "circular list", i.e. the elements of the array are not moved with every new addition, only the pointers to the (n-p)th and (n-q)th value are changed.
In order for the sequence to be random, the initial values of the lag table have to be "random enough". There have not to my knowledge been any quantitive tests of this point, however initializing the lag table using a good linear congruential generator (here I use the Unix/C generator RAND) seems empirically to work fine.
The computation required to generate a new random number, apart from updating 2 pointers, is 1 integer multiplication, which multiplies two 32-bit integers and returns the answer modulo 232 (this is the standard for unsigned ints in C), and 1 floating point multiplication, to convert the integers into a floating point number in [0,1).
The routines in the program are called using a standard interface. These include routines for reading and writing the lag table and lag pointers, to allow checkpointing of the user program.
Note that you MUST call the initialization subroutine before calling the generator.
F. James, "A Review of Pseudo-random Number Generators",
Comput. Phys. Comm. 60 (1990) 329.
D.E. Knuth, "The Art of Computer Programming Vol. 2: Seminumerical
Methods", (Addison-Wesley, Reading, Mass., 1981).
P. L'Ecuyer, "Random numbers for simulation", Comm. ACM 33:10 (1990) 85.