/* random number generators
 *	Contibuted by Jim Wilson (wilson@ucbarpa.berkeley.edu)
 *	November 22, 1988
 */

/* all algorithms use SYSV calling conventions for ease of use
   i.e. lrand48(), seed48(), srand48()
 */

/* Source for two generators have been supplied.  Define the one you want */
/* #define Alg_AS_183 */
#define Minimal_Standard

/* Define this to compile as a standalone test */
/* #define TEST_RNG */


/* both Alg_AS_183 and Minimal_Standard use prime modulus multiplicative
   congruential generators (PMMLCG), also known as Lehmer Grammers, which
   satisfy the following properties 

   (i)   modulus: m - a large prime integer
   (ii)  multiplier: a - an integer in the range 2, 3, ..., m - 1
   (iii) z[n+1] = f(z[n]), for n = 1, 2, ...
   (iv)  f(z) = az mod m
   (v)   u[n] = z[n] / m, for n = 1, 2, ...

   The sequence of z's must be initialized by choosing an initial seed
   z[1] from the range 1, 2, ..., m - 1.  The sequence of z's is a pseudo-
   random sequence drawn without replacement from the set 1, 2, ..., m - 1.
   The u's form a psuedo-random sequence of real numbers between (but not
   including) 0 and 1.

   Both algorithms use Schrage's method to compute the sequence of z's.
   Let m = aq + r, where q = m div a, and r = m mod a.
   Then f(z) = az mod m = az - m * (az div m) =
             = gamma(z) + m * delta(z)
   Where gamma(z) = a(z mod q) - r(z div q)
   and   delta(z) = (z div q) - (az div m)

   If r < q, then for all z in 1, 2, ..., m - 1:
   (1) delta(z) is either 0 or 1
   (2) both a(z mod q) and r(z div q) are in 0, 1, ..., m - 1
   (3) absolute value of gamma(z) <= m - 1
   (4) delta(z) = 1 iff gamma(z) < 0

   Hence each value of z can be computed exactly without overflow as long
   as m can be represented as an integer.
 */


#ifdef Minimal_Standard

/* a good random number generator, correct on any machine with 32 bit integers,
   this algorithm is from:

Stephen K. Park and Keith W. Miller, "Random Number Generators:
	Good ones are hard to find", Communications of the ACM, October 1988,
	vol 31, number 10, pp. 1192-1201.

   If this algorithm is implemented correctly, then if z[1] = 1, then
   z[10001] will equal 1043618065

   Has a full period of 2^31 - 1.
   Returns integers in the range 1 to 2^31-1.
 */

#define RNG_M 2147483647L  /* m = 2^31 - 1 */
#define RNG_A 16807L
#define RNG_Q 127773L      /* m div a */
#define RNG_R 2836L        /* m mod a */

/* 32 bit seed */
static long xsubi;

unsigned short *
seed48 (seed16v)
unsigned short seed16v[3];
{
  static unsigned short old_xsubi[3];
  union seed { unsigned short s[3]; long l; } seed;
  int i;

  seed.l = xsubi;

  for (i = 0; i < 3; i++)
    {
      old_xsubi[i] = seed.s[i];
      seed.s[i] = seed16v[i];
    }

  xsubi = seed.l;

  return (old_xsubi);
}

void
srand48 (seedval)
unsigned long seedval;
{
  /* set seed to value between 1 and m-1 */

  xsubi = (seedval % (RNG_M - 1)) + 1;
}

/* returns a pseudo-random number from set 1, 2, ..., RNG_M - 1 */
long
lrand48 ()
{
  long low, high, test;

  high = xsubi / RNG_Q;
  low = xsubi % RNG_Q;
  test = RNG_A * low - RNG_R * high;
  if (test > 0)
    xsubi = test;
  else
    xsubi = test + RNG_M;
  return xsubi;
}

#ifdef TEST_RNG

main ()
{
  long i, random;

  srand48(0L);

  for (i = 1; i < 10000; i++)
    (void) lrand48 ();

  random = lrand48 ();
  printf ("z[10001] = %ld, should be 1043618065\n", random);
  if (random == 1043618065)
    printf ("success!!!\n");
}

#endif

#endif 

#ifdef Alg_AS_183

/* a very good random number generator, but slow on machines without
   floating point accelerators, this algorithm is correct on any machine
   with 16 bit integers, and at least 16 bits of precision in the
   floating point format, this algorithm is from:

   B. A. Wichmann and I. D. Hill, "An efficient and Portable Pseudo-random
	Number Generator", Applied Statistics, 1982, p 188ff, Algorithm AS 183.

   This algorithm uses three Lehmer grammers, which generate random numbers
   between 0 and 1.  It adds them together and then discards the integer
   part. This gives a better result than a single Lehmer grammar.

   I do not know the period, but it seems to be 2^45.
   Returns integers in the range 1 to 32767.
 */

#define M_0 30269
#define A_0 171
#define Q_0 177
#define R_0 2

#define M_1 30307
#define A_1 172
#define Q_1 176
#define R_1 35

#define M_2 30323
#define A_2 170
#define Q_2 178
#define R_2 63

#define LARGEST_INT 32767

static short xsubi[3];

unsigned short *
seed48 (seed16v)
unsigned short seed16v[3];
{
  static unsigned short old_xsubi[3];
  int i;

  for (i = 0; i < 3; i++)
    {
      old_xsubi[i] = xsubi[i];
      xsubi[i] = seed16v[i];
    }

  return (old_xsubi);
}

void
srand48 (seedval)
unsigned long seedval;
{
  /* set each seed to value between 1 and m-1 */

  xsubi[0] = (seedval % (M_0 - 1)) + 1;
  xsubi[1] = (seedval % (M_1 - 1)) + 1;
  xsubi[2] = (seedval % (M_2 - 1)) + 1;
}

long
lrand48 ()
{
  float random;
  long result;

  xsubi[0] = A_0 * (xsubi[0] % Q_0) - R_0 * (xsubi[0] / Q_0);
  xsubi[1] = A_1 * (xsubi[1] % Q_1) - R_1 * (xsubi[1] / Q_1);
  xsubi[2] = A_2 * (xsubi[2] % Q_2) - R_2 * (xsubi[2] / Q_2);
  
  if (xsubi[0] < 0) xsubi[0] += M_0;
  if (xsubi[1] < 0) xsubi[1] += M_1;
  if (xsubi[2] < 0) xsubi[2] += M_2;

  random = ((float)xsubi[0]/(float)M_0) + ((float)xsubi[1]/(float)M_1) +
    ((float)xsubi[2]/(float)M_2);

  /* 0 < random < 3, discard integer part, original Fortran alg used
     AMOD(random, 1.0), this should accomplish the same thing */
  if (random > 2)      random -= 2.0;
  else if (random > 1) random -= 1.0;

  /* original alg returns random, but we want an integer, so multiply
     by MAX_INT and cast into a long, note that MAX_INT must be
     representable exactly in the floating point format */
  result = (long) (random * (float)LARGEST_INT);
  
  return result;
}

#ifdef TEST_RNG

main ()
{
  long i, random;

  srand48(1L);

  for (i = 1; i < 10000; i++)
    (void) lrand48 ();

  random = lrand48 ();
  printf ("z[10001] = %ld, should be 22773\n", random);
  if (random == 22773)
    printf ("success!!!\n");
}

#endif

#endif