-
Notifications
You must be signed in to change notification settings - Fork 0
/
rngs.C
executable file
·180 lines (159 loc) · 7.09 KB
/
rngs.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
/* -------------------------------------------------------------------------
* 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 <stdio.h>
#include <time.h>
#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");
}