Fast, unbiased, integer pseudo random generator with arbitrary bounds

Question

For a monte carlo integration process, I need to pull a lot of random samples from a histogram that has N buckets, and where N is arbitrary (i.e. not a power of two) but doesn't change at all during the course of the computation.

By a lot, I mean something on the order of 10^10, 10 billions, so pretty much any kind of lengthy precomputation is likely worth it in the face of the sheer number of samples).

I have at my disposal a very fast uniform pseudo random number generator that typically produces unsigned 64 bits integers (all the ints in the discussion below are unsigned).

The naive way to pull a sample : histogram[ prng() % histogram.size() ]

The naive way is very slow: the modulo operation is using an integer division (IDIV) which is terribly expensive and the compiler, not knowing the value of histogram.size() at compile time, can't be up to its usual magic (i.e. http://www.azillionmonkeys.com/qed/adiv.html)

As a matter of fact, the bulk of my computation time is spent extracting that darn modulo.

The slightly less naive way: I use libdivide (http://libdivide.com/) which is capable of pulling off a very fast "divide by a constant not known at compile time".

That gives me a very nice win (25% or so), but I have a nagging feeling that I can do better, here's why:

• First intuition: libdivide computes a division. What I need is a modulo, and to get there I have to do an additional mult and a sub : mod = dividend - divisor*(uint64_t)(dividend/divisor). I suspect there might be a small win there, using libdivide-type techniques that produce the modulo directly.

• Second intuition: I am actually not interested in the modulo itself. What I truly want is to efficiently produce a uniformly distributed integer value that is guaranteed to be strictly smaller than N.

The modulo is a fairly standard way of getting there, because of two of its properties:

• A) mod(prng(), N) is guaranteed to be uniformly distributed if prng() is

• B) mod(prgn(), N) is guaranteed to belong to [0,N[

But the modulo is/does much more that just satisfy the two constraints above, and in fact it does probably too much work.

All need is a function, any function that obeys constraints A) and B) and is fast.

So, long intro, but here comes my two questions:

• Is there something out there equivalent to libdivide that computes integer modulos directly ?

• Is there some function F(X, N) of integers X and N which obeys the following two constraints:

  • If X is a random variable uniformly distributed then F(X,N) is also unirformly distributed
  • F(X, N) is guranteed to be in [0, N[

(PS : I know that if N is small, I do not need to cunsume all the 64 bits coming out of the PRNG. As a matter of fact, I already do that. But like I said, even that optimization is a minor win when compare to the big fat loss of having to compute a modulo).

Edit : prng() % N is indeed not exactly uniformly distributed. But for N large enough, I don't think it's much of problem (or is it ?)

Edit 2 : prng() % N is indeed potentially very badly distributed. I had never realized how bad it could get. Ouch. I found a good article on this : http://ericlippert.com/2013/12/16/how-much-bias-is-introduced-by-the-remainder-technique

Answer 1

Under the circumstances, the simplest approach may work the best. One extremely simple approach that might work out if your PRNG is fast enough would be to pre-compute one less than the next larger power of 2 than your N to use as a mask. I.e., given some number that looks like 0001xxxxxxxx in binary (where x means we don't care if it's a 1 or a 0) we want a mask like 000111111111.

From there, we generate numbers as follows:

1. Generate a number
2. and it with your mask
3. if result > n, go to 1

The exact effectiveness of this will depend on how close N is to a power of 2. Each successive power of 2 is (obviously enough) double its predecessor. So, in the best case N is exactly one less than a power of 2, and our test in step 3 always passes. We've added only a mask and a comparison to the time taken for the PRNG itself.

In the worst case, N is exactly equal to a power of 2. In this case, we expect to throw away roughly half the numbers we generated.

On average, N ends up roughly halfway between powers of 2. That means, on average, we throw away about one out of four inputs. We can nearly ignore the mask and comparison themselves, so our speed loss compared to the "raw" generator is basically equal to the number of its outputs that we discard, or 25% on average.

Answer 2

If you have fast access to the needed instruction, you could 64-bit multiply prng() by N and return the high 64 bits of the 128-bit result. This is sort of like multiplying a uniform real in [0, 1) by N and truncating, with bias on the order of the modulo version (i.e., practically negligible; a 32-bit version of this answer would have small but perhaps noticeable bias).

Another possibility to explore would be use word parallelism on a branchless modulo algorithm operating on single bits, to get random numbers in batches.

Answer 3

Libdivide, or any other complex ways to optimize that modulo are simply overkill. In a situation as yours, the only sensible approach is to

1. ensure that your table size is a power of two (add padding if you must!)

2. replace the modulo operation with a bitmask operation. Like this:

   size_t tableSize = 1 << 16;
   size_t tableMask = tableSize - 1;
   
   ...
   
   histogram[prng() & tableMask]
   

A bitmask operation is a single cycle on any CPU that is worth its money, you can't beat its speed.

--

Note:
I don't know about the quality of your random number generator, but it may not be a good idea to use the last bits of the random number. Some RNGs produce poor randomness in the last bits and better randomness in the upper bits. If that is the case with your RNG, use a bitshift to get the most significant bits:

size_t bitCount = 16;

...

histogram[prng() >> (64 - bitCount)]

This is just as fast as the bitmask, but it uses different bits.

Answer 4

You could extend your histogram to a "large" power of two by cycling it, filling in the trailing spaces with some dummy value (guaranteed to never occur in the real data). E.g. given a histogram

[10, 5, 6]

extend it to length 16 like so (assuming -1 is an appropriate sentinel):

[10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, -1]

Then sampling can be done via a binary mask histogram[prng() & mask] where mask = (1 << new_length) - 1, with a check for the sentinel value to retry, that is,

int value;
do {
    value = histogram[prng() & mask];
} while (value == SENTINEL);

// use `value` here

The extension is longer than necessary to make retries unlikely by ensuring that the vast majority of the elements are valid (e.g. in the example above only 1/16 lookups will "fail", and this rate can be reduced further by extending it to e.g. 64). You could even use a "branch prediction" hint (e.g. __builtin_expect in GCC) on the check so that the compiler orders code to be optimal for the case when value != SENTINEL, which is hopefully the common case.

This is very much a memory vs. speed trade-off.

Answer 5

Just a few ideas to complement the other good answers:

1. What percent of time is spent in the modulo operation, and how do you know what that percent is? I only ask because sometimes people say something is terribly slow when in fact it is less than 10% of the time and they only think it's big because they're using a silly self-time-only profiler. (I have a hard time envisioning a modulo operation taking a lot of time compared to a random number generator.)

2. When does the number of buckets become known? If it doesn't change too frequently, you can write a program-generator. When the number of buckets changes, automatically print out a new program, compile, link, and use it for your massive execution. That way, the compiler will know the number of buckets.

3. Have you considered using a quasi-random number generator, as opposed to a pseudo-random generator? It can give you higher precision of integration in much fewer samples.

4. Could the number of buckets be reduced without hurting the accuracy of the integration too much?

Answer 6

The non-uniformity dbaupp cautions about can be side-stepped by rejecting&redrawing values no less than M*(2^64/M) (before taking the modulus).
If M can be represented in no more than 32 bits, you can get more than one value less than M by repeated multiplication (see David Eisenstat's answer) or divmod; alternatively, you can use bit operations to single out bit patterns long enough for M, again rejecting values no less than M.
(I'd be surprised at modulus not being dwarfed in time/cycle/energy consumption by random number generation.)

Answer 7

To feed the bucket, you may use std::binomial_distribution to directly feed each bucket instead of feeding the bucket one sample by one sample:

Following may help:

int nrolls = 60; // number of experiments
const std::size_t N = 6;
unsigned int bucket[N] = {};

std::mt19937 generator(time(nullptr));

for (int i = 0; i != N; ++i) {
    double proba = 1. / static_cast<double>(N - i);
    std::binomial_distribution<int> distribution (nrolls, proba);
    bucket[i] = distribution(generator);
    nrolls -= bucket[i];
}

Live example

Answer 8

Instead of integer division you can use fixed point math, i.e integer multiplication & bitshift. Say if your prng() returns values in range 0-65535 and you want this quantized to range 0-99, then you do (prng()*100)>>16. Just make sure that the multiplication doesn't overflow your integer type, so you may have to shift the result of prng() right. Note that this mapping is better than modulo since it's retains the uniform distribution.

Answer 9

Thanks everyone for you suggestions.

First, I am now thoroughly convinced that modulo is really evil.
It is both very slow and yields incorrect results in most cases.

After implementing and testing quite a few of the suggestions, what
seems to be the best speed/quality compromise is the solution proposed
by @Gene:

1. pre-compute normalizer as:

   auto normalizer = histogram.size() / (1.0+urng.max());

2. draw samples with:

   return histogram[ (uint32_t)floor(urng() * normalizer);

It is the fastest of all methods I've tried so far, and as far as I can tell,
it yields a distribution that's much better, even if it may not be as perfect
as the rejection method.

Edit: I implemented David Eisenstat's method, which is more or less the same as Jarkkol's suggestion : index = (rng() * N) >> 32. It works as well as the floating point normalization and it is a little faster (9% faster in fact). So it is my preferred way now.