Generating Random Sequences (part I)

Every once in a while, we need a random sequence. Whether to test a data structure’s performance or to run probabilistic unit tests, the provided rand primitive from your favorite programming language has several limitations. First, it’s been known for a while that if most implementations of the C standard library rand() function are not very random, despite being “good enough” in great many cases. Second, and more importantly, it does not allow you to easily control the period nor to generate a permutation on , for example.

There are many methods of generating pseudo-random number sequences. Not all exhibit the same properties and, accordingly, a method may be more useful in one case and perfectly useless in another. High quality pseudo-random number generation is a notoriously hard endeavor, but there is a number of very simple algorithms that will get you out of trouble for certain specific tasks. Let us consider, for example, the example where the pseudo-random generator must generate the numbers in exactly once, in a random order, of course, in exactly draws.

Let’s say, then, that for some reason, you have to generate a random-looking permutation of the numbers . The first approach one may think of is to generate an array with the numbers in it and randomize its contents. The array is then scanned from left to right and all the numbers are drawn one by one.

This approach, although simple, suffers from two important drawbacks. The first is that you must allocate the space for numbers and initialize the array. This is all cute when , but a lot less when or more: you may not be able to fit the list in memory at all. The second is that after the allocation of the array, you still need to shuffle the list quite well so that it looks random enough. This, in turn, asks for a good pseudo-random number generator and an efficient shuffling algorithm. Even the Fisher-Yates Shuffle asks for operations. A C99 implementation would look like:

void swap_int(int * a, int * b)
 {
  int t=*a;
  *a=*b;
  *b=t;
 }

bool fisher_yates(int n, int ** array)
 {
  (*array)=(int *)malloc(n);

  if (*array)
   {
    // let us fill the array
    for (int i=0;i<n;i++) (*array)[i]=i;

    // fisher-yates loop
    for (int i=n-1;i;i--)
     {
      int j=rand() % i; // 0...i-1
      swap_int( *array+i, *array+j);
     }
    return true;
   }
  return false;
 }

This seems nice, but 1) calls on rand() quite a lot and 2) must allocate a size array of integers. Shrinking the integer type to uint16_t if is small enough is not much of an optimization. What we need, is a method that generates a permutation without holding such a large internal state. Better yet, we need something that has storage state in or less.

Turns out that we can do that, provided that we add a bit of initialization time—and not much of it. We chose , the size of the array—which should now be called the length of the sequence. We chose a number that is relatively prime to . That is, and share no common factors. Having just prime also works, but that’s not necessary to use that strong a condition; relatively prime suffice quite enough. Third, choose , the seed, to be in . Maybe initialise it with time() % n or some similar method¹.

The sequence is generated by the following iteration:

So this algorithm will generate a permutation of using only and , bothintegers being bits long.

The C99 implementation of this algorithm would look like:

typedef
struct
 {
  int x,   // the current state
      n,   // the size of the array
      rho; // and the relatively prime number
 } generator_t;

int get_next( generator_t * gen )
 {
  return gen->x=(gen->x+gen->rho) % gen->n;
 }

and it could be called:

  generator_t gen={4,20,13}; // seed, size, rho

  for (int i=0;i<20;i++)
   printf(" %d",get_next(&gen));
  printf("\n");

producing the most cromulent output:

17 10 3 16 9 2 15 8 1 14 7 0 13 6 19 12 5 18 11 4

*
* *

So, how does this generator works exactly? The fact that is relatively prime to plays a key role in the generator. The generator has been presented as a recurrence but can be rewritten as , where is the seeded value.

Let us show (with rather much handwaving) that the generator has a period of and that it covers all the values . Let us consider the seed to be zero for now; the generator becomes . Since is relatively prime to , the values of that solve are multiples of (including 0). Therefore, you’ll have to iterate times to fall back on 0; values all give non-zero residuals. Thus demonstrating that the period of the generator is indeed .

Now, let us show that the values generated by varying from 0 to are all distinct. Suppose we have and , . We can show that if , , exactly because . Indeed: if for some , then , which implies that , so that .

If you’re using a non zero seed, you can redo the above adding in the proof; the net resulting being of “rotating” the sequence by an offset given by the seed .

*
* *

The choice of is a difficult one. First, if , the sequence generated is merely a count (if ) from the seed wrapping at and finishing at the seed minus one, or a decount (if ), which isn’t really random. Same if , except that we’re (de)counting in step of two. For the sequence to look random enough, must be such that , where is Phidias‘ constant, the golden ratio. This celebrated constant is given by:

Can you see why? I’ll let you think this one out for a while² :)

---

¹ Personally, I favor the Linux-specific function gettimeofday that returns the time in seconds and in microseconds. Using the complete time with microsecond gives a rather good time-dependant seed.

² Or maybe does quite the job?