133965Sjdp/* 233965Sjdp * Copyright (c) 1983 Regents of the University of California. 333965Sjdp * All rights reserved. 433965Sjdp * 560484Sobrien * Redistribution and use in source and binary forms, with or without 660484Sobrien * modification, are permitted provided that the following conditions 760484Sobrien * are met: 860484Sobrien * 1. Redistributions of source code must retain the above copyright 960484Sobrien * notice, this list of conditions and the following disclaimer. 1060484Sobrien * 2. Redistributions in binary form must reproduce the above copyright 1160484Sobrien * notice, this list of conditions and the following disclaimer in the 1260484Sobrien * documentation and/or other materials provided with the distribution. 1360484Sobrien * 3. [rescinded 22 July 1999] 1460484Sobrien * 4. Neither the name of the University nor the names of its contributors 1560484Sobrien * may be used to endorse or promote products derived from this software 1660484Sobrien * without specific prior written permission. 1760484Sobrien * 1860484Sobrien * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 1960484Sobrien * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 2060484Sobrien * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 2160484Sobrien * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 2260484Sobrien * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 2360484Sobrien * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 2460484Sobrien * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 2560484Sobrien * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 2660484Sobrien * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 2760484Sobrien * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 2860484Sobrien * SUCH DAMAGE. 2933965Sjdp */ 3033965Sjdp 3133965Sjdp/* 3233965Sjdp * This is derived from the Berkeley source: 3333965Sjdp * @(#)random.c 5.5 (Berkeley) 7/6/88 3433965Sjdp * It was reworked for the GNU C Library by Roland McGrath. 3533965Sjdp */ 3633965Sjdp 3789857Sobrien/* 3889857Sobrien 3989857Sobrien@deftypefn Supplement {long int} random (void) 4089857Sobrien@deftypefnx Supplement void srandom (unsigned int @var{seed}) 4189857Sobrien@deftypefnx Supplement void* initstate (unsigned int @var{seed}, void *@var{arg_state}, unsigned long @var{n}) 4289857Sobrien@deftypefnx Supplement void* setstate (void *@var{arg_state}) 4389857Sobrien 4489857SobrienRandom number functions. @code{random} returns a random number in the 4589857Sobrienrange 0 to @code{LONG_MAX}. @code{srandom} initializes the random 4689857Sobriennumber generator to some starting point determined by @var{seed} 4789857Sobrien(else, the values returned by @code{random} are always the same for each 4889857Sobrienrun of the program). @code{initstate} and @code{setstate} allow fine-grained 4989857Sobriencontrol over the state of the random number generator. 5089857Sobrien 5189857Sobrien@end deftypefn 5289857Sobrien 5389857Sobrien*/ 5489857Sobrien 5533965Sjdp#include <errno.h> 5633965Sjdp 5733965Sjdp#if 0 5833965Sjdp 5933965Sjdp#include <ansidecl.h> 6033965Sjdp#include <limits.h> 6133965Sjdp#include <stddef.h> 6233965Sjdp#include <stdlib.h> 6333965Sjdp 6433965Sjdp#else 6533965Sjdp 6633965Sjdp#define ULONG_MAX ((unsigned long)(~0L)) /* 0xFFFFFFFF for 32-bits */ 6733965Sjdp#define LONG_MAX ((long)(ULONG_MAX >> 1)) /* 0x7FFFFFFF for 32-bits*/ 6833965Sjdp 6933965Sjdp#ifdef __STDC__ 7033965Sjdp# define PTR void * 7160484Sobrien# ifndef NULL 7260484Sobrien# define NULL (void *) 0 7360484Sobrien# endif 7433965Sjdp#else 7533965Sjdp# define PTR char * 7660484Sobrien# ifndef NULL 7760484Sobrien# define NULL (void *) 0 7860484Sobrien# endif 7933965Sjdp#endif 8033965Sjdp 8133965Sjdp#endif 8233965Sjdp 83218822Sdimlong int random (void); 8433965Sjdp 8533965Sjdp/* An improved random number generation package. In addition to the standard 8633965Sjdp rand()/srand() like interface, this package also has a special state info 8733965Sjdp interface. The initstate() routine is called with a seed, an array of 8833965Sjdp bytes, and a count of how many bytes are being passed in; this array is 8933965Sjdp then initialized to contain information for random number generation with 9033965Sjdp that much state information. Good sizes for the amount of state 9133965Sjdp information are 32, 64, 128, and 256 bytes. The state can be switched by 9233965Sjdp calling the setstate() function with the same array as was initiallized 9333965Sjdp with initstate(). By default, the package runs with 128 bytes of state 9433965Sjdp information and generates far better random numbers than a linear 9533965Sjdp congruential generator. If the amount of state information is less than 9633965Sjdp 32 bytes, a simple linear congruential R.N.G. is used. Internally, the 9733965Sjdp state information is treated as an array of longs; the zeroeth element of 9833965Sjdp the array is the type of R.N.G. being used (small integer); the remainder 9933965Sjdp of the array is the state information for the R.N.G. Thus, 32 bytes of 10033965Sjdp state information will give 7 longs worth of state information, which will 10133965Sjdp allow a degree seven polynomial. (Note: The zeroeth word of state 10233965Sjdp information also has some other information stored in it; see setstate 10333965Sjdp for details). The random number generation technique is a linear feedback 10433965Sjdp shift register approach, employing trinomials (since there are fewer terms 10533965Sjdp to sum up that way). In this approach, the least significant bit of all 10633965Sjdp the numbers in the state table will act as a linear feedback shift register, 10733965Sjdp and will have period 2^deg - 1 (where deg is the degree of the polynomial 10833965Sjdp being used, assuming that the polynomial is irreducible and primitive). 10933965Sjdp The higher order bits will have longer periods, since their values are 11033965Sjdp also influenced by pseudo-random carries out of the lower bits. The 11133965Sjdp total period of the generator is approximately deg*(2**deg - 1); thus 11233965Sjdp doubling the amount of state information has a vast influence on the 11333965Sjdp period of the generator. Note: The deg*(2**deg - 1) is an approximation 11433965Sjdp only good for large deg, when the period of the shift register is the 11533965Sjdp dominant factor. With deg equal to seven, the period is actually much 11633965Sjdp longer than the 7*(2**7 - 1) predicted by this formula. */ 11733965Sjdp 11833965Sjdp 11933965Sjdp 12033965Sjdp/* For each of the currently supported random number generators, we have a 12133965Sjdp break value on the amount of state information (you need at least thi 12233965Sjdp bytes of state info to support this random number generator), a degree for 12333965Sjdp the polynomial (actually a trinomial) that the R.N.G. is based on, and 12433965Sjdp separation between the two lower order coefficients of the trinomial. */ 12533965Sjdp 12633965Sjdp/* Linear congruential. */ 12733965Sjdp#define TYPE_0 0 12833965Sjdp#define BREAK_0 8 12933965Sjdp#define DEG_0 0 13033965Sjdp#define SEP_0 0 13133965Sjdp 13233965Sjdp/* x**7 + x**3 + 1. */ 13333965Sjdp#define TYPE_1 1 13433965Sjdp#define BREAK_1 32 13533965Sjdp#define DEG_1 7 13633965Sjdp#define SEP_1 3 13733965Sjdp 13833965Sjdp/* x**15 + x + 1. */ 13933965Sjdp#define TYPE_2 2 14033965Sjdp#define BREAK_2 64 14133965Sjdp#define DEG_2 15 14233965Sjdp#define SEP_2 1 14333965Sjdp 14433965Sjdp/* x**31 + x**3 + 1. */ 14533965Sjdp#define TYPE_3 3 14633965Sjdp#define BREAK_3 128 14733965Sjdp#define DEG_3 31 14833965Sjdp#define SEP_3 3 14933965Sjdp 15033965Sjdp/* x**63 + x + 1. */ 15133965Sjdp#define TYPE_4 4 15233965Sjdp#define BREAK_4 256 15333965Sjdp#define DEG_4 63 15433965Sjdp#define SEP_4 1 15533965Sjdp 15633965Sjdp 15733965Sjdp/* Array versions of the above information to make code run faster. 15833965Sjdp Relies on fact that TYPE_i == i. */ 15933965Sjdp 16033965Sjdp#define MAX_TYPES 5 /* Max number of types above. */ 16133965Sjdp 16233965Sjdpstatic int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 }; 16333965Sjdpstatic int seps[MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 }; 16433965Sjdp 16533965Sjdp 16633965Sjdp 16733965Sjdp/* Initially, everything is set up as if from: 16833965Sjdp initstate(1, randtbl, 128); 16933965Sjdp Note that this initialization takes advantage of the fact that srandom 17033965Sjdp advances the front and rear pointers 10*rand_deg times, and hence the 17133965Sjdp rear pointer which starts at 0 will also end up at zero; thus the zeroeth 17233965Sjdp element of the state information, which contains info about the current 17333965Sjdp position of the rear pointer is just 17433965Sjdp (MAX_TYPES * (rptr - state)) + TYPE_3 == TYPE_3. */ 17533965Sjdp 17633965Sjdpstatic long int randtbl[DEG_3 + 1] = 17733965Sjdp { TYPE_3, 17833965Sjdp 0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342, 17933965Sjdp 0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb, 18033965Sjdp 0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd, 18133965Sjdp 0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86, 18233965Sjdp 0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7, 18333965Sjdp 0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc, 18433965Sjdp 0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b, 18533965Sjdp 0xf5ad9d0e, 0x8999220b, 0x27fb47b9 18633965Sjdp }; 18733965Sjdp 18833965Sjdp/* FPTR and RPTR are two pointers into the state info, a front and a rear 18933965Sjdp pointer. These two pointers are always rand_sep places aparts, as they 19033965Sjdp cycle through the state information. (Yes, this does mean we could get 19133965Sjdp away with just one pointer, but the code for random is more efficient 19233965Sjdp this way). The pointers are left positioned as they would be from the call: 19333965Sjdp initstate(1, randtbl, 128); 19433965Sjdp (The position of the rear pointer, rptr, is really 0 (as explained above 19533965Sjdp in the initialization of randtbl) because the state table pointer is set 19633965Sjdp to point to randtbl[1] (as explained below).) */ 19733965Sjdp 19833965Sjdpstatic long int *fptr = &randtbl[SEP_3 + 1]; 19933965Sjdpstatic long int *rptr = &randtbl[1]; 20033965Sjdp 20133965Sjdp 20233965Sjdp 20333965Sjdp/* The following things are the pointer to the state information table, 20433965Sjdp the type of the current generator, the degree of the current polynomial 20533965Sjdp being used, and the separation between the two pointers. 20633965Sjdp Note that for efficiency of random, we remember the first location of 20733965Sjdp the state information, not the zeroeth. Hence it is valid to access 20833965Sjdp state[-1], which is used to store the type of the R.N.G. 20933965Sjdp Also, we remember the last location, since this is more efficient than 21033965Sjdp indexing every time to find the address of the last element to see if 21133965Sjdp the front and rear pointers have wrapped. */ 21233965Sjdp 21333965Sjdpstatic long int *state = &randtbl[1]; 21433965Sjdp 21533965Sjdpstatic int rand_type = TYPE_3; 21633965Sjdpstatic int rand_deg = DEG_3; 21733965Sjdpstatic int rand_sep = SEP_3; 21833965Sjdp 21933965Sjdpstatic long int *end_ptr = &randtbl[sizeof(randtbl) / sizeof(randtbl[0])]; 22033965Sjdp 22133965Sjdp/* Initialize the random number generator based on the given seed. If the 22233965Sjdp type is the trivial no-state-information type, just remember the seed. 22333965Sjdp Otherwise, initializes state[] based on the given "seed" via a linear 22433965Sjdp congruential generator. Then, the pointers are set to known locations 22533965Sjdp that are exactly rand_sep places apart. Lastly, it cycles the state 22633965Sjdp information a given number of times to get rid of any initial dependencies 22733965Sjdp introduced by the L.C.R.N.G. Note that the initialization of randtbl[] 22833965Sjdp for default usage relies on values produced by this routine. */ 22933965Sjdpvoid 230218822Sdimsrandom (unsigned int x) 23133965Sjdp{ 23233965Sjdp state[0] = x; 23333965Sjdp if (rand_type != TYPE_0) 23433965Sjdp { 23533965Sjdp register long int i; 23633965Sjdp for (i = 1; i < rand_deg; ++i) 23733965Sjdp state[i] = (1103515145 * state[i - 1]) + 12345; 23833965Sjdp fptr = &state[rand_sep]; 23933965Sjdp rptr = &state[0]; 24033965Sjdp for (i = 0; i < 10 * rand_deg; ++i) 24133965Sjdp random(); 24233965Sjdp } 24333965Sjdp} 24433965Sjdp 24533965Sjdp/* Initialize the state information in the given array of N bytes for 24633965Sjdp future random number generation. Based on the number of bytes we 24733965Sjdp are given, and the break values for the different R.N.G.'s, we choose 24833965Sjdp the best (largest) one we can and set things up for it. srandom is 24933965Sjdp then called to initialize the state information. Note that on return 25033965Sjdp from srandom, we set state[-1] to be the type multiplexed with the current 25133965Sjdp value of the rear pointer; this is so successive calls to initstate won't 25233965Sjdp lose this information and will be able to restart with setstate. 25333965Sjdp Note: The first thing we do is save the current state, if any, just like 25433965Sjdp setstate so that it doesn't matter when initstate is called. 25533965Sjdp Returns a pointer to the old state. */ 25633965SjdpPTR 257218822Sdiminitstate (unsigned int seed, PTR arg_state, unsigned long n) 25833965Sjdp{ 25933965Sjdp PTR ostate = (PTR) &state[-1]; 26033965Sjdp 26133965Sjdp if (rand_type == TYPE_0) 26233965Sjdp state[-1] = rand_type; 26333965Sjdp else 26433965Sjdp state[-1] = (MAX_TYPES * (rptr - state)) + rand_type; 26533965Sjdp if (n < BREAK_1) 26633965Sjdp { 26733965Sjdp if (n < BREAK_0) 26833965Sjdp { 26933965Sjdp errno = EINVAL; 27033965Sjdp return NULL; 27133965Sjdp } 27233965Sjdp rand_type = TYPE_0; 27333965Sjdp rand_deg = DEG_0; 27433965Sjdp rand_sep = SEP_0; 27533965Sjdp } 27633965Sjdp else if (n < BREAK_2) 27733965Sjdp { 27833965Sjdp rand_type = TYPE_1; 27933965Sjdp rand_deg = DEG_1; 28033965Sjdp rand_sep = SEP_1; 28133965Sjdp } 28233965Sjdp else if (n < BREAK_3) 28333965Sjdp { 28433965Sjdp rand_type = TYPE_2; 28533965Sjdp rand_deg = DEG_2; 28633965Sjdp rand_sep = SEP_2; 28733965Sjdp } 28833965Sjdp else if (n < BREAK_4) 28933965Sjdp { 29033965Sjdp rand_type = TYPE_3; 29133965Sjdp rand_deg = DEG_3; 29233965Sjdp rand_sep = SEP_3; 29333965Sjdp } 29433965Sjdp else 29533965Sjdp { 29633965Sjdp rand_type = TYPE_4; 29733965Sjdp rand_deg = DEG_4; 29833965Sjdp rand_sep = SEP_4; 29933965Sjdp } 30033965Sjdp 30133965Sjdp state = &((long int *) arg_state)[1]; /* First location. */ 30233965Sjdp /* Must set END_PTR before srandom. */ 30333965Sjdp end_ptr = &state[rand_deg]; 30433965Sjdp srandom(seed); 30533965Sjdp if (rand_type == TYPE_0) 30633965Sjdp state[-1] = rand_type; 30733965Sjdp else 30833965Sjdp state[-1] = (MAX_TYPES * (rptr - state)) + rand_type; 30933965Sjdp 31033965Sjdp return ostate; 31133965Sjdp} 31233965Sjdp 31333965Sjdp/* Restore the state from the given state array. 31433965Sjdp Note: It is important that we also remember the locations of the pointers 31533965Sjdp in the current state information, and restore the locations of the pointers 31633965Sjdp from the old state information. This is done by multiplexing the pointer 31733965Sjdp location into the zeroeth word of the state information. Note that due 31833965Sjdp to the order in which things are done, it is OK to call setstate with the 31933965Sjdp same state as the current state 32033965Sjdp Returns a pointer to the old state information. */ 32133965Sjdp 32233965SjdpPTR 323218822Sdimsetstate (PTR arg_state) 32433965Sjdp{ 32533965Sjdp register long int *new_state = (long int *) arg_state; 32633965Sjdp register int type = new_state[0] % MAX_TYPES; 32733965Sjdp register int rear = new_state[0] / MAX_TYPES; 32833965Sjdp PTR ostate = (PTR) &state[-1]; 32933965Sjdp 33033965Sjdp if (rand_type == TYPE_0) 33133965Sjdp state[-1] = rand_type; 33233965Sjdp else 33333965Sjdp state[-1] = (MAX_TYPES * (rptr - state)) + rand_type; 33433965Sjdp 33533965Sjdp switch (type) 33633965Sjdp { 33733965Sjdp case TYPE_0: 33833965Sjdp case TYPE_1: 33933965Sjdp case TYPE_2: 34033965Sjdp case TYPE_3: 34133965Sjdp case TYPE_4: 34233965Sjdp rand_type = type; 34333965Sjdp rand_deg = degrees[type]; 34433965Sjdp rand_sep = seps[type]; 34533965Sjdp break; 34633965Sjdp default: 34733965Sjdp /* State info munged. */ 34833965Sjdp errno = EINVAL; 34933965Sjdp return NULL; 35033965Sjdp } 35133965Sjdp 35233965Sjdp state = &new_state[1]; 35333965Sjdp if (rand_type != TYPE_0) 35433965Sjdp { 35533965Sjdp rptr = &state[rear]; 35633965Sjdp fptr = &state[(rear + rand_sep) % rand_deg]; 35733965Sjdp } 35833965Sjdp /* Set end_ptr too. */ 35933965Sjdp end_ptr = &state[rand_deg]; 36033965Sjdp 36133965Sjdp return ostate; 36233965Sjdp} 36333965Sjdp 36433965Sjdp/* If we are using the trivial TYPE_0 R.N.G., just do the old linear 36533965Sjdp congruential bit. Otherwise, we do our fancy trinomial stuff, which is the 36633965Sjdp same in all ther other cases due to all the global variables that have been 36733965Sjdp set up. The basic operation is to add the number at the rear pointer into 36833965Sjdp the one at the front pointer. Then both pointers are advanced to the next 36933965Sjdp location cyclically in the table. The value returned is the sum generated, 37033965Sjdp reduced to 31 bits by throwing away the "least random" low bit. 37133965Sjdp Note: The code takes advantage of the fact that both the front and 37233965Sjdp rear pointers can't wrap on the same call by not testing the rear 37333965Sjdp pointer if the front one has wrapped. Returns a 31-bit random number. */ 37433965Sjdp 37533965Sjdplong int 376218822Sdimrandom (void) 37733965Sjdp{ 37833965Sjdp if (rand_type == TYPE_0) 37933965Sjdp { 38033965Sjdp state[0] = ((state[0] * 1103515245) + 12345) & LONG_MAX; 38133965Sjdp return state[0]; 38233965Sjdp } 38333965Sjdp else 38433965Sjdp { 38533965Sjdp long int i; 38633965Sjdp *fptr += *rptr; 38733965Sjdp /* Chucking least random bit. */ 38833965Sjdp i = (*fptr >> 1) & LONG_MAX; 38933965Sjdp ++fptr; 39033965Sjdp if (fptr >= end_ptr) 39133965Sjdp { 39233965Sjdp fptr = state; 39333965Sjdp ++rptr; 39433965Sjdp } 39533965Sjdp else 39633965Sjdp { 39733965Sjdp ++rptr; 39833965Sjdp if (rptr >= end_ptr) 39933965Sjdp rptr = state; 40033965Sjdp } 40133965Sjdp return i; 40233965Sjdp } 40333965Sjdp} 404