Fast Random Integer Generation in an Interval

capabilities of the underlying platform). Related question: For the ultimate in minimal state, would a CBPRNG such as Philox or Threefry (http://www.thesalmons.org/john/random123/papers/random123sc1...) be safe to use here or would the fact that it can be invoked multiple times for a single call (and thus state sequences might end up overlapping between calls) be likely to introduce subtle statistical issues? Slightly off topic - from Lemire's paper: > may not be applicable to specialized processors such as Graphics Processing Units (GPUs) that lack support for the computation of the full multiplication This statement seems to be at odds with the fact that implementations of Philox are provided for both AMD and Nvidia GPUs but Philox relies on efficient mulhi and mullo being available. I didn't bother to look into it further yet though.

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.