1*10d565efSmrg /*
2*10d565efSmrg  * Copyright (c) 1983 Regents of the University of California.
3*10d565efSmrg  * All rights reserved.
4*10d565efSmrg  *
5*10d565efSmrg  * Redistribution and use in source and binary forms, with or without
6*10d565efSmrg  * modification, are permitted provided that the following conditions
7*10d565efSmrg  * are met:
8*10d565efSmrg  * 1. Redistributions of source code must retain the above copyright
9*10d565efSmrg  *    notice, this list of conditions and the following disclaimer.
10*10d565efSmrg  * 2. Redistributions in binary form must reproduce the above copyright
11*10d565efSmrg  *    notice, this list of conditions and the following disclaimer in the
12*10d565efSmrg  *    documentation and/or other materials provided with the distribution.
13*10d565efSmrg  * 3. [rescinded 22 July 1999]
14*10d565efSmrg  * 4. Neither the name of the University nor the names of its contributors
15*10d565efSmrg  *    may be used to endorse or promote products derived from this software
16*10d565efSmrg  *    without specific prior written permission.
17*10d565efSmrg  *
18*10d565efSmrg  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19*10d565efSmrg  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20*10d565efSmrg  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21*10d565efSmrg  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22*10d565efSmrg  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23*10d565efSmrg  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24*10d565efSmrg  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25*10d565efSmrg  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26*10d565efSmrg  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27*10d565efSmrg  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28*10d565efSmrg  * SUCH DAMAGE.
29*10d565efSmrg  */
30*10d565efSmrg 
31*10d565efSmrg /*
32*10d565efSmrg  * This is derived from the Berkeley source:
33*10d565efSmrg  *	@(#)random.c	5.5 (Berkeley) 7/6/88
34*10d565efSmrg  * It was reworked for the GNU C Library by Roland McGrath.
35*10d565efSmrg  */
36*10d565efSmrg 
37*10d565efSmrg /*
38*10d565efSmrg 
39*10d565efSmrg @deftypefn Supplement {long int} random (void)
40*10d565efSmrg @deftypefnx Supplement void srandom (unsigned int @var{seed})
41*10d565efSmrg @deftypefnx Supplement void* initstate (unsigned int @var{seed}, @
42*10d565efSmrg   void *@var{arg_state}, unsigned long @var{n})
43*10d565efSmrg @deftypefnx Supplement void* setstate (void *@var{arg_state})
44*10d565efSmrg 
45*10d565efSmrg Random number functions.  @code{random} returns a random number in the
46*10d565efSmrg range 0 to @code{LONG_MAX}.  @code{srandom} initializes the random
47*10d565efSmrg number generator to some starting point determined by @var{seed}
48*10d565efSmrg (else, the values returned by @code{random} are always the same for each
49*10d565efSmrg run of the program).  @code{initstate} and @code{setstate} allow fine-grained
50*10d565efSmrg control over the state of the random number generator.
51*10d565efSmrg 
52*10d565efSmrg @end deftypefn
53*10d565efSmrg 
54*10d565efSmrg */
55*10d565efSmrg 
56*10d565efSmrg #include <errno.h>
57*10d565efSmrg 
58*10d565efSmrg #if 0
59*10d565efSmrg 
60*10d565efSmrg #include <ansidecl.h>
61*10d565efSmrg #include <limits.h>
62*10d565efSmrg #include <stddef.h>
63*10d565efSmrg #include <stdlib.h>
64*10d565efSmrg 
65*10d565efSmrg #else
66*10d565efSmrg 
67*10d565efSmrg #define	ULONG_MAX  ((unsigned long)(~0L))     /* 0xFFFFFFFF for 32-bits */
68*10d565efSmrg #define	LONG_MAX   ((long)(ULONG_MAX >> 1))   /* 0x7FFFFFFF for 32-bits*/
69*10d565efSmrg 
70*10d565efSmrg #ifdef __STDC__
71*10d565efSmrg #  define PTR void *
72*10d565efSmrg #  ifndef NULL
73*10d565efSmrg #    define NULL (void *) 0
74*10d565efSmrg #  endif
75*10d565efSmrg #else
76*10d565efSmrg #  define PTR char *
77*10d565efSmrg #  ifndef NULL
78*10d565efSmrg #    define NULL (void *) 0
79*10d565efSmrg #  endif
80*10d565efSmrg #endif
81*10d565efSmrg 
82*10d565efSmrg #endif
83*10d565efSmrg 
84*10d565efSmrg long int random (void);
85*10d565efSmrg 
86*10d565efSmrg /* An improved random number generation package.  In addition to the standard
87*10d565efSmrg    rand()/srand() like interface, this package also has a special state info
88*10d565efSmrg    interface.  The initstate() routine is called with a seed, an array of
89*10d565efSmrg    bytes, and a count of how many bytes are being passed in; this array is
90*10d565efSmrg    then initialized to contain information for random number generation with
91*10d565efSmrg    that much state information.  Good sizes for the amount of state
92*10d565efSmrg    information are 32, 64, 128, and 256 bytes.  The state can be switched by
93*10d565efSmrg    calling the setstate() function with the same array as was initiallized
94*10d565efSmrg    with initstate().  By default, the package runs with 128 bytes of state
95*10d565efSmrg    information and generates far better random numbers than a linear
96*10d565efSmrg    congruential generator.  If the amount of state information is less than
97*10d565efSmrg    32 bytes, a simple linear congruential R.N.G. is used.  Internally, the
98*10d565efSmrg    state information is treated as an array of longs; the zeroeth element of
99*10d565efSmrg    the array is the type of R.N.G. being used (small integer); the remainder
100*10d565efSmrg    of the array is the state information for the R.N.G.  Thus, 32 bytes of
101*10d565efSmrg    state information will give 7 longs worth of state information, which will
102*10d565efSmrg    allow a degree seven polynomial.  (Note: The zeroeth word of state
103*10d565efSmrg    information also has some other information stored in it; see setstate
104*10d565efSmrg    for details).  The random number generation technique is a linear feedback
105*10d565efSmrg    shift register approach, employing trinomials (since there are fewer terms
106*10d565efSmrg    to sum up that way).  In this approach, the least significant bit of all
107*10d565efSmrg    the numbers in the state table will act as a linear feedback shift register,
108*10d565efSmrg    and will have period 2^deg - 1 (where deg is the degree of the polynomial
109*10d565efSmrg    being used, assuming that the polynomial is irreducible and primitive).
110*10d565efSmrg    The higher order bits will have longer periods, since their values are
111*10d565efSmrg    also influenced by pseudo-random carries out of the lower bits.  The
112*10d565efSmrg    total period of the generator is approximately deg*(2**deg - 1); thus
113*10d565efSmrg    doubling the amount of state information has a vast influence on the
114*10d565efSmrg    period of the generator.  Note: The deg*(2**deg - 1) is an approximation
115*10d565efSmrg    only good for large deg, when the period of the shift register is the
116*10d565efSmrg    dominant factor.  With deg equal to seven, the period is actually much
117*10d565efSmrg    longer than the 7*(2**7 - 1) predicted by this formula.  */
118*10d565efSmrg 
119*10d565efSmrg 
120*10d565efSmrg 
121*10d565efSmrg /* For each of the currently supported random number generators, we have a
122*10d565efSmrg    break value on the amount of state information (you need at least thi
123*10d565efSmrg    bytes of state info to support this random number generator), a degree for
124*10d565efSmrg    the polynomial (actually a trinomial) that the R.N.G. is based on, and
125*10d565efSmrg    separation between the two lower order coefficients of the trinomial.  */
126*10d565efSmrg 
127*10d565efSmrg /* Linear congruential.  */
128*10d565efSmrg #define	TYPE_0		0
129*10d565efSmrg #define	BREAK_0		8
130*10d565efSmrg #define	DEG_0		0
131*10d565efSmrg #define	SEP_0		0
132*10d565efSmrg 
133*10d565efSmrg /* x**7 + x**3 + 1.  */
134*10d565efSmrg #define	TYPE_1		1
135*10d565efSmrg #define	BREAK_1		32
136*10d565efSmrg #define	DEG_1		7
137*10d565efSmrg #define	SEP_1		3
138*10d565efSmrg 
139*10d565efSmrg /* x**15 + x + 1.  */
140*10d565efSmrg #define	TYPE_2		2
141*10d565efSmrg #define	BREAK_2		64
142*10d565efSmrg #define	DEG_2		15
143*10d565efSmrg #define	SEP_2		1
144*10d565efSmrg 
145*10d565efSmrg /* x**31 + x**3 + 1.  */
146*10d565efSmrg #define	TYPE_3		3
147*10d565efSmrg #define	BREAK_3		128
148*10d565efSmrg #define	DEG_3		31
149*10d565efSmrg #define	SEP_3		3
150*10d565efSmrg 
151*10d565efSmrg /* x**63 + x + 1.  */
152*10d565efSmrg #define	TYPE_4		4
153*10d565efSmrg #define	BREAK_4		256
154*10d565efSmrg #define	DEG_4		63
155*10d565efSmrg #define	SEP_4		1
156*10d565efSmrg 
157*10d565efSmrg 
158*10d565efSmrg /* Array versions of the above information to make code run faster.
159*10d565efSmrg    Relies on fact that TYPE_i == i.  */
160*10d565efSmrg 
161*10d565efSmrg #define	MAX_TYPES	5	/* Max number of types above.  */
162*10d565efSmrg 
163*10d565efSmrg static int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 };
164*10d565efSmrg static int seps[MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
165*10d565efSmrg 
166*10d565efSmrg 
167*10d565efSmrg 
168*10d565efSmrg /* Initially, everything is set up as if from:
169*10d565efSmrg 	initstate(1, randtbl, 128);
170*10d565efSmrg    Note that this initialization takes advantage of the fact that srandom
171*10d565efSmrg    advances the front and rear pointers 10*rand_deg times, and hence the
172*10d565efSmrg    rear pointer which starts at 0 will also end up at zero; thus the zeroeth
173*10d565efSmrg    element of the state information, which contains info about the current
174*10d565efSmrg    position of the rear pointer is just
175*10d565efSmrg 	(MAX_TYPES * (rptr - state)) + TYPE_3 == TYPE_3.  */
176*10d565efSmrg 
177*10d565efSmrg static long int randtbl[DEG_3 + 1] =
178*10d565efSmrg   { TYPE_3,
179*10d565efSmrg       0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342,
180*10d565efSmrg       0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb,
181*10d565efSmrg       0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
182*10d565efSmrg       0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86,
183*10d565efSmrg       0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7,
184*10d565efSmrg       0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
185*10d565efSmrg       0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b,
186*10d565efSmrg       0xf5ad9d0e, 0x8999220b, 0x27fb47b9
187*10d565efSmrg     };
188*10d565efSmrg 
189*10d565efSmrg /* FPTR and RPTR are two pointers into the state info, a front and a rear
190*10d565efSmrg    pointer.  These two pointers are always rand_sep places aparts, as they
191*10d565efSmrg    cycle through the state information.  (Yes, this does mean we could get
192*10d565efSmrg    away with just one pointer, but the code for random is more efficient
193*10d565efSmrg    this way).  The pointers are left positioned as they would be from the call:
194*10d565efSmrg 	initstate(1, randtbl, 128);
195*10d565efSmrg    (The position of the rear pointer, rptr, is really 0 (as explained above
196*10d565efSmrg    in the initialization of randtbl) because the state table pointer is set
197*10d565efSmrg    to point to randtbl[1] (as explained below).)  */
198*10d565efSmrg 
199*10d565efSmrg static long int *fptr = &randtbl[SEP_3 + 1];
200*10d565efSmrg static long int *rptr = &randtbl[1];
201*10d565efSmrg 
202*10d565efSmrg 
203*10d565efSmrg 
204*10d565efSmrg /* The following things are the pointer to the state information table,
205*10d565efSmrg    the type of the current generator, the degree of the current polynomial
206*10d565efSmrg    being used, and the separation between the two pointers.
207*10d565efSmrg    Note that for efficiency of random, we remember the first location of
208*10d565efSmrg    the state information, not the zeroeth.  Hence it is valid to access
209*10d565efSmrg    state[-1], which is used to store the type of the R.N.G.
210*10d565efSmrg    Also, we remember the last location, since this is more efficient than
211*10d565efSmrg    indexing every time to find the address of the last element to see if
212*10d565efSmrg    the front and rear pointers have wrapped.  */
213*10d565efSmrg 
214*10d565efSmrg static long int *state = &randtbl[1];
215*10d565efSmrg 
216*10d565efSmrg static int rand_type = TYPE_3;
217*10d565efSmrg static int rand_deg = DEG_3;
218*10d565efSmrg static int rand_sep = SEP_3;
219*10d565efSmrg 
220*10d565efSmrg static long int *end_ptr = &randtbl[sizeof(randtbl) / sizeof(randtbl[0])];
221*10d565efSmrg 
222*10d565efSmrg /* Initialize the random number generator based on the given seed.  If the
223*10d565efSmrg    type is the trivial no-state-information type, just remember the seed.
224*10d565efSmrg    Otherwise, initializes state[] based on the given "seed" via a linear
225*10d565efSmrg    congruential generator.  Then, the pointers are set to known locations
226*10d565efSmrg    that are exactly rand_sep places apart.  Lastly, it cycles the state
227*10d565efSmrg    information a given number of times to get rid of any initial dependencies
228*10d565efSmrg    introduced by the L.C.R.N.G.  Note that the initialization of randtbl[]
229*10d565efSmrg    for default usage relies on values produced by this routine.  */
230*10d565efSmrg void
srandom(unsigned int x)231*10d565efSmrg srandom (unsigned int x)
232*10d565efSmrg {
233*10d565efSmrg   state[0] = x;
234*10d565efSmrg   if (rand_type != TYPE_0)
235*10d565efSmrg     {
236*10d565efSmrg       register long int i;
237*10d565efSmrg       for (i = 1; i < rand_deg; ++i)
238*10d565efSmrg 	state[i] = (1103515145 * state[i - 1]) + 12345;
239*10d565efSmrg       fptr = &state[rand_sep];
240*10d565efSmrg       rptr = &state[0];
241*10d565efSmrg       for (i = 0; i < 10 * rand_deg; ++i)
242*10d565efSmrg 	random();
243*10d565efSmrg     }
244*10d565efSmrg }
245*10d565efSmrg 
246*10d565efSmrg /* Initialize the state information in the given array of N bytes for
247*10d565efSmrg    future random number generation.  Based on the number of bytes we
248*10d565efSmrg    are given, and the break values for the different R.N.G.'s, we choose
249*10d565efSmrg    the best (largest) one we can and set things up for it.  srandom is
250*10d565efSmrg    then called to initialize the state information.  Note that on return
251*10d565efSmrg    from srandom, we set state[-1] to be the type multiplexed with the current
252*10d565efSmrg    value of the rear pointer; this is so successive calls to initstate won't
253*10d565efSmrg    lose this information and will be able to restart with setstate.
254*10d565efSmrg    Note: The first thing we do is save the current state, if any, just like
255*10d565efSmrg    setstate so that it doesn't matter when initstate is called.
256*10d565efSmrg    Returns a pointer to the old state.  */
257*10d565efSmrg PTR
initstate(unsigned int seed,PTR arg_state,unsigned long n)258*10d565efSmrg initstate (unsigned int seed, PTR arg_state, unsigned long n)
259*10d565efSmrg {
260*10d565efSmrg   PTR ostate = (PTR) &state[-1];
261*10d565efSmrg 
262*10d565efSmrg   if (rand_type == TYPE_0)
263*10d565efSmrg     state[-1] = rand_type;
264*10d565efSmrg   else
265*10d565efSmrg     state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
266*10d565efSmrg   if (n < BREAK_1)
267*10d565efSmrg     {
268*10d565efSmrg       if (n < BREAK_0)
269*10d565efSmrg 	{
270*10d565efSmrg 	  errno = EINVAL;
271*10d565efSmrg 	  return NULL;
272*10d565efSmrg 	}
273*10d565efSmrg       rand_type = TYPE_0;
274*10d565efSmrg       rand_deg = DEG_0;
275*10d565efSmrg       rand_sep = SEP_0;
276*10d565efSmrg     }
277*10d565efSmrg   else if (n < BREAK_2)
278*10d565efSmrg     {
279*10d565efSmrg       rand_type = TYPE_1;
280*10d565efSmrg       rand_deg = DEG_1;
281*10d565efSmrg       rand_sep = SEP_1;
282*10d565efSmrg     }
283*10d565efSmrg   else if (n < BREAK_3)
284*10d565efSmrg     {
285*10d565efSmrg       rand_type = TYPE_2;
286*10d565efSmrg       rand_deg = DEG_2;
287*10d565efSmrg       rand_sep = SEP_2;
288*10d565efSmrg     }
289*10d565efSmrg   else if (n < BREAK_4)
290*10d565efSmrg     {
291*10d565efSmrg       rand_type = TYPE_3;
292*10d565efSmrg       rand_deg = DEG_3;
293*10d565efSmrg       rand_sep = SEP_3;
294*10d565efSmrg     }
295*10d565efSmrg   else
296*10d565efSmrg     {
297*10d565efSmrg       rand_type = TYPE_4;
298*10d565efSmrg       rand_deg = DEG_4;
299*10d565efSmrg       rand_sep = SEP_4;
300*10d565efSmrg     }
301*10d565efSmrg 
302*10d565efSmrg   state = &((long int *) arg_state)[1];	/* First location.  */
303*10d565efSmrg   /* Must set END_PTR before srandom.  */
304*10d565efSmrg   end_ptr = &state[rand_deg];
305*10d565efSmrg   srandom(seed);
306*10d565efSmrg   if (rand_type == TYPE_0)
307*10d565efSmrg     state[-1] = rand_type;
308*10d565efSmrg   else
309*10d565efSmrg     state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
310*10d565efSmrg 
311*10d565efSmrg   return ostate;
312*10d565efSmrg }
313*10d565efSmrg 
314*10d565efSmrg /* Restore the state from the given state array.
315*10d565efSmrg    Note: It is important that we also remember the locations of the pointers
316*10d565efSmrg    in the current state information, and restore the locations of the pointers
317*10d565efSmrg    from the old state information.  This is done by multiplexing the pointer
318*10d565efSmrg    location into the zeroeth word of the state information. Note that due
319*10d565efSmrg    to the order in which things are done, it is OK to call setstate with the
320*10d565efSmrg    same state as the current state
321*10d565efSmrg    Returns a pointer to the old state information.  */
322*10d565efSmrg 
323*10d565efSmrg PTR
setstate(PTR arg_state)324*10d565efSmrg setstate (PTR arg_state)
325*10d565efSmrg {
326*10d565efSmrg   register long int *new_state = (long int *) arg_state;
327*10d565efSmrg   register int type = new_state[0] % MAX_TYPES;
328*10d565efSmrg   register int rear = new_state[0] / MAX_TYPES;
329*10d565efSmrg   PTR ostate = (PTR) &state[-1];
330*10d565efSmrg 
331*10d565efSmrg   if (rand_type == TYPE_0)
332*10d565efSmrg     state[-1] = rand_type;
333*10d565efSmrg   else
334*10d565efSmrg     state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
335*10d565efSmrg 
336*10d565efSmrg   switch (type)
337*10d565efSmrg     {
338*10d565efSmrg     case TYPE_0:
339*10d565efSmrg     case TYPE_1:
340*10d565efSmrg     case TYPE_2:
341*10d565efSmrg     case TYPE_3:
342*10d565efSmrg     case TYPE_4:
343*10d565efSmrg       rand_type = type;
344*10d565efSmrg       rand_deg = degrees[type];
345*10d565efSmrg       rand_sep = seps[type];
346*10d565efSmrg       break;
347*10d565efSmrg     default:
348*10d565efSmrg       /* State info munged.  */
349*10d565efSmrg       errno = EINVAL;
350*10d565efSmrg       return NULL;
351*10d565efSmrg     }
352*10d565efSmrg 
353*10d565efSmrg   state = &new_state[1];
354*10d565efSmrg   if (rand_type != TYPE_0)
355*10d565efSmrg     {
356*10d565efSmrg       rptr = &state[rear];
357*10d565efSmrg       fptr = &state[(rear + rand_sep) % rand_deg];
358*10d565efSmrg     }
359*10d565efSmrg   /* Set end_ptr too.  */
360*10d565efSmrg   end_ptr = &state[rand_deg];
361*10d565efSmrg 
362*10d565efSmrg   return ostate;
363*10d565efSmrg }
364*10d565efSmrg 
365*10d565efSmrg /* If we are using the trivial TYPE_0 R.N.G., just do the old linear
366*10d565efSmrg    congruential bit.  Otherwise, we do our fancy trinomial stuff, which is the
367*10d565efSmrg    same in all ther other cases due to all the global variables that have been
368*10d565efSmrg    set up.  The basic operation is to add the number at the rear pointer into
369*10d565efSmrg    the one at the front pointer.  Then both pointers are advanced to the next
370*10d565efSmrg    location cyclically in the table.  The value returned is the sum generated,
371*10d565efSmrg    reduced to 31 bits by throwing away the "least random" low bit.
372*10d565efSmrg    Note: The code takes advantage of the fact that both the front and
373*10d565efSmrg    rear pointers can't wrap on the same call by not testing the rear
374*10d565efSmrg    pointer if the front one has wrapped.  Returns a 31-bit random number.  */
375*10d565efSmrg 
376*10d565efSmrg long int
random(void)377*10d565efSmrg random (void)
378*10d565efSmrg {
379*10d565efSmrg   if (rand_type == TYPE_0)
380*10d565efSmrg     {
381*10d565efSmrg       state[0] = ((state[0] * 1103515245) + 12345) & LONG_MAX;
382*10d565efSmrg       return state[0];
383*10d565efSmrg     }
384*10d565efSmrg   else
385*10d565efSmrg     {
386*10d565efSmrg       long int i;
387*10d565efSmrg       *fptr += *rptr;
388*10d565efSmrg       /* Chucking least random bit.  */
389*10d565efSmrg       i = (*fptr >> 1) & LONG_MAX;
390*10d565efSmrg       ++fptr;
391*10d565efSmrg       if (fptr >= end_ptr)
392*10d565efSmrg 	{
393*10d565efSmrg 	  fptr = state;
394*10d565efSmrg 	  ++rptr;
395*10d565efSmrg 	}
396*10d565efSmrg       else
397*10d565efSmrg 	{
398*10d565efSmrg 	  ++rptr;
399*10d565efSmrg 	  if (rptr >= end_ptr)
400*10d565efSmrg 	    rptr = state;
401*10d565efSmrg 	}
402*10d565efSmrg       return i;
403*10d565efSmrg     }
404*10d565efSmrg }
405