Bug fix for 4.2bsd random.c

[D. Mitchell]

1,348c
/*
 *	Additive Random Number Generator  (see: Knuth 3.2.2)
 *	D. P. Mitchell  84/04/16.
 */
#define LENGTH	97
#define TAP	12
#define MANTISSA 56

static long    y[LENGTH],  *tap =  &y[0],  *feedback =  &y[LENGTH - TAP];
static double fy[LENGTH], *ftap = &fy[0], *ffeedback = &fy[LENGTH - TAP];
static unsigned long x = 0xc0393e8d;

random()
{
	if (--tap < y) {
		if (x)		/* effect of srandom might be delayed */
			setup();
		tap += LENGTH;
	}
	if (--feedback < y)
		feedback += LENGTH;
	return (*feedback += *tap) >> 1 & 0x7fffffff;
}

double
frandom()
{
	double sum;

	if (--ftap < fy) {
		if (x)		/* fy cannot be preset portably! */
			setup();
		ftap += LENGTH;
	}
	if (--ffeedback < fy)
		ffeedback += LENGTH;
	sum = *ftap + *ffeedback;
	if (sum >= 1.0) {
		sum = *ftap - 0.5;	/* must not overflow! */
		sum += *ffeedback - 0.5;
	}
	*ffeedback = sum;
	return sum;
}

srandom(seed)
long seed;
{
	x = seed;
}

static
setup()
{
	long mask;
	double divisor;
	int i;

	for (i = 0; i < LENGTH; i++) {
		x = 177988*x*x + 505360173*x + 907633385;
		y[i]  = x;
	}
	mask = 0x7fffffff >> (63 - MANTISSA);
	divisor = (double)mask + 1.0;
	for (i = 0; i < LENGTH; i++) {
		x = 177988*x*x + 505360173*x + 907633385;
		fy[i] = (double)x / 4294967296.0;
		x = 177988*x*x + 505360173*x + 907633385;
		fy[i] = ((double)(x & mask) + fy[i]) / divisor;
	}
	x = 0;
}
.