Does this function result in uniformly distributed integers?

Question

Recently, on a project, I encountered the need for uniformly distributed integers within an arbitrary range [a, b] from random bytes, a problem that is quite common and is usually solved using rejection sampling. I was curious, however, if it would be possible to achieve this selection in constant O(1) time without the need for any branching. My thought process led me to the idea of using modular arithmetic.

Imagine a room with a clock that ticks through the range 0 to n - 1, where n is any positive integer that is >= 2. The clock ticks at a constant rate and never stops. Now, if you leave the room and come back after a random amount of time, chosen uniformly from the range 0 to m - 1 (where m >= n), and record the time shown on the clock each visit, you’ll be left with a list of random integers within the desired range.

I am not entirely sure why this works, or at least appears to, but empirical testing suggests that it produces level histograms with no modulo bias.

Both images were created with 10,000,000 random bytes fitted into the range 0 to 251.

bytes % 252

My method

Here is my implementation in Python. Please let me know if this has been considered before, and if there is any real mathematical backing to it.

Thanks.

import os
STATE = 0
def random_int(a: int, b: int):
 global STATE
 int_range = b - a + 1
 bytes_needed = (int_range.bit_length() + 7) // 8
 value = int.from_bytes(os.urandom(bytes_needed))
 STATE = (STATE + value) % int_range
 return STATE + a

• Luke Woodward: I think you have to be careful if the range of integers isn't fixed. Calling random_int with a small range will bias STATE to a lower value. I plotted the distribution of your function when called with moduli 5 and 252 randomly interspersed, and there was a clear bias towards the low end on the mod-252 plot. The effect is even more marked if you move the higher modulus further away from a power of 2.
• Marz: @LukeWoodward That's a very good point, actually. My implementation is not well-designed, so calling with alternating moduli will cause bias. The uniformity only arises with uninterrupted successive calls with a single modulus. If you were to use this in a real project, it would be necessary to have separate states for each modulus if you wanted to use the singular function in this way. Additionally, If the state has sufficient entropy, generating a sequence of ints in one modulus, then switching to another should not cause bias.
• Robert Dodier: Interesting problem, although the proposed solution seems more complex than it needs to be; I don't see that STATE is helping here. On the face of it, it seems like return a + value % int_range would be enough, and I think that without STATE in the picture, successive calls with different ranges would be OK. Also, although I think the analogy with the ticking clock is intriguing, there doesn't seem to be a clear relation to the proposed solution. Hope this helps in some way.
• Severin Pappadeux: @RobertDodier there is a paper by D.Lemire which is best method as of now. Put it into answer

Answer 1

This algorithm produces samples that are highly correlated with each other. While the cumulative distribution over time might be uniform, each sample is selected from a non-uniform distribution that depends on the previous sample.

Consider int_range = 192. For this, bytes_needed will be 1, and value will be, if os.urandom works ideally, uniformly distributed in [0, 256).

Consider what happens in (STATE + value) % int_range. When value is in [0, 192), STATE is cyclically advanced by an increment ranging uniformly from zero to the full cycle. When value is in [192, 255), STATE is cyclically advanced by an increment ranging uniformly from a full cycle (192) to a full cycle plus a third (256 − 192 = 64, one-third of 192).

Thus, over the value values in [0, 256), each residual advance in [0, 64) is produced twice, while each residual advance in [64, 192) is produced only once.

So, while this algorithm might produce a uniform distribution over time, its samples are not uncorrelated with each other. After a sample value of 11, the next sample is twice as likely to be 21 as it is to be 111.

It is simple to prove a pseudo-random number generator cannot generate its next sample with a uniform distribution over m values when it is using all the values in [0, B) for its selection unless m divides B: If it did, then each output value would come from an equal number, say t, of values in [0, B). If each of m values is mapped to from t values and that uses each value in [0, B) once and only once, then tm = B, so m is a factor of B. Since your m does not generally divide the B you are using, some power of 2, no algorithm can generate the next value with a uniform distribution except by not using all of the values in [0, B) (i.e., it must reject some values and generate a new state).

(One can, of course, make the deviation from uniformity arbitrarily small by making B sufficiently large.)