/* ------------------------------------------------------------------------- * This is an ANSI C library for multi-stream random number generation. * The use of this library is recommended as a replacement for the ANSI C * rand() and srand() functions, particularly in simulation applications * where the statistical 'goodness' of the random number generator is * important. The library supplies 256 streams of random numbers; use * SelectStream(s) to switch between streams indexed s = 0,1,...,255. * * The streams must be initialized. The recommended way to do this is by * using the function PlantSeeds(x) with the value of x used to initialize * the default stream and all other streams initialized automatically with * values dependent on the value of x. The following convention is used * to initialize the default stream: * if x > 0 then x is the state * if x < 0 then the state is obtained from the system clock * if x = 0 then the state is to be supplied interactively. * * The generator used in this library is a so-called 'Lehmer random number * generator' which returns a pseudo-random number uniformly distributed * 0.0 and 1.0. The period is (m - 1) where m = 2,147,483,647 and the * smallest and largest possible values are (1 / m) and 1 - (1 / m) * respectively. For more details see: * * "Random Number Generators: Good Ones Are Hard To Find" * Steve Park and Keith Miller * Communications of the ACM, October 1988 * * Name : rngs.c (Random Number Generation - Multiple Streams) * Authors : Steve Park & Dave Geyer * Language : ANSI C * Latest Revision : 09-22-98 * ------------------------------------------------------------------------- */ #include #include #include "rngs.h" #define MODULUS 2147483647 /* DON'T CHANGE THIS VALUE */ #define MULTIPLIER 48271 /* DON'T CHANGE THIS VALUE */ #define CHECK 399268537 /* DON'T CHANGE THIS VALUE */ #define STREAMS 256 /* # of streams, DON'T CHANGE THIS VALUE */ #define A256 22925 /* jump multiplier, DON'T CHANGE THIS VALUE */ #define DEFAULT 123456789 /* initial seed, use 0 < DEFAULT < MODULUS */ static long seed[STREAMS] = {DEFAULT}; /* current state of each stream */ static int stream = 0; /* stream index, 0 is the default */ static int initialized = 0; /* test for stream initialization */ double Random(void) /* ---------------------------------------------------------------- * Random returns a pseudo-random real number uniformly distributed * between 0.0 and 1.0. * ---------------------------------------------------------------- */ { const long Q = MODULUS / MULTIPLIER; const long R = MODULUS % MULTIPLIER; long t; t = MULTIPLIER * (seed[stream] % Q) - R * (seed[stream] / Q); if (t > 0) seed[stream] = t; else seed[stream] = t + MODULUS; return ((double) seed[stream] / MODULUS); } void PlantSeeds(long x) /* --------------------------------------------------------------------- * Use this function to set the state of all the random number generator * streams by "planting" a sequence of states (seeds), one per stream, * with all states dictated by the state of the default stream. * The sequence of planted states is separated one from the next by * 8,367,782 calls to Random(). * --------------------------------------------------------------------- */ { const long Q = MODULUS / A256; const long R = MODULUS % A256; int j; int s; initialized = 1; s = stream; /* remember the current stream */ SelectStream(0); /* change to stream 0 */ PutSeed(x); /* set seed[0] */ stream = s; /* reset the current stream */ for (j = 1; j < STREAMS; j++) { x = A256 * (seed[j - 1] % Q) - R * (seed[j - 1] / Q); if (x > 0) seed[j] = x; else seed[j] = x + MODULUS; } } void PutSeed(long x) /* --------------------------------------------------------------- * Use this function to set the state of the current random number * generator stream according to the following conventions: * if x > 0 then x is the state (unless too large) * if x < 0 then the state is obtained from the system clock * if x = 0 then the state is to be supplied interactively * --------------------------------------------------------------- */ { char ok = 0; if (x > 0) x = x % MODULUS; /* correct if x is too large */ if (x < 0) x = ((unsigned long) time((time_t *) NULL)) % MODULUS; if (x == 0) while (!ok) { printf("\nEnter a positive integer seed (9 digits or less) >> "); scanf("%ld", &x); ok = (0 < x) && (x < MODULUS); if (!ok) printf("\nInput out of range ... try again\n"); } seed[stream] = x; } void GetSeed(long *x) /* --------------------------------------------------------------- * Use this function to get the state of the current random number * generator stream. * --------------------------------------------------------------- */ { *x = seed[stream]; } void SelectStream(int index) /* ------------------------------------------------------------------ * Use this function to set the current random number generator * stream -- that stream from which the next random number will come. * ------------------------------------------------------------------ */ { stream = ((unsigned int) index) % STREAMS; if ((initialized == 0) && (stream != 0)) /* protect against */ PlantSeeds(DEFAULT); /* un-initialized streams */ } void TestRandom(void) /* ------------------------------------------------------------------ * Use this (optional) function to test for a correct implementation. * ------------------------------------------------------------------ */ { long i; long x; double u; char ok = 0; SelectStream(0); /* select the default stream */ PutSeed(1); /* and set the state to 1 */ for(i = 0; i < 10000; i++) u = Random(); GetSeed(&x); /* get the new state value */ ok = (x == CHECK); /* and check for correctness */ SelectStream(1); /* select stream 1 */ PlantSeeds(1); /* set the state of all streams */ GetSeed(&x); /* get the state of stream 1 */ ok = ok && (x == A256); /* x should be the jump multiplier */ if (ok) printf("\n The implementation of rngs.c is correct.\n\n"); else printf("\n\a ERROR -- the implementation of rngs.c is not correct.\n\n"); }