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