Hacker Newsnew | past | comments | ask | show | jobs | submitlogin

That results in many rejections when the inclusive limit is a power of two (or slightly larger). For instance, for [0, 4], one has to reject 38% of the time, and this gets closer to 50% as the limit increases. This requires doing more samples from the (P)RNG.

The "Faster Threshold-Based Discarding" section of that article is the technique that Lemire actually suggests, to avoid the %, and the benchmarks indicate Lemire's method is equal or faster in all but the "Large Shuffle" case (and the "Optimizing Modulo" section reduces that difference significantly).



One way to pick an integer in [0, 4] given a source of uniform random integers in [0, 7] (which I shall call "the d8") is to divide the interval [0, 1] into 5 equal parts, by splitting it at 1/5, 2/5, 3/5, and 4/5. Number the 5 intervals 0 through 4.

Write those four split positions in base 8:

  1/5 = 0.1463(1463)
  2/5 = 0.3146(3146)
  3/5 = 0.4631(4631)
  4/5 = 0.6314(6314)
The parenthesis mark groups of digits that repeat infinitely.

Use the d8 to generate an octal fraction in [0, 1]. Find which of the 5 intervals it falls in, and output that interval's number.

For example, if the d8 comes up 0, 2, 5, or 7, you would output 0, 1, 3, or 4, respectively. If you get anything else, you will need at least one more roll. For example, suppose the first roll was 1. Then you know that the output is going to be 0 or 1. Roll again. On 0, 1, 2, or 3, output 0. On 5, 6, or 7, output a 1. On 4, you need to roll a third time. On the third roll, 1-5 => 0, 7 => 1, and 6 means you'll need a 4th roll. The 4th roll goes 0-2 => 0, 4-7 => 1, and 3 means roll again. The pattern then repeats for potential rolls 5, 6, 7, and 8, and so on.

Another way to work out the same thing, without explicit fractions, is to work out a state diagram, naming each state with a string that lists each reachable output, with some outputs repeated so that they appear in the name in proportion to their probability from that state. That's probably confusion, so here is an example, where we will use a d2 (i.e., a random bit stream) to generate integers in [0, 2]. In other words, we are using a d2 to simulate a d3.

The first state should be able to lead to 0, 1, or 2, an they should be equally probable, so we name the first state 012. There should be two next states from that, one for rolling 0 and one for rolling 1.

The way we get the names for those states is by splitting 012 in two. But we need a name with an even number of characters for that, so lets double all the characters: 001122. Now we can split it in two to get the two next states: 001 and 122.

Same thing to get the next states for 001 (make it even--000011--and then split it: 000 and 011) and 122 (=> 112222 => 112 and 222). So far we have this (draw root down):

  000 011  112 222
    \  /     \  /
     001     122
        \    /
          012
That says if we roll 00 we output 0, and 11 outputs 2, and 01 and 10 need at least another roll. If you work out the next states from 011 you get 001 and 111. 001 was 011's parent, so we've got a loop there. Same on the other side. Putting those in gives:

     A 111 111 B
      \ /   \ /
  000 011  112 222
    \  /     \  /
    A:001    B:122
        \    /
          012
To roll a d3 with a d2, start at the root, use the d2 to traverse the tree, and when you reach a node whose name is 000, 111, or 222, output 0, 1, or 2, respectively.

This all easily generalizes to using a dN to simulate a dM.

The base-M fraction approach can be done without actually building up the table of partition points, and the digits of the one partition point you need for each simulated dM roll can be worked out on demand as you need them, using just integers.

For example, if you need a/b in base M, set r = a. To get successive digits, iterate this:

  r = r * M
  d = floor(r/b)
  r = r % b
The sequence of d that generates is the digits of the base M expansion of a/b.


I'm not entirely sure how that relates to my comment...

However, I think that approach is conceptually very similar to "FP Multiply" (using floating-point arithmetic) or "Integer Multiplication" (using fixed-point arithmetic) in the GP's linked blog post. Those approaches package the whole "generate digits" idea into a single step, by handling the "split positions" directly as native fractions, rather than having to think of them as sequences in a particular base. These are likely to be significantly more efficient too, since they do not have to do (multiple) division/modulos, which are expensive as the main article points out (let alone computing/storing the fraction representation).

Your approach has the benefit that it is not biased (if the entropy in `r` is tracked appropriately, so that it can be resampled), but Lemire's "Debiased Integer Multiplication" version (as well as the optimised versions) are likely a faster way to remove the bias.




Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact

Search: