Random Numbers in a Range Using Generic Programming

Consider the generation of random numbers in an interval, say the generation of an integer in {0,1,...,m-1} where the probability of each integer is 1/m. Usually, you have a source of randomness that provides random bits in chunks of architecture word size—something like the wrapper around the mtrand library (bedaux.net/mtrand); see Listing One.

Mtrand implements the Mersenne Twister random-number generator (www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html), which is a high-quality source of randomness. By "high quality," I mean that you don't need to consider low-order bits of words produced by the generator to be of inferior quality than high-order bits. In principle, you can take the value returned by next() modulo m to produce an integer in the range [0,m-1].

However, this approach is problematic when m is comparable in size to the number of values representable by value_t. For instance, assume that the type of value_t is an unsigned 10-bit integer, and can thus contain values in the [0,1023] range. Consider m=684, which is approximately 2/3 of 1024. Taking a random value_t modulo m results in the lower half of [0,m-1] range having twice the probability of its upper half.

The regular approach to this problem is to ignore it (effectively assuming that m is low enough that the difference in probabilities is insignificant), or to reject the values in the "remainder" portion of the value_t range; see Listing Two.

Taking the previous example with m=684, the constants in Listing Two are:

• M=1023, the maximal value representable by value_t.
• rem=340, the remainder size.
• last=683, the last value in the range of "safe" values returned by next().

Is the remainder values rejection approach costly? It is, especially for values of m close to, but above, the half-range of value_t. For such values, the expected number of calls to next() for each call to next(max) is 2 (this follows from the probability of rejection being 1/2). Because the cost of random bits varies between "moderate" (for Mersenne Twister) to "expensive" (for cryptography-oriented generators, especially hardware implementations), the running time of the relevant program portions doubles in such cases (unless you settle for the low-quality linear congruential generators).

Is it wasteful? The rejected values are thrown away, but perhaps their inherent entropy can still be used—maybe to reduce the probability of rejection on the next iteration of the do/while loop? For instance, in the example, once rnd is rejected, and just before the next loop iteration, rnd-last-1 is uniformly distributed in the [0,rem-1] range. You could partition the [0,m-1] range into exact rem buckets, and in the next iteration draw a value in the corresponding bucket using m=rem. What's problematic in this solution is that the conditions for its application are quite rare because it requires that rem divides m. It can happen—with a value_t type of 10 bits, two possibilities are m=768 and m=320—however, it is possible to do a little better using the Greatest Common Divisor (GCD) algorithm.

Recall that the GCD algorithm computes the maximal number that divides its two arguments, in logarithmic number of arithmetic operations. Let's return to m=684 and rem=340. Both are divisible by 4 (which is the GCD). Thus, if rnd is rejected in the first loop iteration, you can find out in which quarter of the remainder values rnd is located, map it to the corresponding quarter of values in the [0,m-1] range, and let m=684/4=171 in the next iteration. Then, the probability of rejection becomes 1/6 instead of 1/3 (and can in general drop to a much lower value).

It turns out, however, that for randomness sources that are only moderately expensive, the GCD computation cancels any running time benefit attained by reducing the expected number of calls to next(). If you are writing this code for a library, you could specify that it assumes usage with only "expensive" random-number generators, or implement memoization. However, in the first case, the code's generic nature suffers, and in the second case, a penalty is incurred on users who don't need perfect distribution of values in the [0,m-1] range.

This gloomy situation does have a win-win solution if you are using a language that allows compile-time computations (C++, for instance), and if you are willing to limit the approach outlined earlier to values of m known at compile time.

class Random
{
public:
   typedef unsigned value_t;
   Random(value_t seed) : mtrand(seed) { }
   virtual ~Random() { }
   value_t next() { return mtrand(); }
private:
   MTRand_int32 mtrand;
};

class SimpleRange : public Random
{
public:
   SimpleRange(value_t seed) : Random(seed) { }
   using Random::next;
   value_t next(value_t max);
};
// Implementation
SimpleRange::value_t SimpleRange::next(value_t max)
{
   const value_t M    = std::numeric_limits<value_t>::max();
   const value_t rem  = (M - (max - 1)) % max;
   const value_t last = M - rem;
   value_t rnd;
   // Once rnd is in the safe range of
   // [0 .. last], return rnd % max
   do
      // rnd <- [0 .. M]
      rnd = next();
   while (rnd > last);
   return rnd % max;
}