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 7 #if defined(LIBC_SCCS) && !defined(lint) 8 static char sccsid[] = "@(#)random.c 5.3 (Berkeley) 04/06/87"; 9 #endif LIBC_SCCS and not lint 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 151 static long *state = &randtbl[ 1 ]; 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", 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 355