Random distributions are not one-size-fits-all (part 1)

Recently, I was once again on CppCast, and once again^[1], I talked about Catch2 and C++'s standard library support for random numbers. One of the things I said was that random distributions are not one-size-fits-all.

This post is about how this applies to Lemire's algorithm, which is the state of the art for generating uniformly distributed integers in specific intervals.

Introduction

In C++, <random> separates the concern of generating random numbers between Random Number Engines, which generate integers in some predefined range, and distributions, which turn these into other numbers in the range requested by the user. This is a powerful, expert-friendly design, even though the actual standardized library leaves a lot to be desired.

Notably, the distributions are specified as types and are allowed to be stateful. This means that a distribution could be written to optimize for repeatedly generating numbers from the same instance rather than generating a single number as fast as possible^[2].

This post will look specifically at Lemire's algorithm for generating uniformly distributed integers. This is because I recently implemented reproducible uniform integer distribution for Catch2 using Lemire's algorithm, so I can use my own code for the benchmarks^[3].

Lemire's algorithm is the current state of the art for generating uniformly distributed random integers in an interval [0, max)^[4]. I won't go over the proof of the algorithm in this post; if you want to know why it works, go read the paper, but we will look at a pseudo-implementation.

uint32_t lemire( uint32_t max ) {
  uint64_t long_mult = random32() * uint64_t(max);
  uint32_t low_bits = uint32_t(long_mult);
  // Cheap probabilistic check
  if (low_bits < max) {
    uint32_t rejection_threshold = -max % max;
    // Precise check
    while (low_bits < rejection_threshold) {
      long_mult = random32() * uint64_t(max);
      low_bits = uint32_t(long_mult);
    }
  }
  return long_mult >> 32;
}

The core trick is in the probabilistic check guarding the expensive modulo operation used to calculate precise rejection threshold. The chance that we will need the modulo is \( \frac{\textbf{max}}{2^{32}} \), which translates to "basically zero" for a lot of the common use cases^[5].

But we can also end up computing rejection_threshold on every invocation. And if you look at how we compute it, you will notice that it does not depend on the generated random number, only on max. This is where a stateful object can optimize the code for generating multiple numbers; it can compute rejection_threshold once and use it every time it generates a new number.

Will precomputing the threshold pay off? Let's write some benchmarks and see.

Benchmarks

Please note that while this post has benchmark results for uint32_t and uint64_t, the latter are not representative of an optimized implementation. Lemire's algorithm uses extended multiplication, but the used implementation from Catch2 uses naive long-multiplication. This means that the results for uint64_t are much slower than they would be with optimized implementation. Thus, they can show the trend, but the performance difference between the two implementations is smaller than it should be.

Setup

All the code for the benchmarks is publicly available on GitHub. I took the code for the RNG and Lemire's algorithm from Catch2, the benchmarks themselves are done using Catch2's micro benchmarking support.

There are three benchmarks:

1. Stacked heavily for reuse.
   The larger the target interval, the greater the chance Lemire's algorithm will need to compute modulo. Thus, we will generate numbers in interval [0, numeric_limits<type> - 1)^[6], and generate lot of them.
2. Stacked heavily against reuse.
   As long as we never reuse the distribution object instance, there is no advantage to precomputing the modulo. If we also always create the distribution object with a new bound, this corresponds to how it would be used during a Fisher-Yates shuffle.
3. Benchmark based on use in Catch2.
   Catch2 internally uses uniform int distribution in two places. Both places keep reusing the same distribution but with significantly different bounds and expected number of invocations.

For every benchmark, we have three contenders:

1. Just the RNG
   This does not generate the required numbers, but it acts as a baseline to tell us what the overhead of the distribution algorithm is.
2. Plain implementation of Lemire's algorithm
3. Implementation of Lemire's algorithm with precomputed thresholds

The numbers in this post were taken from my personal machine with a Ryzen 7 4750 CPU, using MSVC in version 19.38.33133 (this MSVC version comes installed with VS 17.8.3) and compiling for CMake-generated Release configuration^[7].

Results

Benchmark stacked for reuse

uint32_t

uint64_t

As we can see, when we stack the deck towards the reuse-optimized implementation, the reuse-optimized implementation performs …

…

…

…

… exactly the same as the basic implementation.

Wait what?

Benchmarking is hard: the optimizer is good at its job

As it turns out, when the optimizer sees code like this

constexpr uint32_t high = std::numeric_limits<uint32_t>::max()-1;
BENCHMARK("noreuse, iters=" + std::to_string(iters)) {
	uint32_t sum = 0;
	lemire_algorithm_no_reuse<uint32_t> dist(0, high);
	for (size_t n = 0; n < iters; ++n) {
		sum += dist(rng);
	}
	return sum;
};

it uses the fact that high is statically known to inline everything and precompute the rejection threshold during compilation, so there is no div instruction in the resulting code. Thus, the plain implementation never pays extra cost by recomputing the modulo and cannot lose.

You can play with a simplified example on Compiler Explorer. To keep the generated ASM short, the code only generates a single number, and the RNG implementation is opaque to avoid inlining it.

Since the compiler does more or less the same for the reuse-optimized implementation, neither has an advantage when the target interval is statically known. So, to measure the effect we care about, we have to hide the interval from the compiler. In the benchmark code, this is done by calling an opaque function, which returns the desired value but cannot be inlined due to lack of LTO.

Benchmark results revisited

Benchmark stacked for reuse (again)

uint32_t

uint64_t

As we can see, after we stop the optimizer from optimizing out the important part, precomputing the modulo speeds up the random number generation by \(\approx 20\%\), with the distribution itself being up to \(15\times\) faster.^[8] Remember that this benchmark is the best case for the stateful implementation, as we stacked everything in its favour.

Let's take a look at the opposite benchmark.

Benchmark stacked against reuse

uint32_t

uint64_t

Once again, the difference in the overhead of the two implementations is massive, with plain Lemire being 9 to 12 times faster. This boils down to this benchmark being an inversion of the first one. Where the first benchmark forced the plain Lemire implementation to compute modulo for every number and let the reuse variant compute it exactly once, this benchmark forces the reuse variant to always compute modulo while plain Lemire avoids the computation most of the time.

The difference in the total runtime is much more significant this time around, with the reuse variant being 60-80% slower. This is a terrible result and shows that the reuse variant does not make a good default.

Benchmark(s) from Catch2

The two places in Catch2 needing uniformly generated integers are generating uniformly distributed floating point numbers^[9], and resampling benchmark measurements. The interval for the former depends on the range the user asks for, but we will use \([0, 33554430]\) and \([0, 18014398509481982]\) for the benchmark^[10]. These numbers present the upper bound on what we can expect and thus are most favourable for the reuse implementation. The number of iterations is unknown, but generally, I would expect it to be low (\(<1000\)).

Resampling benchmark measurements defaults to interval \([0, 100)\) and generates 10M random numbers with the one instance it creates^[11]. The user can customize both of these, but generally they won't.

Generating floats

uint32_t

uint64_t

Unsurprisingly, when you only want to get a single random (floating point) number, the plain Lemire algorithm has a large advantage. What is surprising is how few we need to generate for the reuse implementation to catch up and eventually gain an advantage. Even 1000 numbers should be too low to require the plain implementation to recompute the division repeatedly.

For uint64_t, we see the slooow multiplication implementation hide any other difference.

Benchmark resampling

uint32_t

uint64_t

For benchmark resampling, there are minimal differences between the two implementations. Plain Lemire implementation is very unlikely to pay the extra division cost due to the small upper bound, and the reuse implementation easily amortizes the additional cost over the many generated numbers. When generating a lot^[12] of numbers, the reuse implementation starts winning again, but only by small margins.

Conclusion

So, what did we learn from the benchmarks?

First, if the compiler can see the interval bounds during compilation time, it does not matter which of the two shown implementations we pick, because it will optimize the division away in either case. In other words, a function like this

uint32_t roll_D6(RNG& rng) {
	lemire_algorithm<uint32_t> dist(1, 6);    
    return dist(rng);
}

will always perform optimally with either implementation. But when the compiler cannot see the bounds, their performance difference can be massive. If we craft benchmarks specifically to advantage one implementation over the other, we see more than a 10x difference between the overhead imposed by the two implementations.

From how we crafted the edge case benchmarks, we also see that if we can expose only one implementation, the implementation optimized for reuse is a bad default. The benchmark favouring reuse is contrived, and as we make it less so, the advantage of the reuse implementation decreases. On the other hand, the benchmark favouring plain Lemire is based on a real-world usage pattern representing the Fisher-Yates shuffle.

We have also learned that Lemire's algorithm strongly relies on fast extended multiplication. If we do not have optimized extended multiplication for a specific width, other implementation considerations do not matter, and Lemire's algorithm becomes much slower.

From the workload-specific benchmarks, we've learned that Catch2 is indeed better off with the reuse-optimized implementation.^[13] Catch2 does not use random shuffle^[14], and while the performance improvement will be slight, it is functionally free because Catch2 had to implement its own distribution anyway.

That's it for this post. I will write at least one follow-up post about the case where we care about the statistical properties of the generated numbers but not their reproducibility. This means that we can use completely different algorithms, rather than looking into different implementation strategies for the same one.)