Random numbers in C: Some suggestions.
George Marsaglia
This posting ends with 17 lines of
C code that provide eight different
in-line random number generators, six for
random 32-bit integers and two for uniform
reals in (0,1) and (-1,1).
Comments are interspersed with that
code. Various combinations of the six in-line
integer generators may put in C expressions to
provide a wide variety of very fast, long period,
well-tested RNG's. I invite comments, feedback,
verifications and timings.
First, there is narrative giving background
for this posting; you may want to skip it.
Narrative:
Having had experience in producing and
testing for randomness in computers,
I am frequently asked to suggest good
random number generators (RNG's), test
RNG's, or comment on existing RNG's. Many
recent queries have been in two areas:
(1) requests for implementations in C and
(2) comments on generators with immense periods,
particularly the Mersenne Twister.
This posting is mainly for category (1),
for which I suggest a set of C implementations
of RNG's I have developed. C implementations
of my DIEHARD battery of tests will be
discussed elsewhere, and Narasimhan's GUI
version is expected to be released soon.
For (2), I merely repeat what I have said
in response to various queries: the Mersenne
Twister looks good, but it seems to be essentially
a lagged Fibonacci RNG using the exclusive-or
(xor) operation, and experience has shown that
lagged Fibonacci generators using xor provide
unsatisfactory 'randomness' unless the lags are
very long, and even for those with very long lags,
(and even for those using + or - rather than xor),
many people (I among them) are inclined to be
cautious about sequences based on such a simple
operation as: each new integer is the xor, (or sum,
or difference), of two earlier ones. To be sure,
the resulting integers can be "twisted", but not,
I think, as simply or as safely as combining, say
by addition, with members of a sequence from a
(shorter period) generator that has itself passed
extensive tests of randomness.
I also reply that it does not take an immense
program (as, for example, in posted listings
of Twister) to produce a more satisfactory RNG
with an immense period, and give this example,
on which I will expand below: Inclusion of
#define SWB ( t[c+237]=(x=t[c+15])-(y=t[++c]+(x<y)) )
together with suitably initialized seeds in
static unsigned long x,y,t[256]; unsigned char c;
will allow you to put the string SWB in any C
expression and it will provide, in about 100 nanosecs,
a 32-bit random integer with period 2^7578. (Here
and below, ^ means exponent, except in C expressions,
where it means xor (exclusive-or).
Now for the (2) part, in which I suggest a number
of C implementations and invite comment and feedback.
Most of these were previously developed and tested
via Fortran versions. I list eight RNG's, each of
them by means of C's powerful #define device. This
provides fast, compact implementation, allows
insertion of the required random variable directly
into an expression, and, finally, provides a good
selection of RNG's for use individually or in
combination. The latter makes it possible to
further confirm what empirical results suggest:
combining two or more RNG's provides better,
(or no worse) randomness, and for encryption enthusiasts:
combination generators are harder to "crack".
For those wishing to try these eight RNG's:
At the top of your C program, include these
definitions and the static variables that follow.
Everything past this line is either C code or comment.
--------------------------------------------------
#define UL unsigned long
#define znew ((z=36969*(z&65535)+(z>>16))<<16)
#define wnew ((w=18000*(w&65535)+(w>>16))&65535)
#define MWC (znew+wnew)
#define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
#define CONG (jcong=69069*jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
#define LFIB4 (t[c]=t[c]+t[c+58]+t[c+119]+t[++c+178])
#define SWB (t[c+237]=(x=t[c+15])-(y=t[++c]+(x<y)))
#define UNI (KISS*2.328306e-10)
#define VNI ((long) KISS)*4.656613e-10
/* Global static variables: */
static UL z=362436069, w=521288629, jsr=123456789, jcong=380116160;
static UL t[256];
static UL x=0,y=0; static unsigned char c=0;
/* Random seeds must be used to reset z,w,jsr,jcong and
the table t[256] Here is an example procedure, using KISS: */
void settable(UL i1,UL i2,UL i3,UL i4)
{ int i; z=i1;w=i2,jsr=i3; jcong=i4;
for(i=0;i<256;i++) t[i]=KISS; }
/* End of C code; Only comments follow. Stick the above
17 lines in your simulation programs, initialize the table,
and have a variety of promising RNG's at your disposal. */
/* You may want use more complicated names for the
above simple 1-letter variable names: z,w,x,y,t,c,
to avoid clashing with favorites in your code. */
/* Any one of KISS, MWC, LFIB4, SWB, SHR3, or CONG
can be used in an expression to provide a random
32-bit integer, and UNI in an expression will
provide a random uniform in (01), or VNI in (-1,1).
For example, for int i, float v; i=(MWC>>24); will
provide a random byte, while v=4.+3.*UNI; will
provide a uniform v in the interval 4. to 7.
For the super cautious, (KISS+SWB) in an expression
would provide a random 32-bit integer from
a sequence with period > 2^7700, and would only
add some 300 nanoseconds to the computing
time for that expression. */
/* The KISS generator, (Keep It Simple Stupid), is
designed to combine the two multiply-with-carry
generators in MWC with the 3-shift register SHR3
and the congruential generator CONG, using
addition and exclusive-or. Period about 2^123. It
is one of my favorite generators. */
/* The MWC generator concatenates two 16-bit
multiply-with-carry generators, x(n)=36969x(n-1)+carry,
y(n)=18000y(n-1)+carry mod 2^16, has period about
2^60 and seems to pass all tests of randomness. A favorite
stand-alone generator---faster than KISS, which contains it.*/
/* SHR3 is a 3-shift-register generator with
period 2^32-1. It uses
y(n)=y(n-1)(I+L^17)(I+R^13)(I+L^5), with the
y's viewed as binary vectors, L the 32x32
binary matrix that shifts a vector left 1, and
R its transpose. SHR3 seems to pass all except
the binary rank test, since 32 successive
values, as binary vectors, must be linearly
independent, while 32 successive truly random
32-bit integers, viewed as binary vectors, will
be linearly independent only about 29% of the time. */
/* CONG is a congruential generator with the
widely used 69069 as multiplier:
x(n)=69069x(n-1)+1234567. It has period 2^32.
The leading half of its 32 bits seem to pass
all tests, but bits in the last half are too
regular. */
/* LFIB4 is an extension of the class that I have
previously defined as lagged Fibonacci
generators: x(n)=x(n-r) op x(n-s), with the x's
in a finite set over which there is a binary
operation op, such as +,- on integers mod 2^32,
* on odd such integers, exclusive-or (xor) on
binary vectors. Except for those using
multiplication, lagged Fibonacci generators
fail various tests of randomness, unless the
lags are very long. To see if more than two
lags would serve to overcome the problems of 2-
lag generators using +,- or xor, I have
developed the 4-lag generator LFIB4:
x(n)=x(n-256)+x(n-179)+x(n-119)+x(n-55) mod 2^32.
Its period is 2^31*(2^256-1), about 2^287, and
it seems to pass all tests---in particular,
those of the kind for which 2-lag generators
using +,-,xor seem to fail. For even more
confidence in its suitability, LFIB4 can be
combined with KISS, with a resulting period of
about 2^410: just use (KISS+LFIB4) in any C
expression. */
/* SWB is a subtract-with-borrow generator that I
developed to give a simple method for producing
extremely long periods:
x(n)=x(n-222)-x(n-237)-borrow mod 2^32.
The 'borrow' is 0 unless set to 1 if computing
x(n-1) caused overflow in 32-bit integer
arithmetic. This generator has a very long
period, 2^7098(2^480-1), about 2^7578. It seems
to pass all tests of randomness, but,
suspicious of a generator so simple and fast
(62 nanosecs at 300MHz), I would suggest
combining SWB with KISS, MWC, SHR3, or CONG. */
/* Finally, because many simulations call for
uniform random variables in 0<v<1 or -1<v<1, I
use #define statements that permit inclusion of
such variates directly in expressions: using
UNI will provide a uniform random real (float)
in (0,1), while VNI will provide one in (-1,1). */
/* All of these: MWC, SHR3, CONG, KISS, LFIB4,
SWB, UNI and VNI, permit direct insertion of
the desired random quantity into an expression,
avoiding the time and space costs of a function
call. I call these in-line-define functions.
To use them, static variables z,w,jsr and
jcong should be assigned seed values other than
their initial values. If LFIB4 or SWB are
used, the static table t[256] must be
initialized. A sample procedure follows. */
/* A note on timing: It is difficult to provide
exact time costs for inclusion of one of these
in-line-define functions in an expression.
Times may differ widely for different
compilers, as the C operations may be deeply
nested and tricky. I suggest these rough
comparisons, based on averaging ten runs of a
routine that is essentially a long loop:
for(i=1;i<10000000;i++) L=KISS; then with KISS
replaced with SHR3, CONG,... or KISS+SWB, etc.
The times on my home PC, a Pentium 300MHz, in
nanoseconds: LFIB4=64; CONG=90; SWB=100;
SHR3=110; KISS=209; KISS+LFIB4=252; KISS+SWB=310. */
Charles Bond
using
a similar routine in high speed simulations for years. Small departures
from
statistically correct randomness are not a problem for my application,
but
speed is. I believe Knuth briefly discussed the method with guarded
approval -- constrained by the concern that there was no real theory
behind it. Do you know if any theoretical work has been done since
Knuth's
book to justify SWB?
George Marsaglia
Charles Bond wrote:
> Thanks for the post. I want to comment on the SWB routine. I've been
> using
> a similar routine in high speed simulations for years. ...
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
Many years? It hasn't been very many years
since I invented the subtract-with-borrow method,
and developed theory for establishing the periods.
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
> . I believe Knuth briefly discussed the method with guarded
> approval -- constrained by the concern that there was no real theory
> behind it. Do you know if any theoretical work has been done since
> Knuth's book to justify SWB?
> &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
> This business of providing 'theoretical support' for
> RNG's tends to be overdone---perhaps
> because of the influence of Knuth. His marvelous
> expository and poweful mathematical skills have
> justifiably made him the leading authority.
> Many people do not understand that such theory
> is based on the simple fact that congruential
> random numbers "fall mainly in the planes",
> that is, points whose coordinates are succesive
> elements from the sequence lie on a lattice
> with a huge unit-cell volume, (m^(n-1) in n-dimensions),
> compared to the unit-cell volume of 1 for truly random
> integers. So the "lattice test", as I called it,
> applies to congruential generators, although the ideas have
> been extended to the near-lattice-like patterns of
> certain other kinds of generators. But there seem
> to be no such lattice-like patterns for many other
> kinds of generators, and even if there were, it
> is an easy matter to destroy such patterns by
> combining with generators having disparate mathematical
> structures.
>
> The quote from my Keynote Address at the 1984
> Computer Science and Statistics: Symposium on the
> Interface, still applies:
>
> "A random number generator is like sex:
> When it's good, its wonderful;
> And when it's bad, it's still pretty good."
>
> Add to that, in line with my recommendations
> on combination generators;
>
> "And if it's bad, try a twosome or threesome."
>
> George Marsaglia
Mike Oliver
> [...] and experience has shown that
> lagged Fibonacci generators using xor provide
> unsatisfactory 'randomness' unless the lags are
> very long, and even for those with very long lags,
> (and even for those using + or - rather than xor),
Could you give us a pointer to information about
why these RNGs are unsatisfactory and what sort
of test they tend to fail?
--
Disclaimer: I could be wrong -- but I'm not. (Eagles, "Victim of Love")
Finger for PGP public key, or visit http://www.math.ucla.edu/~oliver.
1500 bits, fingerprint AE AE 4F F8 EA EA A6 FB E9 36 5F 9E EA D0 F8 B9
Mike Oliver
> [...] and experience has shown that
> lagged Fibonacci generators using xor provide
> unsatisfactory 'randomness' unless the lags are
> very long, and even for those with very long lags,
> (and even for those using + or - rather than xor),
Could you please give us a pointer to information about
why these RNGs are unsatisfactory and what sort
of test they tend to fail?
Posted and mailed.
Charles Bond
generator was *the same* as your subtract with borrow, only that
it was *similar* (high speed, no multiplications). I gladly acknowledge
your claim on it. But you seem a little skeptical of it yourself, and I
was just curious.
Regards,
C. Bond
George Marsaglia wrote:
> Charles Bond wrote:
>
> > Thanks for the post. I want to comment on the SWB routine. I've been
> > using
> > a similar routine in high speed simulations for years. ...
>
> &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
> Many years? It hasn't been very many years
> since I invented the subtract-with-borrow method,
> and developed theory for establishing the periods.
> &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
>
> > . I believe Knuth briefly discussed the method with guarded
> > approval -- constrained by the concern that there was no real theory
> > behind it. Do you know if any theoretical work has been done since
> > Knuth's book to justify SWB?
>
Eric Backus
> This posting ends with 17 lines of
> C code that provide eight different
> in-line random number generators, six for
> random 32-bit integers and two for uniform
> reals in (0,1) and (-1,1).
> Comments are interspersed with that
> code. Various combinations of the six in-line
> integer generators may put in C expressions to
> provide a wide variety of very fast, long period,
> well-tested RNG's. I invite comments, feedback,
> verifications and timings.
> #define UL unsigned long
> #define znew ((z=36969*(z&65535)+(z>>16))<<16)
> #define wnew ((w=18000*(w&65535)+(w>>16))&65535)
> #define MWC (znew+wnew)
> #define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
> #define CONG (jcong=69069*jcong+1234567)
> #define KISS ((MWC^CONG)+SHR3)
> #define LFIB4 (t[c]=t[c]+t[c+58]+t[c+119]+t[++c+178])
> #define SWB (t[c+237]=(x=t[c+15])-(y=t[++c]+(x<y)))
> #define UNI (KISS*2.328306e-10)
> #define VNI ((long) KISS)*4.656613e-10
> /* Global static variables: */
> static UL z=362436069, w=521288629, jsr=123456789, jcong=380116160;
> static UL t[256];
> static UL x=0,y=0; static unsigned char c=0;
>
> /* Random seeds must be used to reset z,w,jsr,jcong and
> the table t[256] Here is an example procedure, using KISS: */
>
> void settable(UL i1,UL i2,UL i3,UL i4)
> { int i; z=i1;w=i2,jsr=i3; jcong=i4;
> for(i=0;i<256;i++) t[i]=KISS; }
Thank you for providing this extremely useful code. (I'd like to make
use of it, however I see no copyright notice, can I assume you are
making it free for anyone to use?)
I have a small problem with the definition of LFIB4 and SWB. In an
attempt to make these a single line of C code, they both use "++c" in
the same expression as they use "c". A C compiler is free to
rearrange the order in which it calculates the intermediate terms of
these expressions, so the expressions can produce different results
depending on the compiler.
I might propose alternate expressions using the "," operator in
order to remove any ambiguity. With a good compiler, these
expressions probably won't be any slower than your original ones:
#define LFIB4_ALT (t[c]=t[c]+t[c+58]+t[c+119]+t[c+179],t[c++])
#define SWB_ALT (t[c+237]=(x=t[c+15])-(y=t[c+1]+(x<y)),t[c++ +237])
However, these are uglier and harder to understand than your original
expressions, and of course I might have made a mistake in interpreting
where the c++ should go. Any comments?
--
Eric Backus <eric_...@hp.com>
http://labejb.lsid.hp.com/
(425) 335-2495
Bernhard Treutwein
I recently stumbled over the following:
ISAAC: a fast cryptographic random number generator
(see http://ourworld.compuserve.com/homepages/bob_jenkins/isaacafa.htm)
can someone comment on it, is it worth a deeper look ?
--
Bernhard Treutwein Tel. +49-89-5996-642, Fax -615
Institut f. Med. Psychologie Ludwig-Maximilians-Universitaet
bern...@imp.med.uni-muenchen.de http://www.lrz.de/~bernhard
-------------------------------- ---------------------------