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