Bug in rand() and srand()

[Michael Mauldin]

[ generic apology: my mail reply to cubsvax!peters bounced back ]

Rand is a simple linear congruential pseudo-random number generator
(read your Knuth volume II).  The low order bits of a linear
congruential sequence are not very random at all.  The problem is not
rand, but rather your use of rand.  You cannot safely use the result
directly, and you cannot safely use mod (%) to give a smaller range.
You must use div (/) to get the range, or at the very least throw away
the low bits before 'mod'ing.

The right way to use a LCPRNG to get integers frlom 0..n-1 is this:

    result = floor (n * rand() / (MAXRAND+1));

Here is a valid call to the system rand():

    (int) result = (int) range * ((double) rand () / 2147483648.0);

Here is the salient portion of the system rand:

    return ((randx = randx * 1103515245 + 12345) & 0x7fffffff);

Here is my own random number generator which uses Knuth algorithm M and
a few other tricks to make more random numbers.  It also contains a
randint(n) which correctly generates pseudo-random integers from 0..n-1.

/*****************************************************************
 * rand: A very, very random generator, period approx 6.8064e16.
 *
 * Uses algorithm M, "Art of Computer Programming", Vol 2. 1969, D.E.Knuth.
 *
 * Two generators are used to derive the high and low parts of sequence X,
 * and another for sequence Y. These three generators were derived by me.
 *
 * Usage:  initialize by calling srand(seed), then rand() returns a random 
 *         integer from 0..2147483647. srand(0) uses the current time as
 *         the seed. You must call srand() first, because it initializes
 *	   an array which is used by the rand() function.
 *
 * Author: Michael Mauldin, June 14, 1983.
 *****************************************************************/

/* Rand 1, period length 444674 */
# define MUL1 1156
# define OFF1 312342
# define MOD1 1334025
# define RAND1 (seed1=((seed1*MUL1+OFF1)%MOD1))
# define Y      RAND1

/* Rand 2, period length 690709 */
# define MUL2 1366
# define OFF2 827291
# define MOD2 1519572
# define RAND2 (seed2=((seed2*MUL2+OFF2)%MOD2))

/* Rand 3, period length 221605 */
# define MUL3 1156
# define OFF3 198273
# define MOD3 1329657
# define RAND3 (seed3=((seed3*MUL3+OFF3)%MOD3))

/*
 * RAND2 generates 19 random bits, RAND3 generates 17. The X sequence
 * is made up off both, and thus has 31 random bits.
 */

# define X    ((RAND2<<13 ^ RAND3>>3) & 017777777777)

# define AUXLEN 97
static int seed1=872978, seed2=518652, seed3=226543, auxtab[AUXLEN];

/****************************************************************
 * srand: Initialize the seed (0 => use the current time+pid+uid)
 ****************************************************************/

srand (seed)
int seed;
{ register int i;

  if (seed == 0) seed = (getuid() + getpid() + time(0));  

  /* Set the three random number seeds */
  seed1 = (seed1+seed) % MOD1;
  seed2 = (seed2+seed) % MOD2;
  seed3 = (seed3+seed) % MOD3;
  
  for (i=AUXLEN; i--; )
    auxtab[i] = X;
}

/****************************************************************
 * rand: Random integer from 0..2147483647
 ****************************************************************/

int rand ()
{ register int j, result;

  j = AUXLEN * Y / MOD1;	/* j random from 0..AUXLEN-1 */
  result = auxtab[j];
  auxtab[j] = X;
  return (result);
}

/****************************************************************
 * randint: Return a random integer in a given Range
 ****************************************************************/

int randint (range)
int range;
{ double rrand;
  int result;

  rrand = rand();
  result = range * (rrand / 2147483648.);
  return (result);
}

Share and Enjoy!

Michael L. Mauldin (Fuzzy)		Department of Computer Science
Mauldin at CMU-CS-CAD.ARPA			Carnegie-Mellon University
(412) 578-3065				Pittsburgh, PA  15213