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