Most efficient method of generating a random number with a fixed number of bits set

Question

I need to generate a random number, but it needs to be selected from the set of binary numbers with equal numbers of set bits. E.g. choose a random byte value with exactly 2 bits set...

00000000 - no
00000001 - no
00000010 - no
00000011 - YES
00000100 - no
00000101 - YES
00000110 - YES
...

=> Set of possible numbers 3, 5, 6...

Note that this is a simplified set of numbers. Think more along the lines of 'Choose a random 64-bit number with exactly 40 bits set'. Each number from the set must be equally likely to arise.

Answer 1

Do a random selection from the set of all bit positions, then set those bits.

Example in Python:

def random_bits(word_size, bit_count):
    number = 0
    for bit in random.sample(range(word_size), bit_count):
        number |= 1 << bit
    return number

Results of running the above 10 times:

0xb1f69da5cb867efbL
0xfceff3c3e16ea92dL
0xecaea89655befe77L
0xbf7d57a9b62f338bL
0x8cd1fee76f2c69f7L
0x8563bfc6d9df32dfL
0xdf0cdaebf0177e5fL
0xf7ab75fe3e2d11c7L
0x97f9f1cbb1f9e2f8L
0x7f7f075de5b73362L

Answer 2

I have found an elegant solution: random-dichotomy.

Idea is that on average:

• and with a random number is dividing by 2 the number of set bits,
• or is adding 50% of set bits.

C code to compile with gcc (to have __builtin_popcountll):

#include <assert.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
/// Return a random number, with nb_bits bits set out of the width LSB
uint64_t random_bits(uint8_t width, uint8_t nb_bits)
{
    assert(nb_bits <= width);
    assert(width <= 64);
    uint64_t x_min = 0;
    uint64_t x_max = width == 64 ? (uint64_t)-1 : (1UL<<width)-1;
    int n = 0;

    while (n != nb_bits)
    {
        // generate a random value of at least width bits
        uint64_t x = random();
        if (width > 31)
            x ^= random() << 31;
        if (width > 62)
            x ^= random() << 33;

        x = x_min | (x & x_max); // x_min is a subset of x, which is a subset of x_max
        n = __builtin_popcountll(x);
        printf("x_min = 0x%016lX, %d bits\n", x_min, __builtin_popcountll(x_min));
        printf("x_max = 0x%016lX, %d bits\n", x_max, __builtin_popcountll(x_max));
        printf("x     = 0x%016lX, %d bits\n\n", x, n);
        if (n > nb_bits)
            x_max = x;
        else
            x_min = x;
    }

    return x_min;
}

In general less than 10 loops are needed to reach the requested number of bits (and with luck it can take 2 or 3 loops). Corner cases (nb_bits=0,1,width-1,width) are working even if a special case would be faster.

Example of result:

x_min = 0x0000000000000000, 0 bits
x_max = 0x1FFFFFFFFFFFFFFF, 61 bits
x     = 0x1492717D79B2F570, 33 bits

x_min = 0x0000000000000000, 0 bits
x_max = 0x1492717D79B2F570, 33 bits
x     = 0x1000202C70305120, 14 bits

x_min = 0x0000000000000000, 0 bits
x_max = 0x1000202C70305120, 14 bits
x     = 0x0000200C10200120, 7 bits

x_min = 0x0000200C10200120, 7 bits
x_max = 0x1000202C70305120, 14 bits
x     = 0x1000200C70200120, 10 bits

x_min = 0x1000200C70200120, 10 bits
x_max = 0x1000202C70305120, 14 bits
x     = 0x1000200C70201120, 11 bits

x_min = 0x1000200C70201120, 11 bits
x_max = 0x1000202C70305120, 14 bits
x     = 0x1000200C70301120, 12 bits

width = 61, nb_bits = 12, x = 0x1000200C70301120

Of course, you need a good prng. Otherwise you can face an infinite loop.

Answer 3

Say the number of bits to set is b and the word size is w. I would create a vector v of of length w with the first b values set to 1 and the rest set to 0. Then just shuffle v.

Answer 4

Here is another option which is very simple and reasonably fast in practice.

choose a bit at random
if it is already set
    do nothing
else
    set it
    increment count
end if

Repeat until count equals the number of bits you want set.

This will only be slow when the number of bits you want set (call it k) is more than half the word length (call it N). In that case, use the algorithm to set N - k bits instead and then flip all the bits in the result.

I bet the expected running time here is pretty good, although I am too lazy/stupid to compute it precisely right now. But I can bound it as less than 2*k... The expected number of flips of a coin to get "heads" is two, and each iteration here has a better than 1/2 chance of succeeding.

Answer 5

If you don't have the convenience of Python's random.sample, you might do this in C using the classic sequential sampling algorithm:

unsigned long k_bit_helper(int n, int k, unsigned long bit, unsigned long accum) {
  if !(n && k)
    return accum;
  if (k > rand() % n)
    return k_bit_helper(n - 1, k - 1, bit + bit, accum + bit);
  else
    return k_bit_helper(n - 1, k, bit + bit, accum);
}

unsigned long random_k_bits(int k) {
  return k_bit_helper(64, k, 1, 0);
}

The cost of the above will be dominated by the cost of generating the random numbers (true in the other solutions, also). You can optimize this a bit if you have a good prng by batching: for example, since you know that the random numbers will be in steadily decreasing ranges, you could get the random numbers for n through n-3 by getting a random number in the range 0..(n * (n - 1) * (n - 2) * (n - 3)) and then extracting the individual random numbers:

r = randint(0, n * (n - 1) * (n - 2) * (n - 3) - 1);
rn  = r % n; r /= n
rn1 = r % (n - 1); r /= (n - 1);
rn2 = r % (n - 2); r /= (n - 2);
rn3 = r % (n - 3); r /= (n - 3);

The maximum value of n is presumably 64 or 2⁶, so the maximum value of the product above is certainly less than 2²⁴. Indeed, if you used a 64-bit prng, you could extract as many as 10 random numbers out of it. However, don't do this unless you know the prng you use produces independently random bits.

Answer 6

I have another suggestion based on enumeration: choose a random number i between 1 and n choose k, and generate the i-th combination. For example, for n = 6, k = 3 the 20 combinations are:

000111
001011
010011
100011
001101
010101
100101
011001
101001
110001
001110
010110
100110
011010
101010
110010
011100
101100
110100
111000

Let's say we randomly choose combination number 7. We first check whether it has a 1 in the last position: it has, because the first 10 (5 choose 2) combinations have. We then recursively check the remaining positions. Here is some C++ code:

word ithCombination(int n, int k, word i) {
    // i is zero-based
    word x = 0;
    word b = 1;
    while (k) {
        word c = binCoeff[n - 1][k - 1];
        if (i < c) {
            x |= b;
            --k;
        } else {
            i -= c;
        }
        --n;
        b <<= 1;
    }
    return x;
}
word randomKBits(int k) {
    word i = randomRange(0, binCoeff[BITS_PER_WORD][k] - 1);
    return ithCombination(BITS_PER_WORD, k, i);
}

To be fast, we use precalculated binomial coefficients in binCoeff. The function randomRange returns a random integer between the two bounds (inclusively).

I did some timings (source). With the C++11 default random number generator, most time is spent in generating random numbers. Then this solution is fastest, since it uses the absolute minimum number of random bits possible. If I use a fast random number generator, then the solution by mic006 is fastest. If k is known to be very small, it's best to just randomly set bits until k are set.

Answer 7

Not exactly an algorithm suggestion, but just found a really neat solution in JavaScript to get random bits directly from Math.random output bits using ArrayBuffer.

//Swap var out with const and let for maximum performance! I like to use var because of prototyping ease
var randomBitList = function(n){
    var floats = Math.ceil(n/64)+1;
    var buff = new ArrayBuffer(floats*8);
    var floatView = new Float64Array(buff);
    var int8View = new Uint8Array(buff);
    var intView = new Int32Array(buff);
    for(var i = 0; i < (floats-1)*2; i++){
        floatView[floats-1] = Math.random();
        int8View[(floats-1)*8] = int8View[(floats-1)*8+4];
        intView[i] = intView[(floats-1)*2];
    }
    this.get = function(idx){
        var i = idx>>5;//divide by 32
        var j = idx%32;
        return (intView[i]>>j)&1;
        //return Math.random()>0.5?0:1;
    };
    this.getBitList = function(){
        var arr = [];
        for(var idx = 0; idx < n; idx++){
            var i = idx>>5;//divide by 32
            var j = idx%32;
            arr[idx] = (intView[i]>>j)&1;
        }
        return arr;
    }
};