1/*
2 * Copyright (c) 1983, 1993
3 *	The Regents of the University of California.  All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 *    notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 *    notice, this list of conditions and the following disclaimer in the
12 *    documentation and/or other materials provided with the distribution.
13 * 4. Neither the name of the University nor the names of its contributors
14 *    may be used to endorse or promote products derived from this software
15 *    without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
18 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
21 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 * SUCH DAMAGE.
28 */
29
30#if defined(LIBC_SCCS) && !defined(lint)
31static char sccsid[] = "@(#)random.c	8.2 (Berkeley) 5/19/95";
32#endif /* LIBC_SCCS and not lint */
33#include <sys/cdefs.h>
34__FBSDID("$FreeBSD$");
35
36/* #include "namespace.h" */
37/* #include <sys/time.h>          /\* for srandomdev() *\/ */
38/* #include <fcntl.h>             /\* for srandomdev() *\/ */
39#include <stdint.h>
40#include <stdio.h>
41#include <stdlib.h>
42#include "posixcompat.h"
43/* #include <unistd.h>            /\* for srandomdev() *\/ */
44/* #include "un-namespace.h" */
45
46/*
47 * random.c:
48 *
49 * An improved random number generation package.  In addition to the standard
50 * rand()/srand() like interface, this package also has a special state info
51 * interface.  The initstate() routine is called with a seed, an array of
52 * bytes, and a count of how many bytes are being passed in; this array is
53 * then initialized to contain information for random number generation with
54 * that much state information.  Good sizes for the amount of state
55 * information are 32, 64, 128, and 256 bytes.  The state can be switched by
56 * calling the setstate() routine with the same array as was initiallized
57 * with initstate().  By default, the package runs with 128 bytes of state
58 * information and generates far better random numbers than a linear
59 * congruential generator.  If the amount of state information is less than
60 * 32 bytes, a simple linear congruential R.N.G. is used.
61 *
62 * Internally, the state information is treated as an array of uint32_t's; the
63 * zeroeth element of the array is the type of R.N.G. being used (small
64 * integer); the remainder of the array is the state information for the
65 * R.N.G.  Thus, 32 bytes of state information will give 7 ints worth of
66 * state information, which will allow a degree seven polynomial.  (Note:
67 * the zeroeth word of state information also has some other information
68 * stored in it -- see setstate() for details).
69 *
70 * The random number generation technique is a linear feedback shift register
71 * approach, employing trinomials (since there are fewer terms to sum up that
72 * way).  In this approach, the least significant bit of all the numbers in
73 * the state table will act as a linear feedback shift register, and will
74 * have period 2^deg - 1 (where deg is the degree of the polynomial being
75 * used, assuming that the polynomial is irreducible and primitive).  The
76 * higher order bits will have longer periods, since their values are also
77 * influenced by pseudo-random carries out of the lower bits.  The total
78 * period of the generator is approximately deg*(2**deg - 1); thus doubling
79 * the amount of state information has a vast influence on the period of the
80 * generator.  Note: the deg*(2**deg - 1) is an approximation only good for
81 * large deg, when the period of the shift is the dominant factor.
82 * With deg equal to seven, the period is actually much longer than the
83 * 7*(2**7 - 1) predicted by this formula.
84 *
85 * Modified 28 December 1994 by Jacob S. Rosenberg.
86 * The following changes have been made:
87 * All references to the type u_int have been changed to unsigned long.
88 * All references to type int have been changed to type long.  Other
89 * cleanups have been made as well.  A warning for both initstate and
90 * setstate has been inserted to the effect that on Sparc platforms
91 * the 'arg_state' variable must be forced to begin on word boundaries.
92 * This can be easily done by casting a long integer array to char *.
93 * The overall logic has been left STRICTLY alone.  This software was
94 * tested on both a VAX and Sun SpacsStation with exactly the same
95 * results.  The new version and the original give IDENTICAL results.
96 * The new version is somewhat faster than the original.  As the
97 * documentation says:  "By default, the package runs with 128 bytes of
98 * state information and generates far better random numbers than a linear
99 * congruential generator.  If the amount of state information is less than
100 * 32 bytes, a simple linear congruential R.N.G. is used."  For a buffer of
101 * 128 bytes, this new version runs about 19 percent faster and for a 16
102 * byte buffer it is about 5 percent faster.
103 */
104
105/*
106 * For each of the currently supported random number generators, we have a
107 * break value on the amount of state information (you need at least this
108 * many bytes of state info to support this random number generator), a degree
109 * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
110 * the separation between the two lower order coefficients of the trinomial.
111 */
112#define	TYPE_0		0		/* linear congruential */
113#define	BREAK_0		8
114#define	DEG_0		0
115#define	SEP_0		0
116
117#define	TYPE_1		1		/* x**7 + x**3 + 1 */
118#define	BREAK_1		32
119#define	DEG_1		7
120#define	SEP_1		3
121
122#define	TYPE_2		2		/* x**15 + x + 1 */
123#define	BREAK_2		64
124#define	DEG_2		15
125#define	SEP_2		1
126
127#define	TYPE_3		3		/* x**31 + x**3 + 1 */
128#define	BREAK_3		128
129#define	DEG_3		31
130#define	SEP_3		3
131
132#define	TYPE_4		4		/* x**63 + x + 1 */
133#define	BREAK_4		256
134#define	DEG_4		63
135#define	SEP_4		1
136
137/*
138 * Array versions of the above information to make code run faster --
139 * relies on fact that TYPE_i == i.
140 */
141#define	MAX_TYPES	5		/* max number of types above */
142
143#ifdef  USE_WEAK_SEEDING
144#define NSHUFF 0
145#else   /* !USE_WEAK_SEEDING */
146#define NSHUFF 50       /* to drop some "seed -> 1st value" linearity */
147#endif  /* !USE_WEAK_SEEDING */
148
149static const int degrees[MAX_TYPES] =	{ DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 };
150static const int seps [MAX_TYPES] =	{ SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
151
152/*
153 * Initially, everything is set up as if from:
154 *
155 *	initstate(1, randtbl, 128);
156 *
157 * Note that this initialization takes advantage of the fact that srandom()
158 * advances the front and rear pointers 10*rand_deg times, and hence the
159 * rear pointer which starts at 0 will also end up at zero; thus the zeroeth
160 * element of the state information, which contains info about the current
161 * position of the rear pointer is just
162 *
163 *	MAX_TYPES * (rptr - state) + TYPE_3 == TYPE_3.
164 */
165
166static uint32_t randtbl[DEG_3 + 1] = {
167	TYPE_3,
168#ifdef  USE_WEAK_SEEDING
169/* Historic implementation compatibility */
170/* The random sequences do not vary much with the seed */
171	0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342, 0xde3b81e0, 0xdf0a6fb5,
172	0xf103bc02, 0x48f340fb, 0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
173	0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86, 0xda672e2a, 0x1588ca88,
174	0xe369735d, 0x904f35f7, 0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
175	0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b, 0xf5ad9d0e, 0x8999220b,
176	0x27fb47b9,
177#else   /* !USE_WEAK_SEEDING */
178	0x991539b1, 0x16a5bce3, 0x6774a4cd, 0x3e01511e, 0x4e508aaa, 0x61048c05,
179	0xf5500617, 0x846b7115, 0x6a19892c, 0x896a97af, 0xdb48f936, 0x14898454,
180	0x37ffd106, 0xb58bff9c, 0x59e17104, 0xcf918a49, 0x09378c83, 0x52c7a471,
181	0x8d293ea9, 0x1f4fc301, 0xc3db71be, 0x39b44e1c, 0xf8a44ef9, 0x4c8b80b1,
182	0x19edc328, 0x87bf4bdd, 0xc9b240e5, 0xe9ee4b1b, 0x4382aee7, 0x535b6b41,
183	0xf3bec5da
184#endif  /* !USE_WEAK_SEEDING */
185};
186
187/*
188 * fptr and rptr are two pointers into the state info, a front and a rear
189 * pointer.  These two pointers are always rand_sep places aparts, as they
190 * cycle cyclically through the state information.  (Yes, this does mean we
191 * could get away with just one pointer, but the code for random() is more
192 * efficient this way).  The pointers are left positioned as they would be
193 * from the call
194 *
195 *	initstate(1, randtbl, 128);
196 *
197 * (The position of the rear pointer, rptr, is really 0 (as explained above
198 * in the initialization of randtbl) because the state table pointer is set
199 * to point to randtbl[1] (as explained below).
200 */
201static uint32_t *fptr = &randtbl[SEP_3 + 1];
202static uint32_t *rptr = &randtbl[1];
203
204/*
205 * The following things are the pointer to the state information table, the
206 * type of the current generator, the degree of the current polynomial being
207 * used, and the separation between the two pointers.  Note that for efficiency
208 * of random(), we remember the first location of the state information, not
209 * the zeroeth.  Hence it is valid to access state[-1], which is used to
210 * store the type of the R.N.G.  Also, we remember the last location, since
211 * this is more efficient than indexing every time to find the address of
212 * the last element to see if the front and rear pointers have wrapped.
213 */
214static uint32_t *state = &randtbl[1];
215static int rand_type = TYPE_3;
216static int rand_deg = DEG_3;
217static int rand_sep = SEP_3;
218static uint32_t *end_ptr = &randtbl[DEG_3 + 1];
219
220static inline uint32_t good_rand(int32_t);
221
222static inline uint32_t good_rand(int32_t x)
223{
224#ifdef  USE_WEAK_SEEDING
225/*
226 * Historic implementation compatibility.
227 * The random sequences do not vary much with the seed,
228 * even with overflowing.
229 */
230	return (1103515245 * x + 12345);
231#else   /* !USE_WEAK_SEEDING */
232/*
233 * Compute x = (7^5 * x) mod (2^31 - 1)
234 * wihout overflowing 31 bits:
235 *      (2^31 - 1) = 127773 * (7^5) + 2836
236 * From "Random number generators: good ones are hard to find",
237 * Park and Miller, Communications of the ACM, vol. 31, no. 10,
238 * October 1988, p. 1195.
239 */
240	int32_t hi, lo;
241
242	/* Can't be initialized with 0, so use another value. */
243	if (x == 0)
244		x = 123459876;
245	hi = x / 127773;
246	lo = x % 127773;
247	x = 16807 * lo - 2836 * hi;
248	if (x < 0)
249		x += 0x7fffffff;
250	return (x);
251#endif  /* !USE_WEAK_SEEDING */
252}
253
254/*
255 * srandom:
256 *
257 * Initialize the random number generator based on the given seed.  If the
258 * type is the trivial no-state-information type, just remember the seed.
259 * Otherwise, initializes state[] based on the given "seed" via a linear
260 * congruential generator.  Then, the pointers are set to known locations
261 * that are exactly rand_sep places apart.  Lastly, it cycles the state
262 * information a given number of times to get rid of any initial dependencies
263 * introduced by the L.C.R.N.G.  Note that the initialization of randtbl[]
264 * for default usage relies on values produced by this routine.
265 */
266void srandom(unsigned int x)
267{
268	int i, lim;
269
270	state[0] = (uint32_t)x;
271	if (rand_type == TYPE_0)
272		lim = NSHUFF;
273	else {
274		for (i = 1; i < rand_deg; i++)
275			state[i] = good_rand(state[i - 1]);
276		fptr = &state[rand_sep];
277		rptr = &state[0];
278		lim = 10 * rand_deg;
279	}
280	for (i = 0; i < lim; i++)
281		(void)random();
282}
283
284#if 0
285/*
286 * srandomdev:
287 *
288 * Many programs choose the seed value in a totally predictable manner.
289 * This often causes problems.  We seed the generator using the much more
290 * secure random(4) interface.  Note that this particular seeding
291 * procedure can generate states which are impossible to reproduce by
292 * calling srandom() with any value, since the succeeding terms in the
293 * state buffer are no longer derived from the LC algorithm applied to
294 * a fixed seed.
295 */
296void srandomdev(void)
297{
298	int fd, done;
299	size_t len;
300
301	if (rand_type == TYPE_0)
302		len = sizeof state[0];
303	else
304		len = rand_deg * sizeof state[0];
305
306	done = 0;
307	fd = _open("/dev/random", O_RDONLY, 0);
308	if (fd >= 0) {
309		if (_read(fd, (void *) state, len) == (ssize_t) len)
310			done = 1;
311		_close(fd);
312	}
313
314	if (!done) {
315		struct timeval tv;
316		unsigned long junk;
317
318		gettimeofday(&tv, NULL);
319		srandom((getpid() << 16) ^ tv.tv_sec ^ tv.tv_usec ^ junk);
320		return;
321	}
322
323	if (rand_type != TYPE_0) {
324		fptr = &state[rand_sep];
325		rptr = &state[0];
326	}
327}
328#endif
329
330/*
331 * initstate:
332 *
333 * Initialize the state information in the given array of n bytes for future
334 * random number generation.  Based on the number of bytes we are given, and
335 * the break values for the different R.N.G.'s, we choose the best (largest)
336 * one we can and set things up for it.  srandom() is then called to
337 * initialize the state information.
338 *
339 * Note that on return from srandom(), we set state[-1] to be the type
340 * multiplexed with the current value of the rear pointer; this is so
341 * successive calls to initstate() won't lose this information and will be
342 * able to restart with setstate().
343 *
344 * Note: the first thing we do is save the current state, if any, just like
345 * setstate() so that it doesn't matter when initstate is called.
346 *
347 * Returns a pointer to the old state.
348 *
349 * Note: The Sparc platform requires that arg_state begin on an int
350 * word boundary; otherwise a bus error will occur. Even so, lint will
351 * complain about mis-alignment, but you should disregard these messages.
352 */
353char *initstate(unsigned int seed, char *arg_state, size_t n)
354	/* unsigned long seed;		/\* seed for R.N.G. *\/ */
355	/* char *arg_state;		/\* pointer to state array *\/ */
356	/* long n;				/\* # bytes of state info *\/ */
357{
358	char *ostate = (char *)(&state[-1]);
359	uint32_t *int_arg_state = (uint32_t *)arg_state;
360
361	if (rand_type == TYPE_0)
362		state[-1] = rand_type;
363	else
364		state[-1] = MAX_TYPES * (rptr - state) + rand_type;
365	if (n < BREAK_0) {
366		(void)fprintf(stderr,
367		    "random: not enough state (%ld bytes); ignored.\n", n);
368		return(0);
369	}
370	if (n < BREAK_1) {
371		rand_type = TYPE_0;
372		rand_deg = DEG_0;
373		rand_sep = SEP_0;
374	} else if (n < BREAK_2) {
375		rand_type = TYPE_1;
376		rand_deg = DEG_1;
377		rand_sep = SEP_1;
378	} else if (n < BREAK_3) {
379		rand_type = TYPE_2;
380		rand_deg = DEG_2;
381		rand_sep = SEP_2;
382	} else if (n < BREAK_4) {
383		rand_type = TYPE_3;
384		rand_deg = DEG_3;
385		rand_sep = SEP_3;
386	} else {
387		rand_type = TYPE_4;
388		rand_deg = DEG_4;
389		rand_sep = SEP_4;
390	}
391	state = int_arg_state + 1; /* first location */
392	end_ptr = &state[rand_deg];	/* must set end_ptr before srandom */
393	srandom(seed);
394	if (rand_type == TYPE_0)
395		int_arg_state[0] = rand_type;
396	else
397		int_arg_state[0] = MAX_TYPES * (rptr - state) + rand_type;
398	return(ostate);
399}
400
401/*
402 * setstate:
403 *
404 * Restore the state from the given state array.
405 *
406 * Note: it is important that we also remember the locations of the pointers
407 * in the current state information, and restore the locations of the pointers
408 * from the old state information.  This is done by multiplexing the pointer
409 * location into the zeroeth word of the state information.
410 *
411 * Note that due to the order in which things are done, it is OK to call
412 * setstate() with the same state as the current state.
413 *
414 * Returns a pointer to the old state information.
415 *
416 * Note: The Sparc platform requires that arg_state begin on an int
417 * word boundary; otherwise a bus error will occur. Even so, lint will
418 * complain about mis-alignment, but you should disregard these messages.
419 */
420char *setstate(char *arg_state)
421	/* char *arg_state;		/\* pointer to state array *\/ */
422{
423	uint32_t *new_state = (uint32_t *)arg_state;
424	uint32_t type = new_state[0] % MAX_TYPES;
425	uint32_t rear = new_state[0] / MAX_TYPES;
426	char *ostate = (char *)(&state[-1]);
427
428	if (rand_type == TYPE_0)
429		state[-1] = rand_type;
430	else
431		state[-1] = MAX_TYPES * (rptr - state) + rand_type;
432	switch(type) {
433	case TYPE_0:
434	case TYPE_1:
435	case TYPE_2:
436	case TYPE_3:
437	case TYPE_4:
438		rand_type = type;
439		rand_deg = degrees[type];
440		rand_sep = seps[type];
441		break;
442	default:
443		(void)fprintf(stderr,
444		    "random: state info corrupted; not changed.\n");
445	}
446	state = new_state + 1;
447	if (rand_type != TYPE_0) {
448		rptr = &state[rear];
449		fptr = &state[(rear + rand_sep) % rand_deg];
450	}
451	end_ptr = &state[rand_deg];		/* set end_ptr too */
452	return(ostate);
453}
454
455/*
456 * random:
457 *
458 * If we are using the trivial TYPE_0 R.N.G., just do the old linear
459 * congruential bit.  Otherwise, we do our fancy trinomial stuff, which is
460 * the same in all the other cases due to all the global variables that have
461 * been set up.  The basic operation is to add the number at the rear pointer
462 * into the one at the front pointer.  Then both pointers are advanced to
463 * the next location cyclically in the table.  The value returned is the sum
464 * generated, reduced to 31 bits by throwing away the "least random" low bit.
465 *
466 * Note: the code takes advantage of the fact that both the front and
467 * rear pointers can't wrap on the same call by not testing the rear
468 * pointer if the front one has wrapped.
469 *
470 * Returns a 31-bit random number.
471 */
472long random(void)
473{
474	uint32_t i;
475	uint32_t *f, *r;
476
477	if (rand_type == TYPE_0) {
478		i = state[0];
479		state[0] = i = (good_rand(i)) & 0x7fffffff;
480	} else {
481		/*
482		 * Use local variables rather than static variables for speed.
483		 */
484		f = fptr; r = rptr;
485		*f += *r;
486		i = (*f >> 1) & 0x7fffffff;	/* chucking least random bit */
487		if (++f >= end_ptr) {
488			f = state;
489			++r;
490		}
491		else if (++r >= end_ptr) {
492			r = state;
493		}
494
495		fptr = f; rptr = r;
496	}
497	return((long)i);
498}
499