Returning to Perfect Randomness
About a year ago, I published a post about an algorithm for randomness in floating point. Recently, Sebastiano Vigna of the University of Milan reached out about his verification of my algorithm and some tests of my code because he was putting together a Rust crate for random floating-point numbers.
I am going to correct the previous post: Sebastiano found that one of my benchmarks (particularly the one that
looked good) had the wrong baseline because I had aliased the word rand accidentally, benchmarking my
algorithm against itself. The corrected version indicates that there is actually a performance cost when using
a PCG entropy source, aligned with the overhead given a ChaCha8 random source.
Sebastiano also found a very nice algorithm on Github by Pekka Pulkkinen that is faster than mine, using some of the same ideas but reorganizing them to have a quicker fast path. It is my opinion that random number generation in floating point should probably use a floating-point-aware algorithm in the future, and specifically Pulkkinen's algorithm, but it is not free compared to the typical fixed-point method.
Pulkkinen's Algorithm
Pulkkinen's algorithm is below:
1func Float64full() float64 {
2 // Fast path - full-precision division by a power of 2 via bit hacking
3 u := Rand64()
4 z := uint64(leadingZeros(u)) + 1
5 if z <= 12 { // 99.975% of cases
6 return floatFrombits((1023 - z) << 52 | u << z >> 12)
7 }
8 z--
9
10 // Slow path: zoom down until you hit a nonzero number
11 exp := uint64(0)
12 for u == 0 {
13 u = Rand64()
14 z = uint64(leadingZeros(u))
15 exp += 64
16 if exp + z >= 1074 { return 0 }
17 }
18
19 // Non-rounded division combined with a division by 2^exp
20 u = u << z | Rand64() >> (64 - z)
21 exp += z
22 if exp < 1022 {
23 return floatFrombits((1022 - exp) << 52 | u << 1 >> 12)
24 }
25
26 // Subnormal path: zero exponent and undo left-justification of u
27 return floatFrombits(u >> (exp - 1022) >> 12)
28}
Pulkkinen relies on a 64-bit integer source, as is common in random number generators. His fast path starts with a 64-bit
integer and then performs a division by $2^{64}$ without losing precision. The fast path uses a single if
statement that returns in almost all cases. My previous algorithm had a few branches to reach this point and more
micro-ops. There is a clever bit hack here: Pulkkinen uses the leading zero count of the integer to determine
the exponent, z, and then places the remaining bits of the integer into the mantissa (shifted appropriately by z).
His slow path re-draws until you get a non-zero number, dropping the exponent by 64 each time. This is optimized for 64-bit entropy sources, while my previous algorithm drew 53 bits at a time, matching the numerical precision of the floating-point format. Thus, Pulkkinen's algorithm enters its slow path with equivalent probability but loops with much lower probability. This stage tracks the running exponent to use, adding 64 for each loop, and includes a branch for the vanishingly few cases that return a zero. Once you escape the loop, there is the same full-width division, this time modified by the exponent, with a special handler for the subnormal numbers.
Next in Floating-Point Randomness
It has been on my list for the last year to work on correct random functions for exponential and Gaussian distributions, and I want to announce the challenge for those with more time than I. From my survey, these are the next most common random number generation functions behind uniform random floating-point numbers. I believe there is likely to be a very good algorithm for exponentially distributed numbers, as the exponent field of the result is a uniformly distributed integer (although it's not quite this simple). The mantissa field needs some careful handling: It may suffice to use a stochastically rounded exponential function to compute this distribution given a perfect uniform number generator, but this is currently just speculation.
Gaussian random numbers are very hard to generate, and I would not be surprised if a future literature survey reveals large errors in the probability distributions for current algorithms, but it will take some time to produce fully correct Gaussian random generator. The path from uniform randomness to a Gaussian usually transits through both branches and transcendental functions, and I am not so convinced that the simple application of stochastic rounding to these methods will work.
A Note On Humility
The post I made showed a very cool algorithm, and a presentation at FPTalks about this algorithm and this problem does seem to have inspired some discussion and further work, which I am happy about. However, at the time of the post, I wasn't as aware of the state of the art as I probably should have been, and the tone of the post was a bit too bombastic for the level of knowledge (or ignorance) of the subject I had at the time. Many people commented on this but could not find a faster algorithm. The benchmarking mistake was also humbling to see: I try to be very scrupulous about benchmarks for my posts, but the name aliasing mistake was something I didn't catch. In hindsight, this algorithm should not have been free, so a benchmark showing that fact was probably wrong.
Thanks to Sebastiano for catching my benchmarking mistakes and doing a much more thorough search. I am really excited for the crate and the amount of work that has gone into this. Please take a look at his crate if you are using random numbers in floating point!