Dissecting Lemire's nearly divisionless random number generator
After reading the intro of the article, I wanted to see whether I could get the gist of the algorithm just from looking at the first code of it that I come upon. I thus took the first such link from the article (<a href="https:&#x2F;&#x2F;arxiv.org&#x2F;pdf&#x2F;1805.10941.pdf" rel="nofollow">https:&#x2F;&#x2F;arxiv.org&#x2F;pdf&#x2F;1805.10941.pdf</a>) while trying to avoid any spoilers. There was no way to unsee the comment at the top of the article, but I don&#x27;t think it helped much.<p>I did skim the paper to find which of the listed algorithms is the one in question (Algorithm 5), which inevitable primed me a bit (primarily by having the context of the two other rejection sampling algorithms). I also looked at the pseudo-code instead of the real code, so I don&#x27;t think this is a fair comparison with whatever colleagues of the author had to face, but I had no chance to make it rigorous anyway.<p>The algorithm in question, copied straight from the paper (and ASCII-fied):<p><pre><code>    ALGORITHM 5: An efficient algorithm to generate unbiased random integers in an interval.
    Require: source of uniformly-distributed random integers in [0, 2^L)
    Require: target interval [0, s) with s in [0,2^L)

    x := random integer in [0, 2^L)
    m := x * s
    l := m % 2^L
    if (l &lt; s) then &#x2F;&#x2F; Application of the rejection method
      t := (2^L−s) mod s &#x2F;&#x2F; 2^L mod s = (2^L − s) mod s
      while (l &lt; t) do
        x := random integer in [0, 2^L)
        m := x * s
        l := m mod 2^L
      end while
    end if
    return m &#x2F; 2^L
</code></pre>
Looking only at the part outside the if, I felt the overall idea quite intuitive - it&#x27;s essentially the equivalent of<p><pre><code>    return floor(s * rand(0, 1))
</code></pre>
except you&#x27;re working with L bit wide fixed point numbers.<p>The remaining question is &quot;when do you have to actually rejection sample&quot;. This part is not as obvious.<p>Thinking about the ranges of the variables involved, it becomes clear that m can take values in [0, s * 2^L) that are steps of s apart, and we want to bucket this space into segments of length 2^L. Some of these segments might contain one more of those possible m values than the others. What the `if` therefore has to do is &quot;shorten&quot; the segments so that each one covers the same number of possible values (at least that is my assumption for now - see below). The dead space in between is then where we have to actually do the &quot;rejection&quot; part of rejection sampling.<p>This seems to be the core of the idea - trim the fat <i>between</i> the segments, so that they can all end at a 2^L boundary, thereby eliminating the non-power-of-two division for figuring out which segment m fell in.<p>At this point I&#x27;m slightly stumped why the condition is different in the while. We compare with t instead of s, where t &lt; s, but is this required for equiprobability or just for performance reasons?<p>I am left to wonder whether the results for all &quot;not taking the if branch&quot; rolls are equiprobable.<p>If they are, why not rejection-sample on the if directly, completely avoiding any non-trivial remainder calculation in the first place?<p>If they are not equiprobable, that means the results of taking the if branch have to exactly undo the bias of the if-not-taken cases. I assume this is the case, since the goal was to eliminate as many &quot;hard&quot; remainder calculations as possible. On the other hand, it would be very odd for the probabilities involved to work out just right. I mean, if rejection sampling on the if alone has bias, and you only enter the if with a certain probability, and then need to correct that bias while influenced by the conditional probability - that just seems odd.<p>In other words: Does this algorithm have bias? If yes, it&#x27;s indeed freaky how the part within the if above compensates that.<p><pre><code>    x := random integer in [0, 2^L)
    m := x * s
    l := m % 2^L
    if (l &lt; s) then
      start over
    end if
    return m &#x2F; 2^L
</code></pre>
I&#x27;ll investigate this more deeply (probably by actually reading the full article and paper), but only after some sleep.