# Isn't std::bernoulli_distribution inefficient? Designing a bit-parallel Bernoulli generator

C++11 has a convenient Bernoulli RNG, illustrated at http://en.cppreference.com/w/cpp/numeric/random/bernoulli_distribution . However, distilling an entire random integer into a single random bit seems inefficient when the expectation parameter $p$ is rational with a small or power-of-two denominator. Is there a reasonably fast way to generate 32 random Bernoulli bits at once in such cases? My application uses long streams of bits, so I can keep track of statistics if needed (but this would consume runtime).

• Explore the C++ library. Other functions generate more random bits. You can use random_device or one of the actual PRNGs, such as Mersenne Twister, which (in its 64 bit version) should generate 64 bits at once. Or you could uniform_int_distribution with min = 0 and max = 0xffffffffffffffffull. Sep 23 '13 at 6:59
• For the RNG you want to use, does extracting some bits of the RNG output still behave sufficiently randomly? There might be correlations between bit $i$ and bit $i+8$, say, which may make the output less random than you would wish. Sep 23 '13 at 14:15
• @Yuval - the question is how to bias the expectation of individual bits in such a 64-bit integer. For example, if I take two random 64-bit integers and perform bitwise AND, I will get 64 Bernoulli bits with expectation 1/4. Likewise, I can get any $p=1/2^k$ by AND'ing $k$ random integers Sep 24 '13 at 8:14
• @AndrásSalamon - it seems that the bits of random unsigneds (from the full range) should be i.i.d. OTOH, I could tolerate some autocorrelation in the Bernoulli bits, if this buys something useful. Sep 25 '13 at 4:01

Actually, generating a pseudorandom integer is very fast, so the simple approach will likely be quite efficient. As always, beware premature optimization. Benchmark it first before assuming it'll be too slow.

If you have a power-of-two denominator, it's easy to come up with a procedure that doesn't require as many random bits -- but in practice this doesn't matter, because with most PRNGs, it's just as fast to generate an entire word as to generate a few pseudorandom bits.

• Generating pseudo-random numbers is quite fast for many of the common approaches. This doesn't necessarily hold if the random numbers are computed using an actual source of entropy. The OP's notion that if a 32-bit random number is generated then one can do extract more than 1 bit, is sound. Given fast bit operations, it should be possible to extract some number of Bernoulli trials that is a power of 2 with low overhead, saving some time. But agree with your overall point -- beware premature optimization; this part of the code is unlikely to be the bottleneck. Sep 23 '13 at 14:07
• @AndrásSalamon, there is no need for an actual source of entropy for every random number; all you need is a source of entropy for the seed to the PRNG. (If you want especially high-quality pseudorandom numbers, you can use a crypto-strength PRNG, which is still extremely fast.)
– D.W.
Sep 23 '13 at 17:21

If the denominator is a small power of two, say $p = A/2^B$, then you can divide the 64-bit integer into parts of length $B$, and compare each one to $A$. A similar idea works for $p = A/B$, only this time you think of your integer as being written in base $B$, and you have to worry about "extra digits": if $2^{64} = c_kB^k + c_{k-1} B^{k-1} + \cdots + c_0$, you should reject if the input is at least $c_kB^k$, since otherwise the lower-order digits wouldn't be uniform. These methods are not optimal, but might be faster than more randomness-efficient methods.

If you are generating a lot of random bits, each one should cost you about $h(p) = -p\log_2 p -(1-p)\log_2 (1-p)$ bits (less bits the more biased your bits are!). One method to achieve this is "arithmetic decoding" - you think of your random stream as infinite, and encoding a number $x \in [0,1]$. If $x < p$ then the first output bit is $1$, and you continue with $x/p$. If $x \geq p$ then the first output bit is $0$, and you continue with $(x-p)/(1-p)$. Of course, you don't need all the random stream to decide which of these options is true. Just write $p$ in binary and compare it to successive bits of $x$. (The subsequent bits are more difficult this way.)

• Thanks - this is relevant, but for $A/2^B$, there seem to be faster recipes in some cases. For example, if I need expectation 3/4, I would take two 64-bit integers $x$ and $y$, then compute $x$|$y$. Interestingly, I am using 2 random bits per 1/2 bit entropy generated. Sep 25 '13 at 4:23