Fast Random Integer Generation in an Interval

05/28/2018
by   Daniel Lemire, et al.
0

In simulations, probabilistic algorithms and statistical tests, we often generate random integers in an interval (e.g., [0,s)). For example, random integers in an interval are essential to the Fisher-Yates random shuffle. Consequently, popular languages like Java, Python, C++, Swift and Go include ranged random integer generation functions as part of their runtime libraries. Pseudo-random values are usually generated in words of a fixed number of bits (e.g., 32 bits, 64 bits) using algorithms such as a linear congruential generator. We need functions to convert such random words to random integers in an interval ([0,s)) without introducing statistical biases. The standard functions in programming languages such as Java involve integer divisions. Unfortunately, division instructions are relatively expensive. We review an unbiased function to generate ranged integers from a source of random words that avoids integer divisions with high probability. To establish the practical usefulness of the approach, we show that this algorithm can multiply the speed of unbiased random shuffling on x64 processors. Our proposed approach has been adopted by the Go language for its implementation of the shuffle function.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

03/08/2020

The Fast Loaded Dice Roller: A Near-Optimal Exact Sampler for Discrete Probability Distributions

This paper introduces a new algorithm for the fundamental problem of gen...
11/09/2021

An Empirical Study of Automated Unit Test Generation for Python

Various mature automated test generation tools exist for statically type...
04/19/2020

A practical approach to testing random number generators in computer algebra systems

This paper has a practical aim. For a long time, implementations of pseu...
05/18/2018

Towards Taming Java Wildcards and Extending Java with Interval Types

Of the complex features of generic nominally-typed OO type systems, wild...
11/24/2019

Probabilistic methods of bypassing the maze using stones and a random number sensor

In this paper, some open questions that are posed in Ajans' dissertation...
04/22/2019

Interval Algorithm for Random Number Generation: Information Spectrum Approach

The problem of exactly generating a general random process (target proce...
10/25/2018

Random Sampling: Practice Makes Imperfect

The pseudo-random number generators (PRNGs), sampling algorithms, and al...

Code Repositories

fastrange

A fast alternative to the modulo reduction


view repo

fastrand

Fast random number generation in Python


view repo

FastShuffleExperiments

How fast can we shuffle values?


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

There are many efficient techniques to generate high-quality pseudo-random numbers such as Mersenne Twister (Matsumoto:1998:MTE:272991.272995), Xorshift (marsaglia2003xorshift; Panneton:2005:XRN:1113316.1113319), linear congruential generators (l1999tables; LEcuyer:1993:SGM:169702.169698; de1988parallelization; fishman2013monte; LEcuyer:1990:RNS:84537.84555) and so forth (l2017random; l2012random). Many pseudo-random number generators produce 32-bit or 64-bit words that can be interpreted as integers in and respectively: the produced values are practically indistinguishable from truly random numbers in or . In particular, no single value is more likely than any other.

However, we often need random integers selected uniformly from an interval , and this interval may change dynamically. It is useful for selecting an element at random in an array containing  elements, but there are less trivial uses. For example, the Fisher-Yates random shuffle described by Knuth (Knuth1969; Durstenfeld:1964:ARP:364520.364540) (see Algorithm 1) requires one random integer in an interval for each value in an array to be shuffled. Ideally, we would want these values to be generated without bias so that all integer values in are equally likely. Only then are all permutations equally likely. A related algorithm is reservoir sampling (Vitter:1985:RSR:3147.3165) (see Algorithm 2) which randomly selects a subset of values from a possibly very large array, even when the size of the array is initially unknown.

We use random permutations as part of simulation algorithms (Devroye:1997:RVG:268403.268413; Calvin:1998:UPR:280265.280273; Owen:1998:LSS:272991.273010; Osogami:2009:FPB:1540530.1540533; Amrein:2011:VIS:1899396.1899401; Hernandez:2012:CNO:2379810.2379813)

. The performance of randomized permutation algorithms is important. Various non-parametric tests in statistics and machine learning repeatedly permute randomly the original data. In some contexts, Hinrichs et al. found that the computational burden due to random permutations can be prohibitive 

(Hinrichs:2013:SUP:2999611.2999711). Unsurprisingly, much work has been done on parallelizing permutations (Shterev2010; Langr:2014:A9P:2684421.2669372; SANDERS1998305; Gustedt2008; Waechter2012) for greater speed.

0:  array made of elements indexed from to
1:  for  do
2:      random integer in
3:     exchange and
4:  end for
ALGORITHM 1 Fisher-Yates random shuffle: it shuffles an array of size so that  possible permutations are equiprobable.
0:  array made of elements indexed from to
0:  integer ()
1:   array of size
2:  for  do
3:     
4:  end for
5:  for  do
6:      random integer in
7:     if  then
8:        
9:     end if
10:  end for
11:  return
ALGORITHM 2 Reservoir sampling: returns an array containing  distinct elements picked randomly from an array of size so that all  possible samples are equiprobable.

One might think that going from fixed-bit pseudo-random numbers (e.g., 32-bit integers) to pseudo-random numbers in an interval is a minor, inexpensive operation. However, this may not be true, at least when the interval is a changing parameter and we desire a uniformly-distributed result.

  • Let us consider a common but biased approach to the problem of converting numbers from a large interval to numbers in a subinterval (): the modulo reduction . On x64 processors, this could be implemented through the division (div) instruction when both and are parameters. When applied to 32-bit registers, this instruction has a latency of 26 cycles (fog2016instruction). With 64-bit registers, the latency ranges from 35 to 88 cycles, with longer running times for small values of .

  • Another biased but common approach consists in using a fixed-point floating-point representation consisting of the following step:

    • we convert the random word to a floating-point number in the interval ,

    • we convert the integer into a floating-point number,

    • we multiply the two resulting floating-point numbers,

    • and we convert the floating-point result to an integer in .

    When using the typical floating-point standard (IEEE 754), we can at best represent all values in divided by using a 32-bit floating-point number. Thus we do not get the full 32-bit range: we cannot generate all numbers in if . To do so, we must use double precision floating-point numbers, and then we can represent all values in divided by . Moreover converting between floating-point values and integers is not without cost: the corresponding instructions on x64 processors (e.g., cvtsi2ss, cvttss2si) have at least six cycles of latency on Skylake processors (fog2016instruction).

While generating a whole new 64-bit pseudo-random number can take as little as a handful of cycles (Saito2008), transforming it into an integer in an interval ( for ) without bias can take an order of magnitude longer when using division operations.

There is a fast technique that avoids division and does not require floating-point numbers. Indeed, given an integer in the interval , we have that the integer is in for any integer . If the integer is picked at random in , then the result is a random integer in  (overton2011). The division by a power of two () can be implemented by a bit shift instruction, which is inexpensive. A multiplication followed by a shift is much more economical on current processors than a division, as it can be completed in only a handful of cycles. It introduces a bias, but we can correct for it efficiently using the rejection method (see § 4). This multiply-and-shift approach is similar in spirit to the multiplication by a floating-point number in the unit interval () in the sense that can be intuitively compared with where is a random number in .

Though the idea that we can avoid divisions when generating numbers in an interval is not novel, we find that many standard libraries (Java, Go, …) use an approach that incurs at least one integer division per function call. We believe that a better default would be an algorithm that avoids division instructions with high probability. We show experimentally that such an algorithm can provide superior performance.

2. Mathematical Notation

We let be the largest integer smaller than or equal to , we let be the smallest integer greater than or equal to . We let be the integer division of by , defined as . We define the remainder of the division of by as : .

We are interested in the integers in an interval where, typically, or . We refer to these integers as -bit integers.

When we consider the values of as goes from to , we get the values

We have the following lemma by inspection.

Lemma 2.1 ().

Given integers , there are exactly  multiples of in whenever divides . More generally, there are exactly  integers in having a given remainder with whenever divides .

A geometric distribution with success probability

is a discrete distribution taking value with probability . The mean of a geometric distribution is .

3. Existing Unbiased Techniques Found in Common Software Libraries

Assume that we have a source of uniformly-distributed -bit random numbers, i.e., integers in . From such a source of random numbers, we want to produce a uniformly-distributed random integer in for some integer . That is all integers from the interval are equally likely: for any integer . We then say that the result in unbiased.

If divides , i.e., it is a power of two, then we can divide the random integer from by . However, we are interested in the general case where may be any integer value in .

We can achieve an unbiased result by the rejection method (von1961various). For example, we could generate random integers until is in , rejecting all other cases.111For some applications where streams of random numbers need to be synchronized, the rejection method is not applicable (asmussen2007stochastic; bratley2011guide; law2007; fu2017history; l2015random). Rejecting so many values is wasteful. Popular software libraries use more efficient algorithms. We provide code samples in Appendix A.

3.1. The OpenBSD Algorithm

The C standard library in OpenBSD and macOS have an arc4random_uniform function to generate unbiased random integers in an interval . See Algorithm 3. The Go language (e.g., version 1.9) has adopted the same algorithm for its Int63n and Int31n functions, with minor implementation differences (gorand). The GNU C++ standard library (e.g., version 7.2) also relies on the same algorithm (gnucpplib).

The interval has size which is divisible by . From Lemma 2.1, if we generate random integers from integer values in as remainders of a division by (), then each of the integers occur for  integers . To produce integers in , we use the rejection method: we generate integers in but reject the result whenever . If we have a source of unbiased random integers in , then the result of Algorithm 3 is an unbiased random integer in .

The number of random integers consumed by this algorithm follows a geometric distribution, with a success probability . On average, we need  random words. This average is less than two, irrespective of the value of . The algorithm always requires the computation of two remainders.

0:  source of uniformly-distributed random integers in
0:  target interval with
1:  {}
2:   random integer in
3:  while  do {Application of the rejection method}
4:      random integer in
5:  end while
6:  return
ALGORITHM 3 The OpenBSD algorithm.

There is a possible trivial variation on the algorithm where instead of rejecting the integer from when it is part of the first  values (in ), we reject the integer from when it is part of the last  values (in ).

3.2. The Java Approach

It is unfortunate that Algorithm 3 always requires the computation of two remainders, especially because we anticipate such computations to have high latency. The first remainder is used to determine whether a rejection is necessary (), and the second remainder is used to generate the value in as .

The Java language, in its Random class, uses an approach that often requires a single remainder (e.g., as of OpenJDK 9). Suppose we pick a number and we compute its remainder . Having both and , we can determine whether is allowable (is in ) without using another division. When it is the case, we can then return without any additional computation. See Algorithm 4.

0:  source of uniformly-distributed random integers in
0:  target interval with
1:   random integer in
2:   {}
3:  while  do {Application of the rejection method}
4:      random integer in
5:     
6:  end while
7:  return
ALGORITHM 4 The Java algorithm.

The number of random words and remainders used by the Java algorithm follows a geometric distribution, with a success probability . Thus, on average, we need  random words and remainders. Thus when is small (), we need  random words and remainders. This compares favorably to the OpenBSD algorithm that always requires the computation of two remainders. However, the maximal number of divisions required by the OpenBSD algorithm is two, whereas the Java approach could require infinitely many divisions in the worst case.

4. Avoiding Division

Though arbitrary integer divisions are relatively expensive on common processors, bit shifts are less expensive, often requiring just one cycle. When working with unsigned integers, a bit shift is equivalent to a division by a power of two. Thus we can compute quickly for any power of two . Similarly, we can compute the remainder of the division by a power of two as a simple bitwise logical AND: .

Moreover, common general-purpose processors (e.g., x64 or ARM processors) can efficiently compute the full result of a multiplication. That is, when multiplying two 32-bit integers, we can get the full 64-bit result and, in particular, the most significant 32 bits. Similarly, we can multiply two 64-bit integers and get the full 128-bit result or just the most significant 64 bits when using a 64-bit processor. Most modern programming languages (C, C++, Java, Swift, Go…) have native support for 64-bit integers. Thus it is efficient to get the most significant 32 bits of the product of two 32-bit integers. To get the most significant 64 bits of the product of two 64-bit integers in C and C++, we can use either intrinsics (e.g., __umulh when using the Visual Studio compiler) or the __uint128_t extension supported by the GNU GCC and LLVM’s clang compilers. The Swift language has the multipliedFullWidth function that works with both 32-bit and 64-bit integers. It gets compiled to efficient binary code.

Given an integer , we have that . By multiplying by , we take integer values in the range and map them to multiples of in . By dividing by , we map all multiples of in to 0, all multiples of in to one, and so forth. The  interval is . By Lemma 2.1, there are exactly multiples of in intervals since divides the size of the interval (). Thus if we reject the multiples of that appear in , we get that all intervals have exactly multiples of . We can formalize this result as the next lemma.

Lemma 4.1 ().

Given any integer , we have that for any integer , there are exactly values such that and .

Algorithm 5 is a direct application of Lemma 4.1. It generates unbiased random integers in for any integer .

0:  source of uniformly-distributed random integers in
0:  target interval with
1:   random integer in
2:  
3:  
4:  if  then {Application of the rejection method}
5:     {}
6:     while  do
7:         random integer in
8:        
9:        
10:     end while
11:  end if
12:  return
ALGORITHM 5 An efficient algorithm to generate unbiased random integers in an interval.

This algorithm has the same probabilistic distribution of random words as the previously presented algorithms, requiring the same number of -bit uniformly-distributed random words on average. However, the number of divisions (or remainders) is more advantageous. Excluding divisions by a power of two, we have a probability of using no division at all. Otherwise, if the initial value of is less than , then we incur automatically the cost of a single division (to compute ), but no further division (not counting divisions by a power of two). Table 1 and Fig. 1 compare the three algorithms in terms of the number of remainder computations needed.

expected number of remainders per integer in expected number of remainders when the interval is small () maximal number of remainders per integer in
OpenBSD (Algorithm 3) 2 2 2
Java (Algorithm 4) 1
Our approach (Algorithm 5) 0 1
Table 1. Number of remainder computations (excluding those by known powers of two) in the production an unbiased random integer in .
Figure 1. Expected number of remainder computations (excluding those by known powers of two) in the production an unbiased random integer in for 32-bit integers.

5. Experiments

We implemented our software in C++ on a Linux server with an Intel (Skylake) i7-6700 processor running at . This processor has of L1 data cache, of L2 cache per core with of L3 cache, and of RAM (DDR4 2133, double-channel). We use the GNU GCC 5.4 compilers with the “-O3 -march=native” flags. To ensure reproducibility, we make our software freely available.222https://github.com/lemire/FastShuffleExperiments Though we implement and benchmark a Java-like approach (Algorithm 4), all our experiments are conducted using C++ code.

For our experiments, we use a convenient and fast linear congruential generator with the recurrence formula where to update the 128-bit state of the generator ( for ), returning as a 64-bit random word (l1999tables). We start from a 128-bit seed . This well-established generator passes difficult statistical tests such as Big Crush (LEcuyer:2007:TCL:1268776.1268777). It is well suited to x64 processors because they have fast 64-bit multipliers.

We benchmark the time required, per element, to randomly shuffle arrays of integers having different sizes. We can consider array indexes to be either 32-bit or 64-bit values. When working with 64-bit indexes, we require 64-bit integer divisions which are slower than 32-bit integer divisions on x64 processors. We always use the same 64-bit random-number generator, but in the 32-bit case, we only use the least significant 32 bits of each random word. For reference, in both the 32-bit and 64-bit figures, we include the results obtained with the shuffle functions of the standard library (std::shuffle), implemented using our 64-bit random-number generator. For small arrays, the std::shuffle has performance similar to our implementation of the OpenBSD algorithm, but it becomes slightly slower when shuffling larger arrays.

We present our experimental results in Fig. 5. We report the wall-clock time averaged over at least 5 shuffles. When the arrays fit in the cache, we expect them to remain in the cache. The time required is normalized by the number of elements in the array. As long as the arrays fit in the CPU cache, the array size does not affect the performance. As the arrays grow larger, the latency of memory access becomes a factor and the performance decreases.

  • In the 32-bit case, the approach with few divisions can be almost twice as fast as the Java-like approach which itself can be at least 50% faster than the OpenBSD-like approach.

  • When shuffling with 64-bit indexes as opposed to 32-bit indexes, our implementations of the OpenBSD-like and Java-like algorithms become significantly slower (up to three times) due to the higher cost of the 64-bit division. Thus our approach can be more than three times faster than the Java version in the 64-bit case.

The relative speed differences between the different algorithms become less significant when the arrays grow larger. In Fig. 2, we present the ratio of the OpenBSD-like approach with our approach. We see that the relative benefit of our approach diminishes when the array size increases. In the 32-bit case, for very large arrays, our approach is merely 50% faster whereas it is nearly three times faster for small arrays. Using Linux perf

, we estimated the number of cache misses to shuffle an array containing 100 million integers and found that the OpenBSD approach generates about 50% more cache misses than our approach.

Figure 2. Ratio of the timings of the OpenBSD-like approach and of our approach.
(a) 32-bit indexes
(b) 64-bit indexes
Figure 5. Wall-clock time in nanoseconds per element to shuffle arrays of random integers.

To help the processor prefetch memory and reduce the number of cache misses, we can compute the random integers in small blocks, and then shuffle while reading the precomputed integers (see Algorithm 6). The resulting buffered algorithm is equivalent to the conventional Fisher-Yates random shuffle, and it involves the computation of the same number of random indexes, but it differs on how the memory accesses are scheduled. In Fig. 8, we see that the OpenBSD-like approach benefits from the buffering when shuffling large arrays. A significant fraction of the running time of the regular OpenBSD-like implementation is due to caching issues. When applied to our approach, the benefits of the buffering are small, and for small to medium arrays, the buffering is slightly harmful.

0:  array made of elements indexed from to
0:   a small positive constant (the buffer size)
1:  
2:   some array with capacity (the buffer)
3:  while  do
4:     for  do
5:        random integer in
6:     end for
7:     for  do
8:        exchange and
9:     end for
10:     
11:  end while
12:  while  do
13:      random integer in
14:     exchange and
15:     
16:  end while
ALGORITHM 6 Buffered version of the Fisher-Yates random shuffle.
(a) 32-bit indexes
(b) 64-bit indexes
Figure 8. Wall-clock time in nanoseconds per element to shuffle arrays of random integers using either regular shuffles or buffered shuffles with a buffer of size 256 (see Algorithm 6).

We also compare Algorithm 5 to the floating-point approach we described in § 1 where we represent the random number as a floating-point value in which we multiply by to get a number in . As illustrated in Fig. 11, the floating-point approach is slightly slower (10%–30%) whether we use 32-bit or 64-bit indexes. Moreover, it introduces a small bias and it is limited to in the 32-bit case and to in the 64-bit case.

(a) 32-bit indexes
(b) 64-bit indexes
Figure 11. Wall-clock time in nanoseconds per element to shuffle arrays of random integers using either Algorithm 5 or an approach based on floating-point multiplications.

6. Conclusion

We find that the algorithm often used in OpenBSD and macOS (through the arc4random_uniform function) requires two divisions per random number generation. It is also the slowest in our tests. The Java approach that often requires only one division, can be faster. We believe that it should be preferred to the OpenBSD algorithm.

As we have demonstrated, we can use nearly no division at all and at most one in the worst case. Avoiding divisions can multiply the speed of unbiased random shuffling functions on x64 processors.

For its new random shuffle function, the Go programming language adopted our proposed approach (Snyder2017). The Go authors justify this decision by the fact that it results in 30% better speed compared with the application of the OpenBSD approach (Algorithm 3).

Our results are only relevant in the context where the generation of random integers is fast compared with the latency of a division operation. In the case where the generation of random bits is likely the bottleneck, other approaches would be preferable (Bacher:2017:GRP:3040971.3009909). Moreover, our approach may not be applicable to specialized processors such as Graphics Processing Units (GPUs) that lack support for the computation of the full multiplication (Langdon:2009:FHQ:1570256.1570353).333The CUDA API provided by Nvidia offers relevant intrinsic functions such as __umul64hi.

Appendix A Code Samples

// returns value in [0,s)
// random64 is a function returning random 64-bit words
uint64_t openbsd(uint64_t s, uint64_t (*random64)(void)) {
  uint64_t t = (-s) % s;
  uint64_t x;
  do {
    x = random64();
  } while (x < t);
  return x % s;
}\end{lstlisting}\end{minipage}
\begin{minipage}{\linewidth}\begin{lstlisting}
uint64_t java(uint64_t s, uint64_t (*random64)(void)) {
  uint64_t x = random64();
  uint64_t r = x % s;
  while (x - r > UINT64_MAX - s + 1) {
    x = random64();
    r = x % s;
  }
  return r;
}\end{lstlisting}\end{minipage}
\begin{minipage}{\linewidth}
\begin{lstlisting}
uint64_t nearlydivisionless(uint64_t s, uint64_t (*random64)(void)) {
  uint64_t x = random64();
  __uint128_t m = (__uint128_t) x * (__uint128_t) s;
  uint64_t l = (uint64_t) m;
  if (l < s) {
   uint64_t t = -s % s;
    while (l < t) {
      x = random64();
      m = (__uint128_t) x * (__uint128_t) s;
      l = (uint64_t) m;
    }
  }
  return m >> 64;
}\end{lstlisting}
\end{minipage}
\begin{acks}
The work is supported by the \grantsponsor{NSERC}{Natural Sciences and Engineering Research Council of Canada}{}
under grant
 number~\grantnum{NSERC}{RGPIN-2017-03910}. The author would like to thank R.~Startin and J.~Epler for independently reproducing the experimental results and providing valuable feedback.
\end{acks}
%\bibliographystyle{ACM-Reference-Format-Journals}
\bibliographystyle{ACM-Reference-Format}
\bibliography{rangedrng}
\end{document}