Random number generators

[Steve Kargl]

On Tue, Mar 17, 2015 at 07:10:46AM -0500, Pedro Giffuni wrote:
> 
> > One big issue is saving internal state.  IIRC, MT requires 623-bit
> > of internal state.  KISS requires 4 32-bit int.  Thus, if
> > you want to reseed the generator, KISS requires far less effort.
> > 
> 
> Yes, the problem is the that libc requires a single 32 (or 31) bit seed.
> Given that restriction, our existing generator is not bad. Enforcing
> something better breaks the API and is not really practical to get
> crypto-grade randomness for stuff like refreshing a slide in a
> presentation anyways.
> 
> The musl libc approach seemed reasonable but I haven?t looked at the
> base random generator there (I?ve heard the glibc one is not good at all).
> 

GM's kiss generator has a period of 2**123.  The manpage
for random(3) claims a period of about 2**35.  The rand(3)
manpage does not give a period, but I suspect it is well
short of 2**123.

The code that follows my sig uses. a Lehmer LCG generator to
provide the 4 seeds needed for kiss.  The Lehmer LCG takes a
single 32 bit seed.  If reseed() is not called prior to calling
kiss(), then a default set of seeds is used.  If the argument
to reseed() is 0 or 1, then the default set of seeds is used.
It should be straight forward to map reseed() to srand() and
kiss() to rand().  Do note that range kiss() is (0,UINT_MAX],
so one needs to subtract 1 to get [0,UINT_MAX-1) if 0 needs 
to be in the range.

To use this code in preference to random(3) (and in violation of
POSIX?), initstate() and setstate() would need to muck with the
internal state of kiss().  This is certainly possible, but I
don't do it below.

-- 
Steve

#include <sys/types.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <stdint.h>
#include <stdio.h>

/*
 * The KISS generator requires 4 seeds.  Use a simple Lehmer linear
 * congruential generator (LCG) that requires only a single seed to
 * generate the 4 KISS seeds.
 */
static uint64_t lcg_seed = 0;	/* Need to remember the state of LCG. */

static uint32_t
lehmer_lcg(uint32_t a)
{
	return ((uint64_t)a * 279470273UL) % 4294967291UL;
}

/*
 * If reseed() is not called or it is called with an argument of 0 or 1,
 * then use DEFAULT_SEED as the set of seeds for kiss().
 */
#define	DEFAULT_SEED	{123456789, 362436069, 521288629, 916191069}

static const uint32_t default_seed[4] = DEFAULT_SEED;
static uint32_t seed[4] = DEFAULT_SEED;

uint32_t
reseed(uint32_t s)
{
	int i;
	uint32_t kiss_seed;

	/* Copy seed given by user. */
	lcg_seed = s;
	if (lcg_seed == 0 || lcg_seed == 1) {
		for (i = 0; i < 4; i++)
			seed[i] = default_seed[i];
	} else {
		kiss_seed = lcg_seed;
		for (i = 0; i < 4; i++) {
			kiss_seed = lehmer_lcg(kiss_seed);
			seed[i] = kiss_seed;
		}
	}
	return (lcg_seed);
}

/*
 * kiss() returns an integer value in the range of (0, 2**32-1].  The
 * distribution of pseudorandom numbers should be uniform.  The overall
 * perdiod for this generator exceeds 2**123.
 */

#define	LEFT(k, n)	((k)^((k)<<(n)))
#define	RGHT(k, n)	((k)^((k)>>(n)))

uint32_t
kiss(void)
{

	seed[0] = 69069 * seed[0] + 1327217885;
	seed[1] = LEFT(RGHT(LEFT(seed[1],13),17),5);
	seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
	seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
	return (seed[0] + seed[1] + (seed[2] << 16) + seed[3]);
}

#define	NUM	5

int
main(void)
{
	int i;
	uint32_t n;
	struct rusage t0, t1;

	/* Default seeds. */
	printf("Default seeds used.\n");
	for (i = 0; i < NUM; i++)
		printf("%u\n", kiss());

	n = reseed(42);
	printf("\nseed = %u\n", n);
	for (i = 0; i < NUM; i++)
		printf("%u\n", kiss());

	/* Default seeds. */
	n = reseed(0);
	printf("\nseed = %u\n", n);
	for (i = 0; i < NUM; i++)
		printf("%u\n", kiss());

	n = reseed(12345);
	printf("\nseed = %u\n", n);
	for (i = 0; i < NUM; i++)
		printf("%u\n", kiss());

	/* Default seeds. */
	n = reseed(1);
	printf("\nseed = %u\n", n);
	for (i = 0; i < NUM; i++)
		printf("%u\n", kiss());

	getrusage(RUSAGE_SELF, &t0);
	for (i = 0; i < 1000000; i++) n = kiss();
	getrusage(RUSAGE_SELF, &t1);

	{
        	double u;

        	u = (t1.ru_utime.tv_sec - t0.ru_utime.tv_sec);
	        u += (t1.ru_stime.tv_sec - t0.ru_stime.tv_sec);
        	u += 1.e-6*(t1.ru_utime.tv_usec - t0.ru_utime.tv_usec);
        	u += 1.e-6*(t1.ru_stime.tv_usec - t0.ru_stime.tv_usec);

       		printf("\n%u\n", n);
       		printf("%.2f ms for 1 million values\n", n, 1000 * u);
	}

	return 0;
}