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