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