/contrib/ntp/libntp/ntp_random.c

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