Dissecting Lemire's nearly divisionless random number generator

If my understanding is correct, this is an important distinction because it means you can plug different RNGs into this algorithm depending on your needs (in particular the 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.

---

    while (true) {

        u64 fract = random64();

        u128 upto_s = (u128)fract * (u128)s;
        u64  fract_part = (u64)upto_s;
        u64  int_part   = (u64)(upto_s >> 64);

    . X . . X . . X . . X . . X . .

    X . . X . . X . . X . . X . . X

        # Faster early test to avoid modulo
        # s > (2⁶⁴) % s
        if (fract_part >= s) {
            return int_part;
        }

        # 2⁶⁴ % s == (2⁶⁴ - s) % s == (-s) % s
        u64 min_fract = (-s) % s;
        if (fract_part >= min_fract) {
            return int_part;
        }
    }

https://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-...

Absolutely strongly disagree, especially given what's happened to the "average programmer" over time. The only way you can maintain (and improve) skill is by NOT continually pandering to the dumbing-down crowd.

In other words, what someone might call "more readable" can actually be "cruft that gets in the way" to someone else.

Related: https://www.linusakesson.net/programming/kernighans-lever/in...

I've "greyed" bits coming off of random radio waves from the atmosphere in order to test randomness; the drivers i was using for the particular radio didn't allow on-the-fly amplifier adjustments, so if there was, say, a lightning strike within 50 miles, the front-end would get overloaded and "blank" - send mostly zeros, even with incredibly fast AGC control. So entropy wasn't guaranteed unless i could guarantee there wasn't lightning within some predetermined distance of the receiver (as an aside, i have this capability now) - thus for long (96 hour) runs to gather entropy, i'd have to fiddle the output to ensure uniform distribution.

Now, with all that said, there was a side channel attack; if there was a lightning storm or other sort of high energy event near my receiver, an "attacker" could determine that the algorithm to ensure uniformity was contributing more to the bits than the atmospheric noise. Possibly.

As I recall another colleague wrote the original version (I do not remember who.) I bitbashed and optimized it, and added the comments, mainly because I wanted to make sure I understood the reasons it was right.

100% agree.

This is not related only to the importance of code reviews or team work, but directly to the way software is made, that is organic growth over time.

Software is too complex to be designed and built in one step, at least for us mere mortals.

For this reason we need many steps, where software has to be readable, compilable and testable.

...Influencer, Oxford comma enthusiast and Thought Leader Colm MacCárthaigh.

Yikes! One assumes then that he's a Thought leader when it comes to Oxford comma enthusiasts? This is doing my head in.

"Welcome to the very serious blog of Senior Principal Engineer, Keynote presenter, Touring Musician and producer, Cryptography Person, Open Source enthusiast, Philosopher, Influencer, Oxford comma enthusiast and Thought Leader Colm MacCárthaigh."

https://veryseriousblog.com/about

This suggests some questions that might be interesting.

1. For a given application, N is often fixed or at most has a small number of predetermined values. For what values of N can you find a hard coded [0, N) generator that is faster and still reasonably simple?

For example when N is a power of 2, then you can just call the underlying generator and then AND with the appropriate mask (possibly after shifting if the low bits of your underlying generator are not suitable).

If N is a power of 2 and not very large, and the low bits are good, you could take your underlying 64-bit number, split it by masking and shifting into multiple numbers in [0, N), return one of them and save the rest in an array. On subsequent calls return the numbers from the array until you are out.

If you were generating random bytes, say, this would cut your calls to the underlying random number generator by 87.5%.

For small N that are not powers of 2 this could still work. If N were 1000, say, you could treat the 64-bit random number as 6 numbers 10-bit number plus 4 leftover bits. If none of those 6 numbers is in [0, 1000), loop until you get a set of 6 that does have at least one in [0, 1000). (This will happen about once every few billion calls). Return the first one that is in [0, 1000), and store any others in an array to return on subsequent calls.

2. Can you write a random number generator generator? You call it with an N, and it returns a random number generator for [0, N) that is at least as fast as the fast general generator from the article?

I've looked at some special cases, but only for small N and small random number generators. Specifically, generating a dN with a dM, and aimed at something a human could do, by having a state diagram that a human could print out and follow as they roll the dice (or flip coins). For example, here is how to roll a d5 if all you have is a d8, and how to roll a d8 if all you have is a d5 [1]. Kind of neat, and great for a human, but wasteful for a computer, and if you try anything big the number of states grows rapidly. Even just rolling a d19 with a d32 has 325 states.

[1] https://imgur.com/a/XuDprvs

But we don’t generate 1-6 via division, we generate it by multiplying a random number between 0.0 and 1.0 times 6 and add a 1. I’m a bit fuzzy on this division comment. In mathematics division and multiplication are inverse operations, but in CPUs multiply has traditionally been much faster.

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

    return floor(s * rand(0, 1))

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

It is not suggested to use an rng with a known statistical bias, when a random number is required.