# PRNG for generating numbers with n set bits exactly

I'm currently writing some code to generate binary data. I specifically need to generate 64-bit numbers with a given number of set bits; more precisely, the procedure should take some $0 < n < 64$ and return a pseudo-random 64-bit number with exactly $n$ bits set to $1$, and the rest set to 0.

My current approach involves something like this:

1. Generate a pseudorandom 64-bit number $k$.
2. Count the bits in $k$, storing the result in $b$.
3. If $b = n$, output $k$; otherwise go to 1.

This works, but it seems inelegant. Is there some kind of PRNG algorithm which can generate numbers with $n$ set bits more elegantly than this?

What you need is a random number between 0 and ${ 64 \choose n } - 1$. The problem then is to turn this into the bit pattern.

This is known as enumerative coding, and it's one of the oldest deployed compression algorithms. Probably the simplest algorithm is from Thomas Cover. It's based on the simple observation that if you have a word that is $n$ bits long, where the set bits are $x_k \ldots x_1$ in most-significant bit order, then the position of this word in the lexicographic ordering of all words with this property is:

$$\sum_{1 \le i \le k} { x_i \choose i}$$

So, for example, for a 7-bit word:

$$i(0000111) = { 2 \choose 3 } + {1 \choose 2 } + {0 \choose 1} = 0$$ $$i(0001011) = { 3 \choose 3 } + {1 \choose 2 } + {0 \choose 1} = 1$$ $$i(0001101) = { 3 \choose 3 } + {2 \choose 2 } + {0 \choose 1} = 2$$

...and so on.

To get the bit pattern from the ordinal, you just decode each bit in turn. Something like this, in a C-like language:

uint64_t decode(uint64_t ones, uint64_t ordinal)
{
uint64_t bits = 0;
for (uint64_t bit = 63; ones > 0; --bit)
{
uint64_t nCk = choose(bit, ones);
if (ordinal >= nCk)
{
ordinal -= nCk;
bits |= 1 << bit;
--ones;
}
}
return bits;
}

Note that since you only need binomial coefficients up to 64, you can precompute them.

• Cover, T., Enumerative Source Encoding. IEEE Transactions on Information Theory, Vol IT-19, No 1, Jan 1973.
• Beautiful and elegant! Enumerative coding looks like something that's very useful - are there any good resources on it (preferably in textbook form)? Dec 20 '16 at 7:52
• Does this actually give better performance in practice? (Of course it depends on the speed of the RNG.) If not then there's no point in using more complex code. Dec 20 '16 at 10:58
• @Giles I interpreted this as a computer science question, since this is cs.se. I only gave source code because I happened to have it lying around from an RRR array implementation. (See, for example, alexbowe.com/rrr for an explanation of what that means.) Dec 21 '16 at 3:32
• @Gilles To follow up on your question, I implemented both my naive method and the one provided by Pseudonym in Forth. The naive method, even when using a very simple xorshift PRNG, took something in the order of 20 second per number, while Pseudonym's method was almost instantaneous. I used tables of precomputed binomials for this. Dec 22 '16 at 6:30
• @KozRoss If you generate n bit numbers, and look for numbers with k bits set, they would be quite rare if k is far away from n/2; that would explain it. Jan 22 '20 at 20:55

Very similar to Pseudonym's answer, obtained by other means.

The total number of available combinations is approachable by the stars and bars method, so it will have to be $c=\binom{64}{n}$. The total number of 64-bit numbers from which you would be trying to sample your number would be obviously much higher than that.

What you need then is a function that can lead you from a pseudorandom number $k$, ranging from $1$ to $c$, to the corresponding 64-bit combination.

Pascal's triangle can help you with that, because every node's value represents exactly the number of paths from that node to the root of the triangle, and every path can be made to represent one of the strings you are looking for, if all left turns are labeled with a $1$, and every right turn with a $0$.

So let $x$ be the number of bits left to determine, and $y$ be the number of ones left to use.

We know that $\binom{x}{y}=\binom{x-1}{y}+\binom{x-1}{y-1}$, and we can use it to properly determine the next bit of the number at each step:

$\mathtt{while}\;\;\; x>0$

$\quad \mathtt{if}\;\;\; x>y$

$\qquad \mathtt{if}\;\;\;k>\binom{x-1}{y}: \;\;\;s \leftarrow s\; + \mathtt{"1"}, \;k\leftarrow k-\binom{x-1}{y}, \;y \leftarrow y-1$

$\qquad \mathtt{else}:\; \;s \leftarrow s\; + \mathtt{"0"}$

$\quad \mathtt{else}: \;\;s \leftarrow s\; + \mathtt{"1"}, \;y \leftarrow y-1$

$\quad x \leftarrow x-1$

Another quite elegant method is to use bisection as described in this stackoverflow answer. The idea is to keep two words, one known to have at most k bits set and one known to have at least k bits set, and use randomness to move one of these towards having exactly k bits. Here is some source code to illustrate it:

word randomKBits(int k) {
word min = 0;
word max = word(~word(0)); // all 1s
int n = 0;
while (n != k) {
word x = randomWord();
x = min | (x & max);
n = popcount(x);
if (n > k)
max = x;
else
min = x;
}
return min;
}

I made a performance comparison of various methods and this one is typically the fastest unless k is known to be very small.

You can do the following :

1) Generate random number,$k$ between $1$ and $64$.

2) Set $k$ th $0$ to $1$.

3)Repeat steps 1 and 2 $n$ times

$A[]$ is $64$ bit array with all $0$s

for(i=1 to n)
{
k=ran(1,65-i) % random number between 1 and 65-i
for(x=1;x<65;x++)
{
if(A[x]==0)k--;
if(k==0)break;
}
A[x]=1;
}
• The prose doesn't seem to match your code? The code never assigns 1s to the array. Also it doesn't seem to generate a uniform distribution (and not even numbers that satisfy the constraints) when multiple ks collide Dec 20 '16 at 10:13
• @Bergi Ya forgot the $A[x]=1$ line...added it now. And multiple collision of k is handled. See first number is chosen between 1 and 64, second between 1 and "remaining" 63. So it skips the 1 while counting...see the $if(A[x]==0)k--;$ line. And it is uniform distribution. Dec 20 '16 at 12:25
• Ah, I see now. The prose algorithm didn't mention the skipping. Dec 20 '16 at 12:38
• @ArghyaChakraborty Are you using 1-based indexing there? Dec 20 '16 at 19:29
• @KozRoss Start with what happens if $i=1,k=1$ (of course $A$ will be all zeroes) So, it will check $A[1]==0$ and get $true$ meaning $k--;$ which gives $k=0$. So, sets $A[1]=1$ outside the loop. So yeah it is 1-based indexing. To make it 0 based all you have to do is change the inner $for$ to $(x=0;x<64;x++)$ Dec 21 '16 at 3:01

How about setting a random bit, if it is not already set, until k bits are set?

• Worded as it is, this seems to use an unlimited number of simpler operations. Aug 24 at 5:50