# 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. – Yuval Filmus 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. – András Salamon 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 – Igor Markov 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. – Igor Markov Sep 25 '13 at 4:01

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. – Igor Markov Sep 25 '13 at 4:23