#include <stdlib.h>

/*
32bit unsigned is asssumed
traditional bsd random implementation with park-miller seeding
lagged fibonacci generator params: (7,3), (15,1), (31,3), (63,1)
does not have good statistical properties
*/

static unsigned init[] = {
0x00000000,0x991539b1,0x16a5bce3,0x6774a4cd,
0x3e01511e,0x4e508aaa,0x61048c05,0xf5500617,
0x846b7115,0x6a19892c,0x896a97af,0xdb48f936,
0x14898454,0x37ffd106,0xb58bff9c,0x59e17104,
0xcf918a49,0x09378c83,0x52c7a471,0x8d293ea9,
0x1f4fc301,0xc3db71be,0x39b44e1c,0xf8a44ef9,
0x4c8b80b1,0x19edc328,0x87bf4bdd,0xc9b240e5,
0xe9ee4b1b,0x4382aee7,0x535b6b41,0xf3bec5da};

static int n = 31;
static int i = 3;
static int j = 0;
static unsigned *x = init+1;

static unsigned lcg(unsigned x) {
	return (1103515245*x + 12345) & 0x7fffffff;
}

static unsigned parkmiller(unsigned x) {
	/* x = 7^5 x mod (2^31-1) */
	x = 16807*(x%127773) - 2836*(x/127773);
	if (x & 0x80000000)
		x += 0x7fffffff;
	return x;
}

static void *savestate() {
	x[-1] = ((unsigned)n<<16)|(i<<8)|j;
	return x-1;
}

static void loadstate(unsigned *state) {
	x = state+1;
	n = x[-1]>>16;
	i = (x[-1]>>8)&0xff;
	j = x[-1]&0xff;
}

void srandom(unsigned seed) {
	int k;

	i = n == 31 || n == 7 ? 3 : 1;
	j = 0;
	x[0] = seed ? seed : 1;
	for (k = 1; k < n; k++)
		x[k] = parkmiller(x[k-1]);
	for (k = 0; k < 10*n; k++)
		random();
}

char *initstate(unsigned seed, char *state, size_t size) {
	void *old = savestate();
	if (size < 8)
		return 0;
	else if (size < 32)
		n = 0;
	else if (size < 64)
		n = 7;
	else if (size < 128)
		n = 15;
	else if (size < 256)
		n = 31;
	else
		n = 63;
	x = (unsigned*)state + 1;
	srandom(seed);
	return old;
}

char *setstate(char *state) {
	void *old = savestate();
	loadstate((unsigned*)state);
	return old;
}

long random(void) {
	unsigned k;

	if (n == 0)
		return x[0] = lcg(x[0]);
	x[i] += x[j];
	k = x[i]>>1;
	if (++i == n)
		i = 0;
	if (++j == n)
		j = 0;
	return k;
}