_A FAST PSEUDO RANDOM NUMBER GENERATOR_
by W.L. Maier

[LISTING ONE]

/******************************************************************************
*  Module:  r250.cpp   Description: implements R250 random number generator, 
*  from S. Kirkpatrick and E. Stoll, Journal of Computational Physics, 40,
*  p. 517 (1981). Written by:    W. L. Maier
******************************************************************************/

#include <stdlib.h>

/**** Static variables ****/
static unsigned int r250_buffer[250];
static int r250_index;

/**** Function prototypes  ****/
void r250_init(int seed);
unsigned int r250();
double dr250();

/**** Function: r250_init  Description: initializes r250 random number 
generator. ****/
void r250_init(int seed)
{
/*---------------------------------------------------------------------------*/
    int        j, k;
    unsigned int mask;
    unsigned int msb;
/*---------------------------------------------------------------------------*/
    srand(seed);
    r250_index = 0;
    for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
        r250_buffer[j] = rand();
    for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
        if (rand() > 16384)
            r250_buffer[j] |= 0x8000;
    msb = 0x8000;       /* To turn on the diagonal bit   */
    mask = 0xffff;      /* To turn off the leftmost bits */
    for (j = 0; j < 16; j++)
        {
        k = 11 * j + 3;             /* Select a word to operate on        */
        r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
        r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
        mask >>= 1;
        msb >>= 1;
        }
}
/**** Function:  r250  Description: returns a random unsigned integer. ****/
unsigned int r250()
{
/*---------------------------------------------------------------------------*/
    register int    j;
    register unsigned int new_rand;
/*---------------------------------------------------------------------------*/
    if (r250_index >= 147)
        j = r250_index - 147;      /* Wrap pointer around */
    else
        j = r250_index + 103;

    new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
    r250_buffer[r250_index] = new_rand;
    if (r250_index >= 249)      /* Increment pointer for next time */
        r250_index = 0;
    else
        r250_index++;

    return new_rand;
}
/**** Function:  r250  Description: returns a random double in range 0-1. ****/
double dr250()
{
/*---------------------------------------------------------------------------*/
    register int    j;
    register unsigned int new_rand;
/*---------------------------------------------------------------------------*/
    if (r250_index >= 147)
        j = r250_index - 147;     /* Wrap pointer around */
    else
        j = r250_index + 103;

    new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
    r250_buffer[r250_index] = new_rand;
    if (r250_index >= 249)      /* Increment pointer for next time */
        r250_index = 0;
    else
        r250_index++;
    return new_rand / (double)0xffff;   /* Return a number in 0.0 to 1.0 */
}


[LISTING TWO]

/******************************************************************************
*  Module: rtest.c   Description: tests R250 random number generator by 
*  placing data in a set of bins.
******************************************************************************/

#include <stdio.h>
#include <stdlib.h>

/**** Constants ****/
#define NMR_RAND	5000
#define MAX_BINS	500

/**** Function prototypes *****/
unsigned int r250();
void	r250_init(int seed);

/**** Function:	 main  ****/
void main(int argc, char *argv[])
{
/*---------------------------------------------------------------------------*/
	int		j, k;
	int		nmr_bins;
	int		seed;
	int		bins[MAX_BINS];
	double	randm;
	double	bin_limit[MAX_BINS];
	double	bin_inc;
/*---------------------------------------------------------------------------*/
	if (argc != 3)
		{
		printf("Usage -- rtest [nmr_bins] [seed]\n");
		exit(1);
		}
	nmr_bins = atoi(argv[1]);
	if (nmr_bins > MAX_BINS)
		{
		printf("Error -- maximum number of bins is %d\n", MAX_BINS);
		exit(1);
		}
	seed = atoi(argv[2]);
	r250_init(seed);
	bin_inc = 1.0 / nmr_bins;
	for (j = 0; j < nmr_bins; j++)
		{
		bins[j] = 0;		// Initialize bins to zero
		bin_limit[j] = (j + 1) * bin_inc;
		}
	bin_limit[nmr_bins-1] = 1.0e7;	// Make sure all others are in last bin
	for (j = 0; j < NMR_RAND; j++)
		{
		randm = r250() / (double)0xffff;
		for (k = 0; k < nmr_bins; k++)
			if (randm < bin_limit[k])
				{
				(bins[k])++;
				break;
				}
		}
	for (j = 0; j < nmr_bins; j++)
		printf("%d\n", bins[j]);
}