Generating Random Numbers under a Maximum
13 Feb 2024Recently I asked a question about generating random numbers.
My first approach to generating random numbers lower than 1000 was to take ten random bits (which can represent 0..1023).
If that random number was less than 1000, then that was the result. If it was not, I would repeat the process until the ten bits yielded a result less than 1000. (Reusing any of those ten bits introduces bias.) This process is known as sample-and-reject.
Stef pointed out that if the result was greater than 999, then we had log2(24) bits of unused entropy.
This was an awkward amount of entropy to reuse, so I started to think if I could avoid taking it in the first place….
The algorithm below is an entropy efficient way of generating unbiased random numbers, less than a maximum value “max”, from a stream of random bits.
It is more entropy-efficient than sample-and-reject if and only if “max” has a factor which is a positive power of two.
It makes use of the fact that if the maximum has a factor which is a positive power of two, then we can decide whether a number will be too big by taking a number of random bits which are fewer than the length of the maximum (in binary). The difference in the numbers of bits is the exponent in the maximum’s power-of-two factor.
Example:
If max is 1000, it has a factor of two-to-the-power-of-three (8). 1000 is a ten bit number. But we only need to take seven bits to decide if a number 0..1023 is less than 1000.
two_power == 3 log2(8)
bits == 10 log2(1024)
significant_bits == 7 10 - 3
We first generate 7 bits and shift them left three places, to give us:
r == ..._000000??_?????000
We compare this with 1000:
max == ..._00000011_11101000
So long as r >= max, we have to try again with another 7 bits.
Once r < max, we fill in the final three bits of precision:
r == ..._000000??_????????
Entropy Efficient Algorithm for Generating Random Numbers under a Maximum
const BITS_IN_U64: u32 = 64;
fn random_under_max(max: u64) -> u64 {
let mut rng = rand::thread_rng();
// How many bits do we need?
let two_power = max.trailing_zeros();
let bits = BITS_IN_U64 - max.leading_zeros();
let significant_bits = bits - two_power;
let mut r = u64::MAX;
while r >= max {
// First collect the minimal amount of random bits which
// will let us decide whether the random number is going to be less
// than max.
r = 0;
for _ in 0..significant_bits {
// Generate one random bit
let b: bool = rng.gen();
// and shift it into our imprecise random number.
r = r << 1 | b as u64;
}
// Shift the inexact random number into the correct position.
r = r << two_power;
}
// Now fill-in the precise value of this random number.
r = r >> two_power; // Scrub out the insignificant zeros.
for _ in 0..two_power {
// Generate one random bit
let b: bool = rng.gen();
// and shift it into our random number.
r = r << 1 | b as u64;
}
// Return the generated number.
r
}