Add functions to generate integers pseudo-randomly

[xavierleroy]

This PR adds functions to generate pseudo-random integers:

• Z.random_int N returns an integer between 0 (included) and N (excluded)
• Z.random_bits n returns an integer that fits in n bits, i.e. between 0 (included) and 2^n (excluded).

The source of randomness is the Random standard library module, which is not cryptographically secure. For more demanding applications, lower-level functions Z.random_int_gen and Z.random_bits_gen are provided, parameterized by a user-provided function that draws random bytes into a byte array. The documentation explains how these lower-level functions can be used to produce cryptographically-strong random integers.

Addresses: #99

[xavierleroy]

I'm keen on having this new feature in Zarith. Could this PR be reviewed, please? Thanks.

[antoinemine]

Thanks for the PR and sorry for the delay.

This seems very nice! Although, I am not a number theory specialist and I cannot judge whether the way you draw uniform random numbers for non-power-of-two bounds is the optimal way.

Is there a rational behind the 4 in f (nbits1 + 4)? I guess that larger numbers decrease the risk of reject, at the cost of drawing larger random numbers in the first place.

[xavierleroy]

I am not a number theory specialist and I cannot judge whether the way you draw uniform random numbers for non-power-of-two bounds is the optimal way.

I'm not either :-) The algorithm is essentially rejection sampling, with the obvious twist to reject less often. The same algorithm has been used for ages in the Random module of the OCaml standard library.

Is there a rational behind the 4 in f (nbits1 + 4)? I guess that larger numbers decrease the risk of reject, at the cost of drawing larger random numbers in the first place.

Yes, exactly. In the Random module, we draw fixed numbers of bits, and the worst case (Random.int 0x20000001 if I'm not mistaken) is that we reject with probability almost 1/2, so we need two draws on an average. With the + 4, the rejection probability drops to 1/32 (I think), so we need 1 + epsilon draw on an average, but the extra 4 bits can cause an extra call to Random.bits or an extra byte to be read from /dev/random. I think it's a good compromise. Of course, + 4 could be some other small number, I didn't try to measure the difference it would make.

[pascal-cuoq]

To the best of my knowledge, the state of the art for generating random numbers in an interval is https://arxiv.org/pdf/1805.10941.pdf . It avoids the systematic division by replacing it by a systematic multiplication. (This is not a suggestion to change anything in the MR.)

[xavierleroy]

Interesting! Thanks for the reference. Maybe this PR can be merged with the naive division-based algorithm for the time being, but I'd be happy to revisit it later with the multiplication-based algorithm, given enough time to look into it.

[antoinemine]

OK for me! Unless there are any objections, I'll merge.

[xavierleroy]

No objections :-) I just squashed all commits in one.