1*39010Sbostic /* $Header: rnd.c,v 1.7 88/07/28 19:10:55 arnold Exp $ */
230268Sbostic
3*39010Sbostic # include "random.h"
430268Sbostic
530268Sbostic /*
630268Sbostic * code for when the good (berkeley) random number generator is around
730268Sbostic */
830268Sbostic
9*39010Sbostic long
rnd(num)1030268Sbostic rnd(num)
11*39010Sbostic long num;
1230268Sbostic {
13*39010Sbostic extern long random();
14*39010Sbostic
1530268Sbostic return (random() % num);
1630268Sbostic }
1730268Sbostic
18*39010Sbostic void
srnd(num)1930268Sbostic srnd(num)
20*39010Sbostic long num;
2130268Sbostic {
2230268Sbostic srandom(num);
2330268Sbostic }
2430268Sbostic
2530268Sbostic #ifdef NO_RANDOM
2630268Sbostic
2730268Sbostic #ifndef lint
2830268Sbostic static char sccsid[] = "@(#)random.c 4.2 (Berkeley) 83/01/02";
2930268Sbostic #endif
3030268Sbostic
3130268Sbostic #include <stdio.h>
3230268Sbostic
3330268Sbostic /*
3430268Sbostic * random.c:
3530268Sbostic * An improved random number generation package. In addition to the standard
3630268Sbostic * rand()/srand() like interface, this package also has a special state info
3730268Sbostic * interface. The initstate() routine is called with a seed, an array of
3830268Sbostic * bytes, and a count of how many bytes are being passed in; this array is then
3930268Sbostic * initialized to contain information for random number generation with that
4030268Sbostic * much state information. Good sizes for the amount of state information are
4130268Sbostic * 32, 64, 128, and 256 bytes. The state can be switched by calling the
4230268Sbostic * setstate() routine with the same array as was initiallized with initstate().
4330268Sbostic * By default, the package runs with 128 bytes of state information and
4430268Sbostic * generates far better random numbers than a linear congruential generator.
4530268Sbostic * If the amount of state information is less than 32 bytes, a simple linear
4630268Sbostic * congruential R.N.G. is used.
4730268Sbostic * Internally, the state information is treated as an array of longs; the
4830268Sbostic * zeroeth element of the array is the type of R.N.G. being used (small
4930268Sbostic * integer); the remainder of the array is the state information for the
5030268Sbostic * R.N.G. Thus, 32 bytes of state information will give 7 longs worth of
5130268Sbostic * state information, which will allow a degree seven polynomial. (Note: the
5230268Sbostic * zeroeth word of state information also has some other information stored
5330268Sbostic * in it -- see setstate() for details).
5430268Sbostic * The random number generation technique is a linear feedback shift register
5530268Sbostic * approach, employing trinomials (since there are fewer terms to sum up that
5630268Sbostic * way). In this approach, the least significant bit of all the numbers in
5730268Sbostic * the state table will act as a linear feedback shift register, and will have
5830268Sbostic * period 2^deg - 1 (where deg is the degree of the polynomial being used,
5930268Sbostic * assuming that the polynomial is irreducible and primitive). The higher
6030268Sbostic * order bits will have longer periods, since their values are also influenced
6130268Sbostic * by pseudo-random carries out of the lower bits. The total period of the
6230268Sbostic * generator is approximately deg*(2**deg - 1); thus doubling the amount of
6330268Sbostic * state information has a vast influence on the period of the generator.
6430268Sbostic * Note: the deg*(2**deg - 1) is an approximation only good for large deg,
6530268Sbostic * when the period of the shift register is the dominant factor. With deg
6630268Sbostic * equal to seven, the period is actually much longer than the 7*(2**7 - 1)
6730268Sbostic * predicted by this formula.
6830268Sbostic */
6930268Sbostic
7030268Sbostic
7130268Sbostic
7230268Sbostic /*
7330268Sbostic * For each of the currently supported random number generators, we have a
7430268Sbostic * break value on the amount of state information (you need at least this
7530268Sbostic * many bytes of state info to support this random number generator), a degree
7630268Sbostic * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
7730268Sbostic * the separation between the two lower order coefficients of the trinomial.
7830268Sbostic */
7930268Sbostic
8030268Sbostic #define TYPE_0 0 /* linear congruential */
8130268Sbostic #define BREAK_0 8
8230268Sbostic #define DEG_0 0
8330268Sbostic #define SEP_0 0
8430268Sbostic
8530268Sbostic #define TYPE_1 1 /* x**7 + x**3 + 1 */
8630268Sbostic #define BREAK_1 32
8730268Sbostic #define DEG_1 7
8830268Sbostic #define SEP_1 3
8930268Sbostic
9030268Sbostic #define TYPE_2 2 /* x**15 + x + 1 */
9130268Sbostic #define BREAK_2 64
9230268Sbostic #define DEG_2 15
9330268Sbostic #define SEP_2 1
9430268Sbostic
9530268Sbostic #define TYPE_3 3 /* x**31 + x**3 + 1 */
9630268Sbostic #define BREAK_3 128
9730268Sbostic #define DEG_3 31
9830268Sbostic #define SEP_3 3
9930268Sbostic
10030268Sbostic #define TYPE_4 4 /* x**63 + x + 1 */
10130268Sbostic #define BREAK_4 256
10230268Sbostic #define DEG_4 63
10330268Sbostic #define SEP_4 1
10430268Sbostic
10530268Sbostic
10630268Sbostic /*
10730268Sbostic * Array versions of the above information to make code run faster -- relies
10830268Sbostic * on fact that TYPE_i == i.
10930268Sbostic */
11030268Sbostic
11130268Sbostic #define MAX_TYPES 5 /* max number of types above */
11230268Sbostic
11330268Sbostic static int degrees[ MAX_TYPES ] = { DEG_0, DEG_1, DEG_2,
11430268Sbostic DEG_3, DEG_4 };
11530268Sbostic
11630268Sbostic static int seps[ MAX_TYPES ] = { SEP_0, SEP_1, SEP_2,
11730268Sbostic SEP_3, SEP_4 };
11830268Sbostic
11930268Sbostic
12030268Sbostic
12130268Sbostic /*
12230268Sbostic * Initially, everything is set up as if from :
12330268Sbostic * initstate( 1, &randtbl, 128 );
12430268Sbostic * Note that this initialization takes advantage of the fact that srandom()
12530268Sbostic * advances the front and rear pointers 10*rand_deg times, and hence the
12630268Sbostic * rear pointer which starts at 0 will also end up at zero; thus the zeroeth
12730268Sbostic * element of the state information, which contains info about the current
12830268Sbostic * position of the rear pointer is just
12930268Sbostic * MAX_TYPES*(rptr - state) + TYPE_3 == TYPE_3.
13030268Sbostic */
13130268Sbostic
13230268Sbostic static long randtbl[ DEG_3 + 1 ] = { TYPE_3,
13330268Sbostic 0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342,
13430268Sbostic 0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb,
13530268Sbostic 0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
13630268Sbostic 0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86,
13730268Sbostic 0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7,
13830268Sbostic 0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
13930268Sbostic 0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b,
14030268Sbostic 0xf5ad9d0e, 0x8999220b, 0x27fb47b9 };
14130268Sbostic
14230268Sbostic /*
14330268Sbostic * fptr and rptr are two pointers into the state info, a front and a rear
14430268Sbostic * pointer. These two pointers are always rand_sep places aparts, as they cycle
14530268Sbostic * cyclically through the state information. (Yes, this does mean we could get
14630268Sbostic * away with just one pointer, but the code for random() is more efficient this
14730268Sbostic * way). The pointers are left positioned as they would be from the call
14830268Sbostic * initstate( 1, randtbl, 128 )
14930268Sbostic * (The position of the rear pointer, rptr, is really 0 (as explained above
15030268Sbostic * in the initialization of randtbl) because the state table pointer is set
15130268Sbostic * to point to randtbl[1] (as explained below).
15230268Sbostic */
15330268Sbostic
15430268Sbostic static long *fptr = &randtbl[ SEP_3 + 1 ];
15530268Sbostic static long *rptr = &randtbl[ 1 ];
15630268Sbostic
15730268Sbostic
15830268Sbostic
15930268Sbostic /*
16030268Sbostic * The following things are the pointer to the state information table,
16130268Sbostic * the type of the current generator, the degree of the current polynomial
16230268Sbostic * being used, and the separation between the two pointers.
16330268Sbostic * Note that for efficiency of random(), we remember the first location of
16430268Sbostic * the state information, not the zeroeth. Hence it is valid to access
16530268Sbostic * state[-1], which is used to store the type of the R.N.G.
16630268Sbostic * Also, we remember the last location, since this is more efficient than
16730268Sbostic * indexing every time to find the address of the last element to see if
16830268Sbostic * the front and rear pointers have wrapped.
16930268Sbostic */
17030268Sbostic
17130268Sbostic static long *state = &randtbl[ 1 ];
17230268Sbostic
17330268Sbostic static int rand_type = TYPE_3;
17430268Sbostic static int rand_deg = DEG_3;
17530268Sbostic static int rand_sep = SEP_3;
17630268Sbostic
17730268Sbostic static long *end_ptr = &randtbl[ DEG_3 + 1 ];
17830268Sbostic
17930268Sbostic
18030268Sbostic
18130268Sbostic /*
18230268Sbostic * srandom:
18330268Sbostic * Initialize the random number generator based on the given seed. If the
18430268Sbostic * type is the trivial no-state-information type, just remember the seed.
18530268Sbostic * Otherwise, initializes state[] based on the given "seed" via a linear
18630268Sbostic * congruential generator. Then, the pointers are set to known locations
18730268Sbostic * that are exactly rand_sep places apart. Lastly, it cycles the state
18830268Sbostic * information a given number of times to get rid of any initial dependencies
18930268Sbostic * introduced by the L.C.R.N.G.
19030268Sbostic * Note that the initialization of randtbl[] for default usage relies on
19130268Sbostic * values produced by this routine.
19230268Sbostic */
19330268Sbostic
srandom(x)19430268Sbostic srandom( x )
19530268Sbostic
19630268Sbostic unsigned x;
19730268Sbostic {
19830268Sbostic register int i, j;
19930268Sbostic
20030268Sbostic if( rand_type == TYPE_0 ) {
20130268Sbostic state[ 0 ] = x;
20230268Sbostic }
20330268Sbostic else {
20430268Sbostic j = 1;
20530268Sbostic state[ 0 ] = x;
20630268Sbostic for( i = 1; i < rand_deg; i++ ) {
20730268Sbostic state[i] = 1103515245*state[i - 1] + 12345;
20830268Sbostic }
20930268Sbostic fptr = &state[ rand_sep ];
21030268Sbostic rptr = &state[ 0 ];
21130268Sbostic for( i = 0; i < 10*rand_deg; i++ ) random();
21230268Sbostic }
21330268Sbostic }
21430268Sbostic
21530268Sbostic
21630268Sbostic
21730268Sbostic /*
21830268Sbostic * initstate:
21930268Sbostic * Initialize the state information in the given array of n bytes for
22030268Sbostic * future random number generation. Based on the number of bytes we
22130268Sbostic * are given, and the break values for the different R.N.G.'s, we choose
22230268Sbostic * the best (largest) one we can and set things up for it. srandom() is
22330268Sbostic * then called to initialize the state information.
22430268Sbostic * Note that on return from srandom(), we set state[-1] to be the type
22530268Sbostic * multiplexed with the current value of the rear pointer; this is so
22630268Sbostic * successive calls to initstate() won't lose this information and will
22730268Sbostic * be able to restart with setstate().
22830268Sbostic * Note: the first thing we do is save the current state, if any, just like
22930268Sbostic * setstate() so that it doesn't matter when initstate is called.
23030268Sbostic * Returns a pointer to the old state.
23130268Sbostic */
23230268Sbostic
23330268Sbostic char *
initstate(seed,arg_state,n)23430268Sbostic initstate( seed, arg_state, n )
23530268Sbostic
23630268Sbostic unsigned seed; /* seed for R. N. G. */
23730268Sbostic char *arg_state; /* pointer to state array */
23830268Sbostic int n; /* # bytes of state info */
23930268Sbostic {
24030268Sbostic register char *ostate = (char *)( &state[ -1 ] );
24130268Sbostic
24230268Sbostic if( rand_type == TYPE_0 ) state[ -1 ] = rand_type;
24330268Sbostic else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type;
24430268Sbostic if( n < BREAK_1 ) {
24530268Sbostic if( n < BREAK_0 ) {
24630268Sbostic fprintf( stderr, "initstate: not enough state (%d bytes) with which to do jack; ignored.\n" );
24730268Sbostic return;
24830268Sbostic }
24930268Sbostic rand_type = TYPE_0;
25030268Sbostic rand_deg = DEG_0;
25130268Sbostic rand_sep = SEP_0;
25230268Sbostic }
25330268Sbostic else {
25430268Sbostic if( n < BREAK_2 ) {
25530268Sbostic rand_type = TYPE_1;
25630268Sbostic rand_deg = DEG_1;
25730268Sbostic rand_sep = SEP_1;
25830268Sbostic }
25930268Sbostic else {
26030268Sbostic if( n < BREAK_3 ) {
26130268Sbostic rand_type = TYPE_2;
26230268Sbostic rand_deg = DEG_2;
26330268Sbostic rand_sep = SEP_2;
26430268Sbostic }
26530268Sbostic else {
26630268Sbostic if( n < BREAK_4 ) {
26730268Sbostic rand_type = TYPE_3;
26830268Sbostic rand_deg = DEG_3;
26930268Sbostic rand_sep = SEP_3;
27030268Sbostic }
27130268Sbostic else {
27230268Sbostic rand_type = TYPE_4;
27330268Sbostic rand_deg = DEG_4;
27430268Sbostic rand_sep = SEP_4;
27530268Sbostic }
27630268Sbostic }
27730268Sbostic }
27830268Sbostic }
27930268Sbostic state = &( ( (long *)arg_state )[1] ); /* first location */
28030268Sbostic end_ptr = &state[ rand_deg ]; /* must set end_ptr before srandom */
28130268Sbostic srandom( seed );
28230268Sbostic if( rand_type == TYPE_0 ) state[ -1 ] = rand_type;
28330268Sbostic else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type;
28430268Sbostic return( ostate );
28530268Sbostic }
28630268Sbostic
28730268Sbostic
28830268Sbostic
28930268Sbostic /*
29030268Sbostic * setstate:
29130268Sbostic * Restore the state from the given state array.
29230268Sbostic * Note: it is important that we also remember the locations of the pointers
29330268Sbostic * in the current state information, and restore the locations of the pointers
29430268Sbostic * from the old state information. This is done by multiplexing the pointer
29530268Sbostic * location into the zeroeth word of the state information.
29630268Sbostic * Note that due to the order in which things are done, it is OK to call
29730268Sbostic * setstate() with the same state as the current state.
29830268Sbostic * Returns a pointer to the old state information.
29930268Sbostic */
30030268Sbostic
30130268Sbostic char *
setstate(arg_state)30230268Sbostic setstate( arg_state )
30330268Sbostic
30430268Sbostic char *arg_state;
30530268Sbostic {
30630268Sbostic register long *new_state = (long *)arg_state;
30730268Sbostic register int type = new_state[0]%MAX_TYPES;
30830268Sbostic register int rear = new_state[0]/MAX_TYPES;
30930268Sbostic char *ostate = (char *)( &state[ -1 ] );
31030268Sbostic
31130268Sbostic if( rand_type == TYPE_0 ) state[ -1 ] = rand_type;
31230268Sbostic else state[ -1 ] = MAX_TYPES*(rptr - state) + rand_type;
31330268Sbostic switch( type ) {
31430268Sbostic case TYPE_0:
31530268Sbostic case TYPE_1:
31630268Sbostic case TYPE_2:
31730268Sbostic case TYPE_3:
31830268Sbostic case TYPE_4:
31930268Sbostic rand_type = type;
32030268Sbostic rand_deg = degrees[ type ];
32130268Sbostic rand_sep = seps[ type ];
32230268Sbostic break;
32330268Sbostic
32430268Sbostic default:
32530268Sbostic fprintf( stderr, "setstate: state info has been munged; not changed.\n" );
32630268Sbostic }
32730268Sbostic state = &new_state[ 1 ];
32830268Sbostic if( rand_type != TYPE_0 ) {
32930268Sbostic rptr = &state[ rear ];
33030268Sbostic fptr = &state[ (rear + rand_sep)%rand_deg ];
33130268Sbostic }
33230268Sbostic end_ptr = &state[ rand_deg ]; /* set end_ptr too */
33330268Sbostic return( ostate );
33430268Sbostic }
33530268Sbostic
33630268Sbostic
33730268Sbostic
33830268Sbostic /*
33930268Sbostic * random:
34030268Sbostic * If we are using the trivial TYPE_0 R.N.G., just do the old linear
34130268Sbostic * congruential bit. Otherwise, we do our fancy trinomial stuff, which is the
34230268Sbostic * same in all ther other cases due to all the global variables that have been
34330268Sbostic * set up. The basic operation is to add the number at the rear pointer into
34430268Sbostic * the one at the front pointer. Then both pointers are advanced to the next
34530268Sbostic * location cyclically in the table. The value returned is the sum generated,
34630268Sbostic * reduced to 31 bits by throwing away the "least random" low bit.
34730268Sbostic * Note: the code takes advantage of the fact that both the front and
34830268Sbostic * rear pointers can't wrap on the same call by not testing the rear
34930268Sbostic * pointer if the front one has wrapped.
35030268Sbostic * Returns a 31-bit random number.
35130268Sbostic */
35230268Sbostic
35330268Sbostic long
random()35430268Sbostic random()
35530268Sbostic {
35630268Sbostic long i;
35730268Sbostic
35830268Sbostic if( rand_type == TYPE_0 ) {
35930268Sbostic i = state[0] = ( state[0]*1103515245 + 12345 )&0x7fffffff;
36030268Sbostic }
36130268Sbostic else {
36230268Sbostic *fptr += *rptr;
36330268Sbostic i = (*fptr >> 1)&0x7fffffff; /* chucking least random bit */
36430268Sbostic if( ++fptr >= end_ptr ) {
36530268Sbostic fptr = state;
36630268Sbostic ++rptr;
36730268Sbostic }
36830268Sbostic else {
36930268Sbostic if( ++rptr >= end_ptr ) rptr = state;
37030268Sbostic }
37130268Sbostic }
37230268Sbostic return( i );
37330268Sbostic }
37430268Sbostic
375*39010Sbostic #endif /* NO_RANDOM */
376