UPDATE: See this blog post.

Let’s say you need a random number from 0 to 2 (0, 1, or 2). If you do this:

`listing-1.c`

`[GitHub]`

`[Download]`

your code is broken, thanks to modulo bias. You will always get modulo bias (however slight, or great, it may be) whenever `RAND_MAX + 1`

of your `rand()`

function is not evenly divisible by your modulus `n`

.

Continuing with our example, let’s say `rand()`

’s `RAND_MAX`

is 4, so there are 5 different possible values that `rand()`

can generate (0, 1, 2, 3, or 4). Now rand() will return each value 20% of the time, assuming that `rand()`

has very good distribution from 0 to `RAND_MAX`

.^{1} But since we take the modulo by 3, **we get 0 40% of the time (0 and 3)!** This fact alone should make you worried. Visually it is as follows:

The following code is correct:

```
/* Listing 2 */
int n, x, rand_limit, rand_excess;
n = 3;
rand_excess = (RAND_MAX + 1) % n;
rand_limit = RAND_MAX - rand_excess;
while (x = rand() > rand_limit) {};
return x % n;
```

`listing-2.c`

`[GitHub]`

`[Download]`

So, if `RAND_MAX`

is 4 as in the previous example, `rand_excess`

is \[
(4 + 1) \mbox{ mod 3}\\
= 5 \mbox{ mod 3}\\
= 2
\]. Then `rand_limit`

becomes \(4 - 2 = 2\). The `while`

statement then throws out the values 3 and 4 (i.e., x is only allowed to be 0, 1, or 2). Then, we return the modulo expression `x % n`

(which is redundant here onlyu because `RAND_MAX`

is very low for the sake of our example).

A problem might exist if `RAND_MAX`

is already the maximum value allowed by the type that `rand()`

returns. E.g., if `rand()`

returns a 64-bit unsigned integer (in C99, this is the `uint64_t`

type) , and `RAND_MAX`

is already set to \(2^{64} - 1\), then `RAND_MAX - 1`

would wrap back around to 0, and `rand_excess`

would likewise be 0. To avoid this, you can use the alternate expression `rand_excess = (RAND_MAX % n) + 1;`

for setting `rand_excess`

:

```
/* Listing 3 */
int n, x, rand_limit, rand_excess;
n = 3;
rand_excess = (RAND_MAX % n) + 1; /* only difference vs listing-2.c */
rand_limit = RAND_MAX - rand_excess;
while (x = rand() > rand_limit) {};
return x % n;
```

`listing-3.c`

`[GitHub]`

`[Download]`

This way, you ensure that you shrink the value of `RAND_MAX`

to something smaller first, before adding 1 to it as in Listing 2. Unfortunately, there still remains a problem with this workaround: there is a remote chance that `rand_excess`

will be equal to `n`

, needlessly reducing the size of `rand_limit`

(and thus throwing away perfectly bias-free numbers generated by `rand()`

). For example, say `RAND_MAX`

is 8 and `n`

is 3. Then `rand_excess`

is \[
(8 \mbox{ mod 3}) + 1\\
= 2 + 1\\
= 3
\]. Now `rand_limit`

is \(8 - 3 = 5\). But `RAND_MAX`

of 8 is already valid because there are 9 possible values 0 through 8, and there is no modulo bias to begin with (since \(9 \mbox{ mod 3} = 0\))! To get around this, we can do one more modulo operation `rand_excess = ((RAND_MAX % n) + 1) % n;`

to avoid the case where `rand_excess`

could equal `n`

:

```
/* Listing 4 */
int n, x, rand_limit, rand_excess;
n = 3;
rand_excess = ((RAND_MAX % n) + 1) % n; /* only difference vs listing-3.c */
rand_limit = RAND_MAX - rand_excess;
while (x = rand() > rand_limit) {};
return x % n;
```

`listing-4.c`

`[GitHub]`

`[Download]`

So if `RAND_MAX`

is 8 and `n`

is 3, then `rand_excess`

is \[
((8 \mbox{ mod 3}) + 1) \mbox{ mod 3}\\
= (2 + 1) \mbox{ mod 3}\\
= 3 \mbox{ mod 3}\\
= 0
\] and `rand_limit`

would not be reduced by any amount.^{2}

If you actually know beforehand the exact size of `RAND_MAX`

(usually it is `INT_MAX`

or some such, determined by the type of the variable returned by `rand()`

), you should instead just define `rand_limit`

as a constant to save some CPU cycles (write out the mathematical expressions involved by hand to find the number and just hardcode it).