rand() example code, unnecessary check for larger than max?

Question

I've been looking into the int rand() function from <stdlib.h> in C11 when I stumbled over the following cppreference-example for rolling a six sided die.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 
int main(void)
{
    srand(time(NULL)); // use current time as seed for random generator
    int random_variable = rand();
    printf("Random value on [0,%d]: %d\n", RAND_MAX, random_variable);
 
    // roll a 6-sided die 20 times
    for (int n=0; n != 20; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + rand()/((RAND_MAX + 1u)/6); // Note: 1+rand()%6 is biased
        printf("%d ",  x); 
    }
}

Specifically this part:

[...]
        while(x > 6) 
            x = 1 + rand()/((RAND_MAX + 1u)/6); // Note: 1+rand()%6 is biased
[...]

Questions:

Why the addition of + 1u? Since rand() is [0,RAND_MAX] I'm guessing that doing rand()/(RAND_MAX/6) -> [0,RAND_MAX/(RAND_MAX/6)] -> [0,6]? And since it's integer division (LARGE/(LARGE+small)) < 1 -> 0, adding 1u gives it the required range of [0,5]?
Building on the previous question, assuming [0,5], 1 + (rand()/((RAND_MAX+1u)/6)) should only go through [1,6] and never trigger a second loop?

Been poking around to see if rand() has returned float at some point, but that seems like a pretty huge breakage towards old code? I guess the check makes sense if you add 1.0f instead of 1u making it a floating point division?

Trying to wrap my head around this, have a feeling that I might be missing something..

(P.s. This is not a basis for anything security critical, I'm just exploring the standard library. D.s)

"Why the addition of + 1u" How should we know? What's the code supposed to do? — Lundin, Oct 10 '19 at 11:46
The example is for rolling a six-sided dice, I'll add it to the question. — GlassShark, Oct 10 '19 at 11:48
_Why the addition of + 1u_, probably intended to prevent divide-by-zero. — ryyker, Oct 10 '19 at 11:55
rand() returns a number in the range of [0..RAND_MAX] inclusive, there are therefore RAND_MAX+1 possible values. Why the % operator is biased is explained [here](https://stackoverflow.com/questions/10984974/why-do-people-say-there-is-modulo-bias-when-using-a-random-number-generator). — Hans Passant, Oct 10 '19 at 12:06
`while(x > 6)`?? That's going to be true for a long time, surely? — r3mainer, Oct 10 '19 at 12:08
@squeamishossifrage No. With `RAND_MAX =32767` (MSV C) about 1 in 16k times. — chux - Reinstate Monica, Oct 10 '19 at 12:40
"should only go through [1,6] and never trigger a second loop?" --> Pretend `RAND_MAX == 9` and walk through the code. — chux - Reinstate Monica, Oct 10 '19 at 12:44

score 5 · Accepted Answer · edited Jun 20 '20 at 09:12

The code avoids bias by ensuring each possible result in [1, 6] is the output from exactly the same number of return values from rand.

By definition, rand returns int values from 0 to RAND_MAX. So there are 1+RAND_MAX possible values it can return. If 1+RAND_MAX is not a multiple of 6, then it is impossible to partition it into 6 exactly equal intervals of integers. So the code partitions it into 6 equal intervals that are as big as possible and one odd-size fragment interval. Then the results of rand are mapped into these intervals: The first six intervals correspond to results from 1 to 6, and the last interval is rejected, and the code tries again.

When we divide 1+RAND_MAX by 6, there is some quotient q and some remainder r. Now consider the result of rand() / q:

When rand produces a number in [0, q−1], rand() / q will be 0.
When rand produces a number in [q, 2q−1], rand() / q will be 1.
When rand produces a number in [2q, 3q−1], rand() / q will be 2.
When rand produces a number in [3q, 4q−1], rand() / q will be 3.
When rand produces a number in [4q, 5q−1], rand() / q will be 4.
When rand produces a number in [5q, 6q−1], rand() / q will be 5.
When rand produces a number that is 6q or greater, rand() / q will be 6.

Observe that in each of the first six intervals, there are exactly q numbers. In the seventh interval, the possible return values are in [6q, RAND_MAX]. That interval contains r numbers.

This code works by rejecting that last interval:

int x = 7;
while(x > 6) 
    x = 1 + rand()/((RAND_MAX + 1u)/6);

Whenever rand produces a number in that last fragmentary interval, this code rejects it and tries again. When rand produces a number in one of the whole intervals, this code accepts it and exits (after adding 1 so the results in x are 1 to 6 instead of 0 to 5).

Thus, every output from 1 to 6, inclusive, is mapped to from an exactly equal number of rand values.

This is the best way to produce a uniform distribution from rand in the sense that it has the fewest rejections, given we are using a scheme like this.¹ The range of rand has been split into six intervals that are as big as possible. The remaining fragmentary interval cannot be used because the remainder r is less than six, so the r unused values cannot be split evenly over the six desired values for x.

Footnote

¹ This is not necessarily the best way to use rand to generate random numbers in [1, 6] overall. For example, from a single rand call with RAND_MAX equal to 32767, we could view the value as a base-six numeral from 000000 to 411411. If it is under 400000, we can take the last five digits, which are each uniformly distributed in [0, 5], and adding one gts us the desired [1, 6]. If it is in [400000, 410000), we can use the last four digits. If it is in [410000, 411000), we can use the last three, and so on. Additionally, the otherwise discarded information, such as the leading digit, might be pooled over multiple rand calls to increase the average number of outputs we get per call to rand.

I'll also note that some `rand()` implementations are known to have very poor random distributions of the lower-order bits. A straight modulo operation would result in non-random results from such a `rand()` implementation. The bias introduced by a straight modulo implementation is probably completely insignificant in the results from a perfect `rand()` implementation as long as `RAND_MAX` is much greater than the divisor. — Andrew Henle, Oct 10 '19 at 12:43
Ah, splendid answer! I somehow assumed that RAND_MAX was divisable by 6.. This together with @chux comment clarifies it. In summary -> `9/6 > 0.5 -> 1 != 0`. Accepted! — GlassShark, Oct 10 '19 at 13:55
@GlassShark note that this probably isn't "the best way" in general, see http://www.pcg-random.org/posts/bounded-rands.html for some other ways of doing this and the various trade-offs — Sam Mason, Oct 11 '19 at 20:06

rand() example code, unnecessary check for larger than max?

1 Answers1

Footnote