Some context on the Lemire idea if others haven't seen it:

https://lemire.me/blog/2019/06/06/nearly-divisionless-random-integer-generation-on-various-systems/

https://shufflesharding.com/posts/dissecting-lemire

https://sts10.github.io/2020/10/10/lemire-neaarly-divisionless-random.html

For bias validation -- you could do some kind of primitive statistical test (histogram). It can't be perfectly reliable but could give a little confidence it isn't totally broken.

The C code in Lemire's paper uses a uint64_t bound and __uint128_t variable. Maybe it is worth a small note to this effect in the commit message (of course, 32/64 or 64/128 doesn't matter for the algorithm itself).

From @cem's comment I think the point is that the bias can be easily observed with a simple approach like a histogram of results by choosing an appropriate upper bound. For instance, choose an upper bound of 2³²✕⅔; with the naive arc4random() % 2863311531, 1000000 iterations and a 16-bin histogram I get:

0: 125241
1: 124677
2: 124899
3: 125077
4: 125090
5: 83651
6: 62995
7: 62451
8: 61890
9: 62520
10: 41509
11: 0
12: 0
13: 0
14: 0
15: 0

The bias is easily visible. When I perform the experiment with arc4random_uniform(2863311531) the histogram looks as we'd expect:

0: 94402
1: 93591
2: 93493
3: 93299
4: 93506
5: 94024
6: 94089
7: 94252
8: 93504
9: 93060
10: 62780
11: 0
12: 0
13: 0
14: 0
15: 0

But note that bias is apparent when I try the proposed change here -- I get:

0: 125090
1: 124805
2: 124952
3: 124626
4: 124852
5: 83803
6: 62237
7: 62913
8: 62071
9: 62763
10: 41888
11: 0
12: 0
13: 0
14: 0
15: 0

Yeah, you can do a simple bias test with an upper limit of 3 for example. Bucket 0 is more likely than buckets 1 and 2 under a simplistic rand32() % 3 implementation. Or something like that. It is probably easier to measure with a smaller power of 2 generator.

In D47659#1087607, @cem wrote:

Yeah, you can do a simple bias test with an upper limit of 3 for example. Bucket 0 is more likely than buckets 1 and 2 under a simplistic rand32() % 3 implementation. Or something like that. It is probably easier to measure with a smaller power of 2 generator.

Quite on the contrary; with the biased modulo reduction and a limit of 3, 1431655765 numbers each become 1 or 2, but 1431655766 numbers become 0. The bias is very tiny and requires many billion observations to eek out.

IIRC the best way would be to use a threshold around (2/3)*UINT32_MAX, though even that still requires at least 5 observations per bucket for statistical significance of the chi-square test, requiring many billion observations in total. This does not seem suitable for a unit test that is part of a larger test suite.

Yes, yes. It is easier to analyze with e.g. rand2() % 3. But you and Ed have both mentioned a way to measure this for rand32 more tractably -- great!

I am not suggesting a bias test has to be a unit test.

I'll see that I can write a bias test outside of our test suite. It'll probably take quite a while to run.

Here is a draft bias test for arc4random_uniform(). It samples the distribution of arc4random_uniform() but doesn't yet have code to do a chi-square test or similar.

https://gist.github.com/clausecker/1a4e0d0358f26ef5a6fd01d74b187799

Using this test, we get for a biased implementation arc4random_uniform(uint32_t threshold) { return arc4random() % theshold; }:

distribution properties:
        threshold:      3221225472
        observations:   21474836480
        least common:   6719 (0 observations)
        most common:    83532884 (35 observations)
        median:         1073740653 (expected 1610612736)
        average:        1342174728.876399 (expected 1610612735.500000)
        deviation:      942708369.486011 (expected 929887696.689840)

sample properties:
        median:         6
        average:        6.666667
        deviation:      3.496032 (expected 2.581989)

for Lemire's algorithm as proposed in this DR:

distribution properties:
	threshold:	3221225472
	observations:	21474836480
	least common:	353 (0 observations)
	most common:	647634042 (29 observations)
	median:		1610608903 (expected 1610612736)
	average:	1610614338.364706 (expected 1610612735.500000)
	deviation:	929886996.823639 (expected 929887696.689840)

sample properties:
	median:		6
	average:	6.666667
	deviation:	2.582016 (expected 2.581989)

and for the currently shipped implementation from OpenBSD:

distribution properties:
	threshold:	3221225472
	observations:	21474836480
	least common:	1683 (0 observations)
	most common:	1512057592 (28 observations)
	median:		1610621068 (expected 1610612736)
	average:	1610616491.083907 (expected 1610612735.500000)
	deviation:	929889833.268859 (expected 929887696.689840)

sample properties:
	median:		6
	average:	6.666667
	deviation:	2.581999 (expected 2.581989)

As you can see, the numbers are pretty much dead on for the implementation we ship and the one from OpenBSD, but have a high discrepancy for a known biased implementation. I hope that shows that the code yields unbiased results.

@cem Do you think this statistical test instils sufficient confidence to commit this change set?

I have no objection.