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

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 (https://arxiv.org/pdf/1805.10941.pdf https://arxiv.org/pdf/1805.10941.pdf) while trying to avoid any spoilers. There was no way to unsee the comment at the top of the article, but I don't think it helped much. 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'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. The algorithm in question, copied straight from the paper (and ASCII-fied): 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 < s) then // Application of the rejection method t := (2^L−s) mod s // 2^L mod s = (2^L − s) mod s while (l < t) do x := random integer in [0, 2^L) m := x * s l := m mod 2^L end while end if return m / 2^L Looking only at the part outside the if, I felt the overall idea quite intuitive - it's essentially the equivalent of return floor(s * rand(0, 1)) except you're working with L bit wide fixed point numbers. The remaining question is "when do you have to actually rejection sample". This part is not as obvious. 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 "shorten" 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 "rejection" part of rejection sampling. This seems to be the core of the idea - trim the fat between 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. At this point I'm slightly stumped why the condition is different in the while. We compare with t instead of s, where t < s, but is this required for equiprobability or just for performance reasons? I am left to wonder whether the results for all "not taking the if branch" rolls are equiprobable. If they are, why not rejection-sample on the if directly, completely avoiding any non-trivial remainder calculation in the first place? 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 "hard" 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. In other words: Does this algorithm have bias? If yes, it's indeed freaky how the part within the if above compensates that. x := random integer in [0, 2^L) m := x * s l := m % 2^L if (l < s) then start over end if return m / 2^L I'll investigate this more deeply (probably by actually reading the full article and paper), but only after some sleep.