Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!cornell!batcomputer!braner
From: braner@batcomputer.tn.cornell.edu (Moshe Braner)
Newsgroups: comp.lang.c
Subject: Yet another pseudo-random number generator
Summary: fast, passed all tests I submitted it to.
Message-ID: <7383@batcomputer.tn.cornell.edu>
Date: 10 Feb 89 20:58:39 GMT
Reply-To: braner@tcgould.tn.cornell.edu (Moshe Braner)
Organization: Cornell Theory Center, Cornell University, Ithaca NY
Lines: 94

/*
 *  Generate random integers using the shift-register algorithm.
 *  Copyright Moshe Braner,  Nov. 1987.  All rights reserved.
 *  Permission granted for noncommercial usage.
 *
 *  The algorithm here is basically  x[i] = x[i-31]  EOR x[i-7].
 *  (Could possibly improve by using x[i] = x[i-127] EOR x[i-63]?)
 *  Each x[i] here is 32 bits wide, could also be 16 or whatever.
 */

/* "#define FRAND" to include the floating-point stuff */

static long gla[31];			/* array of randoms to shift	*/
static int  glf=1;			/* seeding flag			*/
static long *glp, *glq, *gls, *glt;	/* pointers to shift register	*/
#ifdef FRAND
static double glfm;
#endif

/*
 * Seed ran5 - always do this first.
 * The seed is any positive integer.
 */
void
sran5(seed)
	int seed;
{
	register int i;
	register long x;

	x = (long) seed;
	for(i=0; i<100; i++)
		x = 31821*x+1;
	for(i=0; i<31; i++) {
		x = 31821*x+1;
		gla[i] = x;
	}
	glf = 0;
	gls = &gla[0];
	glp = &gla[0];
	glq = &gla[24];
	glt = &gla[31];
#ifdef FRAND
	glfm = 1.0;
	for (i=0; i<31; i++)
		glfm *= 0.5;		/* EXACTLY 2^(-31) */
#endif
}

/*
 * The basic ran5() generator, which returns
 * a long integer uniform between 0 and 2^31-1.
 */
long
ran5()
{
	register long x, y;

/* short form:

	x = ((*glp) ^ (*glq));
	if (++glp == glt)
		glp = gls;
	if (++glq == glt)
		glq = gls;
	return (x & 0x7FFFFFFF);

   safer form: */

	if (glf)
		sran5(0);
	y = x = ((*glp) ^ (*glq));	/* here FORTRAN dies */
	x <<= 1;
	if (y < 0)
		x += 1;		/* wish C had a 'rotate' operation */
	*glp = x;
	if (++glp == glt)
		glp = gls;
	if (++glq == glt)
		glq = gls;
	return (y & 0x7FFFFFFF);
}

#ifdef FRAND
/*
 * The floating-point ran5() generator.
 * Returns a double uniform on [0.0,1.0).
 */
double
fran5()
{
	return ((double)ran5() * glfm);
}
#endif FRAND