1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id$"; |
3 |
#endif |
4 |
#ifndef lint |
5 |
static char sccsid[] = "@(#)random.c 1.1 87/12/17 SMI"; /* from UCB 4.2 83/01/02 */ |
6 |
#endif |
7 |
|
8 |
#include <stdio.h> |
9 |
|
10 |
/* |
11 |
* random.c: |
12 |
* An improved random number generation package. In addition to the standard |
13 |
* rand()/srand() like interface, this package also has a special state info |
14 |
* interface. The initstate() routine is called with a seed, an array of |
15 |
* bytes, and a count of how many bytes are being passed in; this array is then |
16 |
* initialized to contain information for random number generation with that |
17 |
* much state information. Good sizes for the amount of state information are |
18 |
* 32, 64, 128, and 256 bytes. The state can be switched by calling the |
19 |
* setstate() routine with the same array as was initiallized with initstate(). |
20 |
* By default, the package runs with 128 bytes of state information and |
21 |
* generates far better random numbers than a linear congruential generator. |
22 |
* If the amount of state information is less than 32 bytes, a simple linear |
23 |
* congruential R.N.G. is used. |
24 |
* Internally, the state information is treated as an array of longs; the |
25 |
* zeroeth element of the array is the type of R.N.G. being used (small |
26 |
* integer); the remainder of the array is the state information for the |
27 |
* R.N.G. Thus, 32 bytes of state information will give 7 longs worth of |
28 |
* state information, which will allow a degree seven polynomial. (Note: the |
29 |
* zeroeth word of state information also has some other information stored |
30 |
* in it -- see setstate() for details). |
31 |
* The random number generation technique is a linear feedback shift register |
32 |
* approach, employing trinomials (since there are fewer terms to sum up that |
33 |
* way). In this approach, the least significant bit of all the numbers in |
34 |
* the state table will act as a linear feedback shift register, and will have |
35 |
* period 2^deg - 1 (where deg is the degree of the polynomial being used, |
36 |
* assuming that the polynomial is irreducible and primitive). The higher |
37 |
* order bits will have longer periods, since their values are also influenced |
38 |
* by pseudo-random carries out of the lower bits. The total period of the |
39 |
* generator is approximately deg*(2**deg - 1); thus doubling the amount of |
40 |
* state information has a vast influence on the period of the generator. |
41 |
* Note: the deg*(2**deg - 1) is an approximation only good for large deg, |
42 |
* when the period of the shift register is the dominant factor. With deg |
43 |
* equal to seven, the period is actually much longer than the 7*(2**7 - 1) |
44 |
* predicted by this formula. |
45 |
*/ |
46 |
|
47 |
|
48 |
|
49 |
/* |
50 |
* For each of the currently supported random number generators, we have a |
51 |
* break value on the amount of state information (you need at least this |
52 |
* many bytes of state info to support this random number generator), a degree |
53 |
* for the polynomial (actually a trinomial) that the R.N.G. is based on, and |
54 |
* the separation between the two lower order coefficients of the trinomial. |
55 |
*/ |
56 |
|
57 |
#define TYPE_0 0 /* linear congruential */ |
58 |
#define BREAK_0 8 |
59 |
#define DEG_0 0 |
60 |
#define SEP_0 0 |
61 |
|
62 |
#define TYPE_1 1 /* x**7 + x**3 + 1 */ |
63 |
#define BREAK_1 32 |
64 |
#define DEG_1 7 |
65 |
#define SEP_1 3 |
66 |
|
67 |
#define TYPE_2 2 /* x**15 + x + 1 */ |
68 |
#define BREAK_2 64 |
69 |
#define DEG_2 15 |
70 |
#define SEP_2 1 |
71 |
|
72 |
#define TYPE_3 3 /* x**31 + x**3 + 1 */ |
73 |
#define BREAK_3 128 |
74 |
#define DEG_3 31 |
75 |
#define SEP_3 3 |
76 |
|
77 |
#define TYPE_4 4 /* x**63 + x + 1 */ |
78 |
#define BREAK_4 256 |
79 |
#define DEG_4 63 |
80 |
#define SEP_4 1 |
81 |
|
82 |
extern long random(); |
83 |
|
84 |
/* |
85 |
* Array versions of the above information to make code run faster -- relies |
86 |
* on fact that TYPE_i == i. |
87 |
*/ |
88 |
|
89 |
#define MAX_TYPES 5 /* max number of types above */ |
90 |
|
91 |
static int degrees[ MAX_TYPES ] = { DEG_0, DEG_1, DEG_2, |
92 |
DEG_3, DEG_4 }; |
93 |
|
94 |
static int seps[ MAX_TYPES ] = { SEP_0, SEP_1, SEP_2, |
95 |
SEP_3, SEP_4 }; |
96 |
|
97 |
|
98 |
|
99 |
/* |
100 |
* Initially, everything is set up as if from : |
101 |
* initstate( 1, &randtbl, 128 ); |
102 |
* Note that this initialization takes advantage of the fact that srandom() |
103 |
* advances the front and rear pointers 10*rand_deg times, and hence the |
104 |
* rear pointer which starts at 0 will also end up at zero; thus the zeroeth |
105 |
* element of the state information, which contains info about the current |
106 |
* position of the rear pointer is just |
107 |
* MAX_TYPES*(rptr - state) + TYPE_3 == TYPE_3. |
108 |
*/ |
109 |
|
110 |
static long randtbl[ DEG_3 + 1 ] = { TYPE_3, |
111 |
0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342, |
112 |
0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb, |
113 |
0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd, |
114 |
0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86, |
115 |
0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7, |
116 |
0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc, |
117 |
0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b, |
118 |
0xf5ad9d0e, 0x8999220b, 0x27fb47b9 }; |
119 |
|
120 |
/* |
121 |
* fptr and rptr are two pointers into the state info, a front and a rear |
122 |
* pointer. These two pointers are always rand_sep places aparts, as they cycle |
123 |
* cyclically through the state information. (Yes, this does mean we could get |
124 |
* away with just one pointer, but the code for random() is more efficient this |
125 |
* way). The pointers are left positioned as they would be from the call |
126 |
* initstate( 1, randtbl, 128 ) |
127 |
* (The position of the rear pointer, rptr, is really 0 (as explained above |
128 |
* in the initialization of randtbl) because the state table pointer is set |
129 |
* to point to randtbl[1] (as explained below). |
130 |
*/ |
131 |
|
132 |
static long *fptr = &randtbl[ SEP_3 + 1 ]; |
133 |
static long *rptr = &randtbl[ 1 ]; |
134 |
|
135 |
|
136 |
|
137 |
/* |
138 |
* The following things are the pointer to the state information table, |
139 |
* the type of the current generator, the degree of the current polynomial |
140 |
* being used, and the separation between the two pointers. |
141 |
* Note that for efficiency of random(), we remember the first location of |
142 |
* the state information, not the zeroeth. Hence it is valid to access |
143 |
* state[-1], which is used to store the type of the R.N.G. |
144 |
* Also, we remember the last location, since this is more efficient than |
145 |
* indexing every time to find the address of the last element to see if |
146 |
* the front and rear pointers have wrapped. |
147 |
*/ |
148 |
|
149 |
static long *state = &randtbl[ 1 ]; |
150 |
|
151 |
static int rand_type = TYPE_3; |
152 |
static int rand_deg = DEG_3; |
153 |
static int rand_sep = SEP_3; |
154 |
|
155 |
static long *end_ptr = &randtbl[ DEG_3 + 1 ]; |
156 |
|
157 |
|
158 |
|
159 |
/* |
160 |
* srandom: |
161 |
* Initialize the random number generator based on the given seed. If the |
162 |
* type is the trivial no-state-information type, just remember the seed. |
163 |
* Otherwise, initializes state[] based on the given "seed" via a linear |
164 |
* congruential generator. Then, the pointers are set to known locations |
165 |
* that are exactly rand_sep places apart. Lastly, it cycles the state |
166 |
* information a given number of times to get rid of any initial dependencies |
167 |
* introduced by the L.C.R.N.G. |
168 |
* Note that the initialization of randtbl[] for default usage relies on |
169 |
* values produced by this routine. |
170 |
*/ |
171 |
|
172 |
srandom( x ) |
173 |
|
174 |
unsigned x; |
175 |
{ |
176 |
register int i, j; |
177 |
|
178 |
if( rand_type == TYPE_0 ) { |
179 |
state[ 0 ] = x; |
180 |
} |
181 |
else { |
182 |
j = 1; |
183 |
state[ 0 ] = x; |
184 |
for( i = 1; i < rand_deg; i++ ) { |
185 |
state[i] = 1103515245*state[i - 1] + 12345; |
186 |
} |
187 |
fptr = &state[ rand_sep ]; |
188 |
rptr = &state[ 0 ]; |
189 |
for( i = 0; i < 10*rand_deg; i++ ) random(); |
190 |
} |
191 |
} |
192 |
|
193 |
|
194 |
|
195 |
/* |
196 |
* initstate: |
197 |
* Initialize the state information in the given array of n bytes for |
198 |
* future random number generation. Based on the number of bytes we |
199 |
* are given, and the break values for the different R.N.G.'s, we choose |
200 |
* the best (largest) one we can and set things up for it. srandom() is |
201 |
* then called to initialize the state information. |
202 |
* Note that on return from srandom(), we set state[-1] to be the type |
203 |
* multiplexed with the current value of the rear pointer; this is so |
204 |
* successive calls to initstate() won't lose this information and will |
205 |
* be able to restart with setstate(). |
206 |
* Note: the first thing we do is save the current state, if any, just like |
207 |
* setstate() so that it doesn't matter when initstate is called. |
208 |
* Returns a pointer to the old state. |
209 |
*/ |
210 |
|
211 |
char * |
212 |
initstate( seed, arg_state, n ) |
213 |
|
214 |
unsigned seed; /* seed for R. N. G. */ |
215 |
char *arg_state; /* pointer to state array */ |
216 |
int n; /* # bytes of state info */ |
217 |
{ |
218 |
register char *ostate = (char *)( &state[ -1 ] ); |
219 |
|
220 |
if( rand_type == TYPE_0 ) state[ -1 ] = rand_type; |
221 |
else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type; |
222 |
if( n < BREAK_1 ) { |
223 |
if( n < BREAK_0 ) { |
224 |
fprintf( stderr, "initstate: not enough state (%d bytes) with which to do jack; ignored.\n" ); |
225 |
return; |
226 |
} |
227 |
rand_type = TYPE_0; |
228 |
rand_deg = DEG_0; |
229 |
rand_sep = SEP_0; |
230 |
} |
231 |
else { |
232 |
if( n < BREAK_2 ) { |
233 |
rand_type = TYPE_1; |
234 |
rand_deg = DEG_1; |
235 |
rand_sep = SEP_1; |
236 |
} |
237 |
else { |
238 |
if( n < BREAK_3 ) { |
239 |
rand_type = TYPE_2; |
240 |
rand_deg = DEG_2; |
241 |
rand_sep = SEP_2; |
242 |
} |
243 |
else { |
244 |
if( n < BREAK_4 ) { |
245 |
rand_type = TYPE_3; |
246 |
rand_deg = DEG_3; |
247 |
rand_sep = SEP_3; |
248 |
} |
249 |
else { |
250 |
rand_type = TYPE_4; |
251 |
rand_deg = DEG_4; |
252 |
rand_sep = SEP_4; |
253 |
} |
254 |
} |
255 |
} |
256 |
} |
257 |
state = &( ( (long *)arg_state )[1] ); /* first location */ |
258 |
end_ptr = &state[ rand_deg ]; /* must set end_ptr before srandom */ |
259 |
srandom( seed ); |
260 |
if( rand_type == TYPE_0 ) state[ -1 ] = rand_type; |
261 |
else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type; |
262 |
return( ostate ); |
263 |
} |
264 |
|
265 |
|
266 |
|
267 |
/* |
268 |
* setstate: |
269 |
* Restore the state from the given state array. |
270 |
* Note: it is important that we also remember the locations of the pointers |
271 |
* in the current state information, and restore the locations of the pointers |
272 |
* from the old state information. This is done by multiplexing the pointer |
273 |
* location into the zeroeth word of the state information. |
274 |
* Note that due to the order in which things are done, it is OK to call |
275 |
* setstate() with the same state as the current state. |
276 |
* Returns a pointer to the old state information. |
277 |
*/ |
278 |
|
279 |
char * |
280 |
setstate( arg_state ) |
281 |
|
282 |
char *arg_state; |
283 |
{ |
284 |
register long *new_state = (long *)arg_state; |
285 |
register int type = new_state[0]%MAX_TYPES; |
286 |
register int rear = new_state[0]/MAX_TYPES; |
287 |
char *ostate = (char *)( &state[ -1 ] ); |
288 |
|
289 |
if( rand_type == TYPE_0 ) state[ -1 ] = rand_type; |
290 |
else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type; |
291 |
switch( type ) { |
292 |
case TYPE_0: |
293 |
case TYPE_1: |
294 |
case TYPE_2: |
295 |
case TYPE_3: |
296 |
case TYPE_4: |
297 |
rand_type = type; |
298 |
rand_deg = degrees[ type ]; |
299 |
rand_sep = seps[ type ]; |
300 |
break; |
301 |
|
302 |
default: |
303 |
fprintf( stderr, "setstate: state info has been munged; not changed.\n" ); |
304 |
} |
305 |
state = &new_state[ 1 ]; |
306 |
if( rand_type != TYPE_0 ) { |
307 |
rptr = &state[ rear ]; |
308 |
fptr = &state[ (rear + rand_sep)%rand_deg ]; |
309 |
} |
310 |
end_ptr = &state[ rand_deg ]; /* set end_ptr too */ |
311 |
return( ostate ); |
312 |
} |
313 |
|
314 |
|
315 |
|
316 |
/* |
317 |
* random: |
318 |
* If we are using the trivial TYPE_0 R.N.G., just do the old linear |
319 |
* congruential bit. Otherwise, we do our fancy trinomial stuff, which is the |
320 |
* same in all ther other cases due to all the global variables that have been |
321 |
* set up. The basic operation is to add the number at the rear pointer into |
322 |
* the one at the front pointer. Then both pointers are advanced to the next |
323 |
* location cyclically in the table. The value returned is the sum generated, |
324 |
* reduced to 31 bits by throwing away the "least random" low bit. |
325 |
* Note: the code takes advantage of the fact that both the front and |
326 |
* rear pointers can't wrap on the same call by not testing the rear |
327 |
* pointer if the front one has wrapped. |
328 |
* Returns a 31-bit random number. |
329 |
*/ |
330 |
|
331 |
long |
332 |
random() |
333 |
{ |
334 |
long i; |
335 |
|
336 |
if( rand_type == TYPE_0 ) { |
337 |
i = state[0] = ( state[0]*1103515245 + 12345 )&0x7fffffff; |
338 |
} |
339 |
else { |
340 |
*fptr += *rptr; |
341 |
i = (*fptr >> 1)&0x7fffffff; /* chucking least random bit */ |
342 |
if( ++fptr >= end_ptr ) { |
343 |
fptr = state; |
344 |
++rptr; |
345 |
} |
346 |
else { |
347 |
if( ++rptr >= end_ptr ) rptr = state; |
348 |
} |
349 |
} |
350 |
return( i ); |
351 |
} |
352 |
|