Shuffle Box

January 17, 2014

In the previous exercise we implemented the minimum standard random number generator of Park and Miller. In today’s exercise we implement a similar random number generator, due to Knuth, and show how both random number generators can be improved.

The Park-Miller random number generator is a linear congruential generator of the form Xn+1 = a · Xn (mod m). A disadvantage of this type of random number generator is that it cannot produce a value of zero, since then all future random numbers will also be zero. Although we didn’t mention it, that caused some problems in the generators of the prior exercise, requiring careful code and proper parameters to prevent.

Knuth’s random number generator is a mixed linear congruential generator of the form Xn+1 = a · Xn + c (mod m). Because a constant is added after each multiplication, zero becomes a valid random number. Knuth recommends a = 69069, c = 1234567 and m = 232 for 32-bit arithmetic and a = 6364136223846793005, c = 1442695040888963407 and m = 264 for 64-bit arithmetic, though some people set c = 1.

In addition to regularizing zero, Knuth’s mixed linear congruential random number generator has the advantage of being very fast, at least in languages with fixed-length integer datatypes, because the modulo operation is never performed; the arithmetic just wraps around the fixed length of the integer. Thus, assuming the seed is an unsigned long integer, a C version of Knuth’s generator is just

unsigned long integer knuth()
{
    return ( seed = 69069 * seed + 1234567 );
}

Very simple. Very fast.

One disadvantage of the Knuth random number generator is that the least-significant bit cycles with period 2, the second-least-significant bit cycles with period 4, and so on, because the modulus is a power of two. And both the Park-Miller and Knuth generators suffer from the problem of serial correlation. In many cases this doesn’t matter. But a simple technique called a shuffle box improves the randomness of both the Park-Miller and Knuth generators.

A shuffle box consists of an array of the k most recent seeds, initialized by the first k random numbers following the initial seed. Each time a new random number is requested, an element of the array j is selected by j = ⌊ k · seed / m ⌋, the jth element of the array is returned as the random number, and the next number in the random sequence is stored in the jth element of the array. K is often a power of two, because then j can be quickly computed by a right shift, though it isn’t necessary. This process is described by Knuth in Algorithm 3.2.2B.

Your task is to write a basic version of Knuth’s random number generator, then write versions of both the Park-Miller and Knuth random number generators that use shuffle boxes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

About these ads

Pages: 1 2

One Response to “Shuffle Box”

  1. Paul said

    In Python. I.s.o. creating dedicated shuffle boxes for Park Miller and Knuth prng, I created a shuffle box, that takes any prng.

    from __future__ import division
    import time
    T0 = time. clock()
    
    def get_seed(seed=None):
        if seed is None:
            return (int(time.clock() * 2147483647 * 2147483647))
        return seed
        
    def park_miller(seed=None):
        """  yields floats in range (0, 1)
        """
        M = 2147483647
        A = 16807
        seed = get_seed(seed)
        while 1:
            seed = (A * seed) % M
            yield seed / M
            
    def knuth(seed=None, bits=32):
        """   yields floats in range [0, 1)"""
        a, c, m = ((69069, 1234567, 2 ** 32) if bits == 32 else
            (6364136223846793005, 1442695040888963407, 2 ** 64))
        seed  = get_seed(seed)
        while 1:
            seed = (seed * a + c) % m
            yield seed / m
    
    def prng_shuffle(prng, nboxes=55):
        """generator that yields a uniform random number in range[0,1)
            prng: function that returns a uniform random number in range[0,1)
        """
        boxes = [prng() for _ in range(nboxes)]
        index = int(prng() * nboxes)
        while 1:
            result = boxes[index]
            boxes[index] = prng()
            index = int(result * nboxes)
            yield result
            
    #example
    rand = prng_shuffle(knuth(bits=32).next).next
    print rand()
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 600 other followers

%d bloggers like this: