Lemire’s Fast Alternative to Modulo Reduction

Daniel Lemire’s fast alternative to modulo reduction maps an uniformly distributed unsigned integer into a bounded range \([0, s)\) using fixed-point multiplication rather than expensive hardware division. By treating the random integer as a fractional value between 0 and 1, multiplying it by the target bound, and extracting the most significant bits of the resulting product, modern processors can generate bounded random numbers in only a few clock cycles.

The Cost of Traditional Modulo Reduction

The standard approach to generate a random number within a bound \(s\) uses the modulo operator:

\[\text{result} = x \pmod s\]

While simple, this method relies on integer division instructions (such as DIV or IDIV on x86 architectures). Division is one of the slowest arithmetic operations on modern CPUs, often requiring 10 to 40 clock cycles. Modulo operations also introduce “modulo bias” unless combined with rejection sampling, which traditionally introduces even more division operations to compute rejection thresholds.

The Fixed-Point Concept

Lemire’s algorithm views an \(L\)-bit random integer \(x\) (where \(0 \le x < 2^L\)) as a fractional value in the interval \([0, 1)\):

\[\text{fraction} = \frac{x}{2^L}\]

To scale this value to the desired range \([0, s)\), you multiply the fraction by \(s\):

\[\text{scaled value} = \frac{x \times s}{2^L}\]

The bounded integer is the floor of this expression:

\[\text{result} = \left\lfloor \frac{x \times s}{2^L} \right\rfloor\]

Binary Implementation Without Division

In binary arithmetic, division by \(2^L\) is equivalent to shifting right by \(L\) bits. When working with standard register widths (such as \(L = 32\) bits):

  1. Full-Precision Multiplication: Compute the full 64-bit product \(P = x \times s\), where both \(x\) and \(s\) are 32-bit unsigned integers.
  2. Extracting the Quotient (High Bits): The upper 32 bits of \(P\) correspond to \(\lfloor P / 2^{32} \rfloor\). This represents the scaled integer in the range \([0, s)\).
  3. Extracting the Remainder (Low Bits): The lower 32 bits of \(P\) represent the fractional remainder: \(P \pmod{2^{32}}\).

Modern 64-bit architectures can perform a \(32\text{-bit} \times 32\text{-bit} \to 64\text{-bit}\) multiplication in 3 to 4 clock cycles, after which a right-shift (or reading the upper register) immediately yields the bounded result without invoking a division instruction.

// Fast bounded random generation (biased version)
uint32_t fast_range(uint32_t x, uint32_t s) {
    uint64_t full_product = (uint64_t)x * (uint64_t)s;
    return (uint32_t)(full_product >> 32);
}

Eliminating Bias with Efficient Rejection Sampling

Directly taking the upper bits yields a slight bias if \(2^L\) is not evenly divisible by \(s\). Lemire’s method eliminates this bias by evaluating the lower bits (the fractional portion).

A bias occurs only when the fraction falls into the incomplete zone of size \(2^L \pmod s\).

  1. Calculate the Threshold: The threshold \(t = 2^L \pmod s\) represents the boundary below which results must be rejected. This can also be written as \(t = -s \pmod s\).
  2. Check the Fractional Part: Let \(l\) be the lower \(L\) bits of \(P = x \times s\).
  3. Conditional Rejection:
    • If \(l \ge t\), the result in the upper bits is completely unbiased.
    • If \(l < t\), the number is rejected, a new random \(x\) is drawn, and the check repeats.
// Fast bounded random generation (unbiased version)
uint32_t unbiased_bounded_rand(uint32_t (*rng)(void), uint32_t s) {
    uint32_t x = rng();
    uint64_t full_product = (uint64_t)x * (uint64_t)s;
    uint32_t low_bits = (uint32_t)full_product;

    if (low_bits < s) {
        uint32_t threshold = -s % s;
        while (low_bits < threshold) {
            x = rng();
            full_product = (uint64_t)x * (uint64_t)s;
            low_bits = (uint32_t)full_product;
        }
    }
    return (uint32_t)(full_product >> 32);
}

Because the initial branch low_bits < s is false for almost all random inputs (the probability of entering the branch is only \(s / 2^{32}\)), the division to compute threshold is almost never executed. The algorithm operates nearly 100% of the time purely through multiplication and bit shifts, making it significantly faster than standard modulo-based bounded generators.