[infogulch]

> In 2018, Daniel Lemire found an algorithm that avoids the divisions nearly all the time (see also his 2019 blog post). In math/rand, adopting Lemire’s algorithm would make Intn(1000) 20-30% faster...

I recently found a super simple algorithm that appears to produce a number in the interval [0,N] with a branchless expression with a single multiplication in an extended number size. (Sorry I don't have a reference.)

Say you want to generate a number, G, in interval [0,N] where N<=UInt32Max. The algorithm is:

    G = uint32( uint64(N)*uint64(rand.UInt32())>>32 )

It seems like this should select a number in the range with no bias. Is there something I missed?

[ironhaven]

You have written a deterministic function. If you test this function with all 4 billion uint32 on one odd interval, go and count the number of times you get each result. Now look at your results, are all of the numbers equally likely or is there bias towards some outputs?

Ps: it looks like your function is exclusive like [0,N) not [0,N] Also your function is described in this blog post https://www.pcg-random.org/posts/bounded-rands.html

[infogulch]

Correct on all counts. Thanks!

[Someone]

> It seems like this should select a number in the range with no bias. Is there something I missed?

Yes. There are many values of N that aren’t divisors of UInt32Max.

As the article says: “However, no algorithm can convert 2⁶³ equally likely values into n equally likely values unless 2⁶³ is a multiple of n: otherwise some outputs will necessarily happen more often than others. (As a simpler example, try converting 4 equally likely values into 3.)”

[infogulch]

Are you sure? Maybe this isn't a good test but it seems pretty evenly distributed to me:

https://go.dev/play/p/IeJQEAclBCU

Edit: maybe this shows the bias better: https://go.dev/play/p/3eKJibIlF1a

[cmeacham98]

UInt32Max (i.e. 4294967295) is divisible by 3, so your code actually is perfectly random (or more accurately, as random as go's rand package). It would be biased with N=4, for example.

Regardless, with small values of N, the bias is very slight so you would need many many iterations to see the imperfection in a statically significant way.

[infogulch]

That makes sense.

A quick search didn't reveal any good resources for how to test the quality of a random number generator in a number range. Is what I came up with the best strategy, and you just need to run it for much longer (and compare to a known-good implementation) to see the difference?

[codeflo]

Why would you need to “see” it? Unlike the distribution of the RNG itself, this is trivial to solve analytically.

[josephg]

> (As a simpler example, try converting 4 equally likely values into 3.)

No, but you can convert a RNG that emits 4 equally likely values into an RNG that emits 3 equally likely values. Just - anytime the RNG returns 4, try again.

Here's a fun puzzle / annoying interview question: You have a biased coin. You can flip it as often as you want, but heads and tails are not equally likely. Without figuring out the bias of the coin, how do you produce purely random bits?

[stouset]

My guess after a few glasses of wine and thinking on it for a few minutes:

Flip twice. If both flips are the same discard the result. Output 0 for TH, 1 for HT.

[josephg]

Nice - you got it!

The followup is this: That approach only uses about 50% of your coin flips. The other 50% are discarded. How would you improve the efficiency?

[stouset]

Enlist my friends to flip more coins in parallel :)

At a high level I’d probably try and exploit the fact that every bit sequence with a given number of H and T has equal probability. e.g., HHHT HHTH HTHH THHH are equally probable and so can be mapped to four different values. That still only gets me 2 bits (50%) but other combinations (e.g., variations on HHTT) could get me log2(6) bits. I’m guessing with a higher number of flips I could extract (on average) more and more as a proportion. No clue what the asymptote would be.

[stouset]

Thinking further, for N flips you get 0 bits of entropy for all H or all T. For all other sequences, you get log2(N choose count(H)) bits of entropy, and you can average the sum of these over N.

According to Wolfram Alpha this works as N gets larger but it’s not great. For 16 flips you get 9.5 bits of entropy, but hey at least I beat half! 32 flips gets you about 20 bits of entropy. By 64 flips you get 43 bits, and that’s approaching 2/3 efficiency. Maybe not so bad after all!

Going higher is a little tough since I’m on mobile but it starts crawling reaching only 71% efficiency by 1024 flips. I’m curious if it does actually asymptotically reach 100% efficiency (for a fair coin), even if quite slowly.

Edit: Playing more[1] it really seems to approach 72.1. I wonder if I can figure out the asymptote analytically…

[1] https://www.wolframalpha.com/input?i=sum+from+i=1+to+N-1+of+...

[assbuttbuttass]

This maps UInt32Max input values to N output values, so there is guaranteed to be bias by pigeonhole, unless N divides Uint32Max

[fanf2]

Your algorithm is the first step of Lemire’s algorithm, without the followup check to debias the result. https://dotat.at/@/2020-10-29-nearly-divisionless-random-num...

[flyingmutant]

This algorithm produces biased result with probability 1/2^(32-bitwidth(N)). Using 64 or 128 random bits can make the bias practically undetectable. Comprehensive overview of the approach can be found here: https://github.com/apple/swift/pull/39143

[baruz]

Your results will be biased. It is tiny with small values of N, and absent when N is a power of two, but the skew becomes more obvious when your N is 2^31 + 2^30 + 1, for example.

[funny_falcon]

Lemire uses your found but corrects its biases.