Unbiased Random Integers from a Finite Random Source

Every randomized algorithm that indexes an array, picks a pivot, or chooses an edge needs a uniform integer in a range that is almost never a power of two. The source supplies bits: a generator returns a 32-bit or 64-bit word, each value equally likely. The conversion from that word to a number in [0, n) is where uniformity is usually lost, and the loss is not random noise that averages out. It is a fixed multiplicative distortion of the output distribution, and it survives any amount of repetition.

The remainder leaves part of the range heavier

Take a 32-bit source, uniform on [0, 2^32), and compute x % 3. Since 2^32 = 4294967296 = 3 × 1431655765 + 1, the residue 0 occurs 1431655766 times while residues 1 and 2 occur 1431655765 times each. The relative excess is about 7 × 10^-10. That is small enough to escape notice in a test suite and large enough to matter in a Monte Carlo estimate that sums a billion draws, because the distortion is systematic rather than statistical.

The size of the distortion depends on how close n is to the top of the source range. Write q = floor(2^k / n) and r = 2^k mod n, with k the word width. Every residue class receives q preimages, and the first r classes receive q + 1. The relative excess of those classes is 1/q. For k = 32 and n = 3 × 2^30 = 3221225472, the quotient is 1 and the remainder is 2^30, so a third of the output values are twice as likely as the rest. For n just above 2^31, the same effect appears with roughly a factor of two between the two halves. A one-in-ten-billion skew and a fifty-percent skew are the same bug at different values of n.

The general statement is worth keeping: modulo reduction is exact only when n divides 2^k. Powers of two are safe. Nothing else is.

Rejection restores exactness with a bounded number of draws

The repair is to discard the part of the source range that cannot be divided evenly. Compute limit = 2^k - (2^k mod n). Draw x; if x < limit, return x % n, otherwise draw again. Every residue now has exactly limit / n preimages among the accepted values, so the output is uniform by construction rather than by approximation.

The cost is a random number of draws, and it is bounded in a way that is easy to state and easy to get wrong. The acceptance probability is limit / 2^k = 1 - r/2^k. Two cases cover every n in [1, 2^k]. If n > 2^{k-1} then the quotient q is 1 and r = 2^k - n < 2^{k-1}, so acceptance exceeds 1/2. If n ≤ 2^{k-1} then r < n ≤ 2^{k-1}, so acceptance again exceeds 1/2. Acceptance is therefore always strictly above one half, and the expected number of draws is below 2. The tail is geometric: the probability that a call needs more than t draws falls off as 2^-t.

Two implementation details decide whether the bound holds. The comparison must be against limit, not against n, and the loop must redraw rather than wrap or clamp. Clamping the rejected value to n - 1 reintroduces exactly the bias the rejection was meant to remove, concentrated on the last value of the range.

One multiply instead of a division

Rejection sampling as described needs a 64-bit division to compute limit, and division is the expensive integer operation on common processors: on the order of tens of cycles for 64-bit operands, against a few cycles for a multiply. Daniel Lemire’s multiply-and-reject construction moves the division off the hot path. Multiply the draw x by n in 128-bit arithmetic, take the high 64 bits as the candidate result, and inspect the low 64 bits to decide whether the candidate is biased.

The reasoning is a fixed-point argument. The product x · n / 2^64 distributes the 2^64 possible draws as evenly as possible across [0, n): each output value receives either floor(2^64 / n) or ceil(2^64 / n) preimages. The values that receive the extra preimage are precisely those whose low half falls below the threshold t = 2^64 mod n. Rejecting draws with low < t removes them, and the remaining distribution is exactly uniform. The expected number of draws is again below 2, for the same reason as before.

#include <stdint.h>

/* n must be in [1, 2^64). next_u64() returns a uniform 64-bit word. */
uint64_t bounded_lemire(uint64_t n) {
    uint64_t x = next_u64();
    __uint128_t m = (__uint128_t)x * n;
    uint64_t low = (uint64_t)m;
    if (low < n) {
        /* 2^64 mod n, computed as (-n) % n in unsigned arithmetic */
        uint64_t threshold = (uint64_t)(-n) % n;
        while (low < threshold) {
            x = next_u64();
            m = (__uint128_t)x * n;
            low = (uint64_t)m;
        }
    }
    return (uint64_t)(m >> 64);
}

The guard if (low < n) is not an optimization detail to be dropped: the threshold is only relevant when it can exceed the low half, and when low ≥ n the candidate is already known to be exact. The threshold itself is 2^64 mod n because unsigned negation computes 2^64 - n, and the remainder of that value modulo n is the same as the remainder of 2^64.

For small ranges there is a third option: draw k = ceil(log2 n) bits and reject any value at or above n. For n = 6, three bits give eight outcomes and six are accepted, so the acceptance rate is 3/4 and the expected number of bit-draws is 8/3. The method needs a source that can produce an arbitrary number of bits, which is a different interface from a word-returning generator, but it avoids both division and multiplication.

Bias that survives correct range reduction

A correct [0, n) conversion does not make the output uniform if the source is not. Three failure modes sit outside the conversion step.

The first is floating-point scaling. Computing floor(u * n) for a u in [0, 1) drawn from a double is exact only when the product is representable, and a double carries 53 bits of mantissa. For n above 2^53, distinct integers collapse onto the same double and some outputs become unreachable. For small n, the rounding of u itself leaves a small periodic distortion. A generator that can return exactly 1.0 produces n, which is out of range and indexes one past the end of an array.

The second is a weak source. A linear congruential generator has a period in its low-order bits: bit i of the state alternates with period 2^{i+1}, so the low bit alternates every step. Taking x % n for even n from such a generator inherits that structure. RFC 4086 discusses the entropy requirements for cryptographic use and the reasons a short seed cannot be stretched into an unpredictable stream; the same reasoning applies to any simulation whose result must not be predictable.

The third is compounding. A shuffle performs n - 1 bounded draws, and a pivot choice performs one per partition. A per-draw relative bias of 10^-9 is negligible once. Applied at every step of a procedure whose output is a permutation, it becomes a preference for particular orderings that no amount of additional sampling removes. That is why the conversion routine deserves exactness rather than an argument that the error is small.

How to check a conversion

For a source with a small word width, the output distribution can be enumerated exactly rather than estimated. A 16-bit source has 65536 possible values; iterating all of them and counting the outputs of the conversion gives the exact multiset of results, with no sampling error and no measurement of any machine. A correct conversion produces counts that differ by at most one, and the values that receive the extra count are the ones the rejection step was designed to leave in place. A conversion with a bug produces a count pattern that follows the residue structure of n, which is visible immediately at n = 3 and at n = 2^{15} + 1.

The same enumeration answers the cost question. Counting how many of the 65536 values are rejected gives the exact acceptance rate, which should match limit / 2^k to the last unit. Both numbers are properties of the code, not of a machine, so they can be asserted in a test that runs in microseconds and never flakes.

Downstream, this routine is the single point where uniformity enters. Randomized quicksort’s expected bound assumes a uniformly chosen pivot, which is one bounded draw per partition; the analysis in Randomized Quicksort is a statement about that draw and not about the generator behind it. The shuffle in Uniform Random Permutations needs n - 1 such draws, each with a different upper bound. Get the primitive right and those analyses apply; get it wrong and they describe an algorithm nobody is running.