Random Number Generators

Since random number generators are useful, and are used in many assignments, we will highlight a few aspects of random number generation, starting with the most basic generator: (source file )


public class uniform_random {
  // Basic Lehmer generator - constants
  static final long m = 2147483647L;
  static final long a = 48271L;
  static final long q = 44488L;
  static final long r = 3399L;

  // "Global" variable - the seed, set to some 
  // arbitrary non-zero value.
  static long r_seed = 12345678L; 

  // Basic Lehmer generator - uniform[0,1]
  // For more information see Knuth, Vol. II.
  public static double uniform ()
  {
    long hi = r_seed / q;
    long lo = r_seed - q * hi;
    long t = a * lo - r * hi;
    if (t > 0)
      r_seed = t;
    else
      r_seed = t + m;
    return ( (double) r_seed / (double) m );
  }

  public static void main (String[] argv)
  {
    // Let's test the generator. The mean should be 0.5.
    double sum = 0;
    for (int i=1; i<=10000; i++)
      sum += uniform ();

    System.out.println ("Average of 10000 samples: " + sum/10000);
  }

}

This generator is an example of the well-known Lehmer class of random number generators:

It has been known that certain arithmetic operations on prime numbers, when performed repeatedly, result in almost-random digits being produced. The mod operator is an example.
We will not go into the theory of random generation here. Instead, let's play with the code to understand how to use it.

Consider the following modification to the above code: (source file )


public class uniform_random {
  // Basic Lehmer generator - constants
  static final long m = 2147483647L;
  static final long a = 48271L;
  static final long q = 44488L;
  static final long r = 3399L;

  // "Global" variable - the seed, set to some 
  // arbitrary non-zero value.
  static long r_seed = 12345678L; 

  // Basic Lehmer generator - uniform[0,1]
  // For more information see Knuth, Vol. II.
  public static double uniform ()
  {
    long hi = r_seed / q;
    long lo = r_seed - q * hi;
    long t = a * lo - r * hi;
    if (t > 0)
      r_seed = t;
    else
      r_seed = t + m;

    double returnValue = ( (double) r_seed / (double) m );
    System.out.println ("r_seed=" + r_seed + " m=" + m + " returnValue=" + returnValue);   
    return returnValue;
  }

  public static void main (String[] argv)
  {
    // Print out 10 values.
    for (int i=1; i<=10; i++)
      uniform ();
  }

}

Here's what gets printed out:

r_seed=1085252519 m=2147483647 returnValue=0.5053600852868334
r_seed=508259731 m=2147483647 returnValue=0.23667688073435653
r_seed=1352291773 m=2147483647 returnValue=0.6297099281240767
r_seed=1563240271 m=2147483647 returnValue=0.7279404773041329
r_seed=890733155 m=2147483647 returnValue=0.414779947798131
r_seed=1810028418 m=2147483647 returnValue=0.8428601635819581
r_seed=1509587083 m=2147483647 returnValue=0.7029562647002545
r_seed=862973489 m=2147483647 returnValue=0.4018533459873187
r_seed=1852986660 m=2147483647 returnValue=0.8628641538614706
r_seed=677683663 m=2147483647 returnValue=0.3155710470469534

Notice that the r_seed value is always less than m and therefore r_seed / m is always less than 1.
Also, the operations ensure that r_seed is always positive. Therefore, we always return a value between 0 and 1.
It is the stream of numbers so produced that is random.

Now, it's important to realize that the random numbers so produced have the property of being uniform:

It so happens that the r_seed values produced are always between 1 and m.
There are of course a finite number of values between 1 and m.
After a (long) while, the pattern repeats.
(In fact, after m calls to uniform()).
Between repetitions of the entire cycle, no number is repeated.
Thus, no number is particularly "favored" over any other.
This has an important consequence: the numbers generated between 0 and 1 are uniformly distributed.

Next, an obvious question is: what should the initial value of r_seed be?
Answer:

It can be any positive number between 1 and m.
In the above code, it's arbitrarily chosen to be the long constant 12345678L. It could just as well have been initialized to 5.

If you execute the file once, you get the output:

r_seed=1085252519 m=2147483647 returnValue=0.5053600852868334
r_seed=508259731 m=2147483647 returnValue=0.23667688073435653
r_seed=1352291773 m=2147483647 returnValue=0.6297099281240767
r_seed=1563240271 m=2147483647 returnValue=0.7279404773041329
r_seed=890733155 m=2147483647 returnValue=0.414779947798131
r_seed=1810028418 m=2147483647 returnValue=0.8428601635819581
r_seed=1509587083 m=2147483647 returnValue=0.7029562647002545
r_seed=862973489 m=2147483647 returnValue=0.4018533459873187
r_seed=1852986660 m=2147483647 returnValue=0.8628641538614706
r_seed=677683663 m=2147483647 returnValue=0.3155710470469534

Now if you execute again, you get:

r_seed=1085252519 m=2147483647 returnValue=0.5053600852868334
r_seed=508259731 m=2147483647 returnValue=0.23667688073435653
r_seed=1352291773 m=2147483647 returnValue=0.6297099281240767
r_seed=1563240271 m=2147483647 returnValue=0.7279404773041329
r_seed=890733155 m=2147483647 returnValue=0.414779947798131
r_seed=1810028418 m=2147483647 returnValue=0.8428601635819581
r_seed=1509587083 m=2147483647 returnValue=0.7029562647002545
r_seed=862973489 m=2147483647 returnValue=0.4018533459873187
r_seed=1852986660 m=2147483647 returnValue=0.8628641538614706
r_seed=677683663 m=2147483647 returnValue=0.3155710470469534

It's the same output! So, how can we call this random?

Here's what's going on:

Simulations are never written to always start differently.
The idea is, if there's a bug, you should be able to reproduce the "randomness" exactly, to help in debugging.
To actually use multiple random values for experimental purposes, you can do one of the following:
- Edit the file to use a different initial r_seed.
  This is cumbersome if you want thousands of "runs".
- The better option is to create a loop to repeat your experiments, and never reset the seed in the loop.

For example, here's how we create an array of 5 random values repeatedly:

public class random_example {

  static final long m = 2147483647L;
  static final long a = 48271L;
  static final long q = 44488L;
  static final long r = 3399L;

  static long r_seed = 12345678L; 

  public static double uniform ()
  {
    long hi = r_seed / q;
    long lo = r_seed - q * hi;
    long t = a * lo - r * hi;
    if (t > 0)
      r_seed = t;
    else
      r_seed = t + m;
    return ( (double) r_seed / (double) m );
  }

  public static void main (String[] argv)
  {
    double[] A = new double [5];
    int howManyArrays = 3;
    for (int n=1; n<=howManyArrays; n++) {
      // Generate random values; 
      for (int i=0; i<5; i++) 
        A[i] = uniform();

      // Print out array. 
      System.out.print ("Array #" + n + " : ");
      for (int i=0; i<5; i++) 
        System.out.print ("  " + A[i]);
      System.out.println ();
    }
    
  }

}

Thus, to perform estimation, we will need to run an experiment multiple times and take an average:

As an example, consider the following experiment:

Suppose we generate an array of 5 random values (between 0 and 1).
What is the probability that the values generated in increasing order?
To estimate this probability, let us generate 10,000 such arrays and count how many of them have this property.

Here is the code: (source file )

public class random_example2 {

  static final long m = 2147483647L;
  static final long a = 48271L;
  static final long q = 44488L;
  static final long r = 3399L;

  static long r_seed = 12345678L; 

  public static double uniform ()
  {
    long hi = r_seed / q;
    long lo = r_seed - q * hi;
    long t = a * lo - r * hi;
    if (t > 0)
      r_seed = t;
    else
      r_seed = t + m;
    return ( (double) r_seed / (double) m );
  }

  public static void main (String[] argv)
  {
    double[] A = new double [5];
    int numArrays = 10000;
    int numOccurences = 0;

    for (int n=1; n<=numArrays; n++) {
      for (int i=0; i<5; i++) 
        A[i] = uniform();

      // Check for property.
      boolean property = true;
      for (int i=1; i<5; i++) {
        if (A[i] < A[i-1])
          property = false;
      }

      // Count.
      if (property)
        numOccurences++;
    }
    
    // Now estimate probability.
    double prob = (double) numOccurences / (double) numArrays;
    System.out.println ("Probability estimate: " + prob);
  }

}