49.2.4 Generating integers in a specified range

In applications using a randomness source, it happens to need a pseudo–random integer X in a given range 0 <= X < U, while the PRNG generates a pseudo–random integer N in the range 0 <= N < M. In other words: We need a sample from the range 0 <= X < U having uniform probability distribution, by means of a sample from the range 0 <= N < M having a uniform probability distribution.

We distinguish the two cases U <= M and M < U, because when U <= M there are more integers in the range 0 <= N < M than in the range 0 <= X < U, while when M < U it is the other way around.

• When U <= M, we have to understand that doing a simple N mod U is incorrect, because it breaks uniform distribution of the result in the range 0 <= X < U. We can use the following formulation which gives uniformly distributed results:

  Q  = floor(M / U)
  QU = Q * U
  
  loop:
    N = <generate the next integer>
    if (N < QU)
    then X = floor(N / Q)
    else goto loop
  

  notice that QU = Q * U <= M; also notice that it is statistically possible that the algorithm loops forever, but we can hope that it finds a solution in a reasonable short time. It works like this:

  1. Conceptually partition the range 0 <= N < QU in U intervals each of length Q, indexed by X in the range 0 <= X < U.
  2. Generate an integer N and discard it if QU <= N.
  3. If N is in the range 0 <= N < QU, then it is also in a range X * Q < N < (X+1) * Q. X is the result.
• When M < U, we have to generate pseudo–random integers N' uniformly distributed in a range 0 <= N' < M', with M < U <= M'. We can do it by computing the following formula, which is a polynomial in M and a linear combination with positive coefficients in N(i) with 0 <= i < k:

  N' = N0 + M * N1 + M^2 * N2 + ... +
     + ... + M^(k-2) * N(k-2) + M^(k-1) * N(k-1)
  

  where N0, N1, ..., N(k-1) are all generated integers in the range 0 <= N(i) < M. The maximum value, that is M' - 1, is realised when every N(i) equals M - 1:

  M' - 1 = (M - 1) + M * (M - 1) + M^2 * (M - 1) + ... +
         + ... + M^(k-2) * (M - 1) + M^(k-1) * (M - 1)
         = (M - 1) * [1 + M + M^2 + ... + M^(k-2) + M^(k-1)]
         = (M - 1) * (1 - M^k)/(1 - M)
         = (M - 1) * (M^k - 1)/(M - 1)
         = M^k - 1
  

  which implies M^k = M'. So by selecting k such that U <= M^k we are sure that N' will be in a suitable range.

  Computing this polinomial/combination is like expressing the big pseudo–random integer N' in base M where N(i) are the digits. Each possible N' is uniquely associated to a k–tuple of integers N(i). Being the probability of N(i) uniformly distributed, each possible k–tuple has uniform probability in the set of all the possible k–tuples. So also the probability of N' is uniformly distributed in the range 0 <= N' < M^k.

  Then we proceed like we did before for U <= M:

  Q  = floor(M' / U)
  QU = Q * U
  
  loop:
    N' = <compute the next polynomial>
    if (N' < QU)
    then X = floor(N' / Q)
    else goto loop