3
$\begingroup$

Say I have a random number generator that generates a number within [0, RAND_MAX], and RAND_MAX < UINT_MAX.

How do I generate a random number within [0, i] such that i>RAND_MAX and i<UINT_MAX, while maintaining an even distribution, and without exceeding UINT_MAX in any calculations?

The first attempt is to split the range into n equal parts such that i mod n == 0 by finding the least denominator of i (using one of the integer factorization sieve algorithms). Each sublist has equal distribution, any can be picked at random. This fails if i is prime:

0 1 2 3 4 5 6
|___| |_____|

The second attempt is to split the range into two parts a = [0, floor(i/2)] and b = [floor(i/2)+1, i], and then to calculate the distribution as the ratio of the length of the two sub-ranges:

g = gcd(length(a), length(b))
if random(0, length(a)/g + length(b)/g - 1) >= length(a)/g
    random(0, floor(i/2))
else
    random(floor(i/2)+1, i)

Given the following range, the right sub-range is selected twice as often as the left sub-range:

0 1 2 3 4 5
|_| |_____|

Unfortunately the first call to random() may recurse forever if gcd() is ever one.

The final attempt is to correct the distribution skew introduced by splitting a range of odd length:

if random(1)
    if random(1)
        random(0, floor(i/2))
    else
        random(floor(i/2) + 1, i)
else
    if random(1)
        random(0, floor(i/2) - 1)
    else
        random(floor(i/2), i)

Which probably doesn't work.

Edit:

For anyone with similar question, the sum of two uniform variables follows a triangular distribution (see Irwin-Hall Distribution). This is invalid:

# triangular distribution:
random(a, floor(b/2)) + random(a, b - floor(b/2))

Edit (on accept):

The two answers I considered are both of high quality. This answer considers RNGs on a per-digit or per-bit basis. The accepted answers consider RNGs numerically, without considering digits or bits. In both cases, the comments are important to read.

I've also posted a follow-up question Maintain Uniform Distribution across Subranges.

$\endgroup$
2
$\begingroup$

For the sake of your example, let's say RAND_MAX=12 and i=17. Then do the following procedure: Choose two random numbers $r_1,r_2$ and combine them to a single random number uniformly distributed in $[0,144)$ by computing $r=r_1*12+r_2$. This is of course the wrong interval. You get a uniformly distributed random number in $[0,17)$ by repeating this process until you end up with a number $r$ in this interval; in other words, every time you get a random number $r\ge i$, you discard it and just try again. This guarantees that you end up with a uniform random number.

$\endgroup$
  • $\begingroup$ Thanks, I never specified RAND_MAX=b^n, where b is a base and n the number of digits, so r = r1 * RAND_MAX + r2 answers most of my question. You meant [0,144) instead of [0,24), and during the discard phase the interval can be [0,17*n) such that 17*n < r to reduce discards. However I also wish to avoid overflow, ie r<RAND_MAX<UINT_MAX. From the example (RAND_MAX=12, i=17), let UINT_MAX=31. Overflow can be detected, and discards should happen if r >= i or if r overflows. But given uniform distribution, r1*12 + r2 overflows more than 20% of the time. Can this be minimized? $\endgroup$ – user19087 Dec 23 '15 at 3:23
  • $\begingroup$ This might work - I'm not sure if it remains uniform: $r=RandomRangeInclusive(0,1)*12 + r_2$, to generate a random number in $[0,24)$. $\endgroup$ – user19087 Dec 23 '15 at 3:48
  • $\begingroup$ Oh, is that what you intended in the first place? $\endgroup$ – user19087 Dec 23 '15 at 3:50
  • $\begingroup$ You were right with the interval $[0,144)$. I've fixed my answer. $\endgroup$ – Wolfgang Bangerth Dec 24 '15 at 2:43
  • $\begingroup$ Regarding the overflow, the procedure shown above leads to a uniformly distributed number in $[0,144)$ but this leads to too much discard because $144\gg 17$.. What you can do is take a modulo operation on anything that properly divides 144 but is larger than 17. For example, if you compute $r = (r_1*12+r_2) mod 24$ then you get a uniformly distributed number $r$ in $[0,24)$ which has significantly less discard when you need to throw away everything larger than or equal to 17. Likewise, you could have done $r = (r_1*12+r_2) mod 18$ for even less discard. $\endgroup$ – Wolfgang Bangerth Dec 24 '15 at 2:47
2
$\begingroup$

Instead of fixing the provided generator use a reasonable modern choice. It will be faster and have better statistical quality. Possible examples include the various variants of xorshift+, xorshift* and PCG.

To directly respond to the question asked. You can generate a sample, mask out the maximal power of two bits available, generate another and shift that and bit or in the previous. For a non power of two 'n' perform a rejection method.

EDIT 2: So using a 64-bit xorshift+ to generate a 64-bit uniform sequence might look like this:

// multiple 'state' blocks to allow for friendly multi-threading
typedef struct {
  uint64_t s0;
  uint64_t s1;
} rng_state;

// get 64-bits
inline uint64_t rng_next(rng_state_t* s)
{
  uint64_t s1 = s->s0;
  uint64_t s0 = s->s1;

  s->s0 = s0;
  s1   ^= s1 << 23;
  s->s1 = s1 ^ s0 ^ (s1 >> 18) ^ (s0 >> 5);

  return s->s1 + s0; 
}

Trying to "patch-up" random, given simplifying assumptions that it returns a uniform an N-bit number and extending it to 2N-bits:

inline uint64_t rng_next() {
   uint64_t r0 = random();
   return (r0 << RANDOM_BITS) | random();
 }

Taking a quick peek at source it looks like the current version of (say) glibc will return 31-bits and by default it uses a power-of-two LCG (very low quality). The end result of patching up is significantly lower quality and significantly higher runtime cost.

NOTE: All code is typed in post and probably doesn't even compile.

$\endgroup$
  • $\begingroup$ Unless you name a reasonable modern choice and make a convincing case that implementing this generator is simpler than answering the question, I'm not considering this. Also keep in mind that random() as specified by the POSIX.1-2001 specification is relatively modern, especially platform-specific implementations. Finally, every PRNG has an equivalent to RAND_MAX, so my question still stands. And besides, this is an interesting question - how can the distribution be preserved across an uneven split? $\endgroup$ – user19087 Dec 21 '15 at 21:54
  • $\begingroup$ You can generate a sample $\endgroup$ – MB Reynolds Dec 21 '15 at 21:57
  • $\begingroup$ Trouble typing on cell. Updated my answer. $\endgroup$ – MB Reynolds Dec 21 '15 at 22:08
  • $\begingroup$ Thanks, this is what I'm looking for, but I'm not familiar with the terminology: "generate a sample" -> call random() ? "mask out the maximal power of two bits available" -> every set bit is a power of two, so 0b0 ? "shift that and bit or in the previous" -> not sure what you mean. Re: rejection, I didn't follow the previous steps so I'm not sure what you're getting at. $\endgroup$ – user19087 Dec 22 '15 at 0:56
  • $\begingroup$ Do you mean, something like stats.stackexchange.com/a/31133 ? $\endgroup$ – user19087 Dec 22 '15 at 1:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.