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


//
// Note:
//   This file is slightly modified by Jim Steuert from
//   the original by W. L. Maier
//

/******************************************************************************
*  Module:  r250.c   
*  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
*
*	METHOD...
*		16 parallel copies of a linear shift register with
*		period 2^250 - 1.  FAR longer period than the usual
*		linear congruent generator, and commonly faster as
*		well.  (For details see the above paper, and the
*		article in DDJ referenced below.)
*
*	HISTORY...
*		Sep 92: Number returned by dr250() is in range [0,1) instead 
*			of [0,1], so for example a random angle in the
*			interval [0, 2*PI) can be calculated
*			conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
*		Aug 92: Initialization is optional.  Default condition is 
*			equivalent to initializing with the seed 12345,
*			so that the sequence of random numbers begins:
*			1173, 53403, 52352, 35341...  (9996 numbers
*			skipped) ...57769, 14511, 46930, 11942, 7978,
*			56163, 46506, 45768, 21162, 43113...  Using ^=
*			operator rather than two separate statements. 
*			Initializing with own linear congruent
*			pseudorandom number generator for portability. 
*			Function prototypes moved to a header file. 
*			Implemented r250n() to generate numbers
*			uniformly distributed in a specific range
*			[0,n), because the common expedient rand()%n is
*			incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
*		May 91: Published by W. L. Maier, "A Fast Pseudo Random Number 
*			Generator", Dr. Dobb's Journal #176.
******************************************************************************/

long r250_index = 0;
unsigned long r250_buffer[250] = { 0, 3, };
unsigned long gblseed = 1;
unsigned long xseed = 1;


/*================================================*/

unsigned long myrand()
{
   unsigned long rval;
   gblseed = gblseed*0x015a4e35L + 1;
   rval = ((gblseed>>16) & 0x7fff);
   return( rval );
}

/*================================================*/
void r250_init(long seed)
{
    long  j, k;
    long  i,x;
    unsigned long mask;
    unsigned long msb;
	unsigned long left;
	unsigned long right;

    gblseed = seed;
    r250_index = 0;

    for (j = 0; j < 250; j++)
      r250_buffer[j] = myrand();

    for (j = 0; j < 250; j++)
    {
		left = myrand();
		right = 16384;
        if ( left > right )
            r250_buffer[j] |= 0x8000;
    }
    msb = 0x8000;
    mask = 0xffff;

    for (j = 0; j < 16; j++)
    {
      k = (11*j) + 3;
      r250_buffer[k] &= mask;
      r250_buffer[k] |= msb;
      mask >>= 1;
      msb >>= 1;
    }
}

/*================================================*/

unsigned long r250()
{
    long    j;
    unsigned long new_rand;
    unsigned long a,b,c;

    if (r250_index >= 147)
        j = r250_index - 147;
    else
        j = r250_index + 103;
    new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
    r250_buffer[r250_index] = new_rand;
    if (r250_index >= 249)
        r250_index = 0;
    else
        r250_index++;
    return new_rand;
}

/*================================================*/

unsigned long r250n(unsigned long n)
{
    long    j;
    unsigned long new_rand, limit;
    unsigned long retval;
    unsigned long a,b,c;

    limit = (65535U/n)*n;
    do 
    {
        new_rand = r250();
        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++;
    }
    while(new_rand >= limit);
    retval = new_rand % n;
    return( retval );
}