4 ms·

Higher quality random floats

116 points by ingve 3y ago

raymondh 3y ago
FWIW, there is a recipe for this in the Python docs[0]. It was developed in collaboration with Alan Downey and Tim Peters: from random import Random from math import ldexp class FullRandom(Random): def random(self): mantissa = 0x10_0000_0000_0000 | self.getrandbits(52) exponent = -53 x = 0 while not x: x = self.getrandbits(32) exponent += x.bit_length() - 32 return ldexp(mantissa, exponent) [0] https://docs.python.org/3/library/random.html#recipes https://docs.python.org/3/library/random.html#recipes
FabHK 3y ago
Very, very nice. Julia currently uses the same trick as Lua (generate uniform Float64 bit pattern between 1 and 2 by masking a random UInt64 and reinterpreting, then subtract 1). I wonder what the speed is. Having the “other” half open unit interval in addition to the traditional one might come in handy for several applications.
- adgjlsfhk1 3y ago
  For rngs that natively generate 64 bits (some only give you 52 bits at a time), the speed is really good.
orlp 3y ago
Readers may also be interested in this old comment of mine: https://news.ycombinator.com/item?id=30895669 https://news.ycombinator.com/item?id=30895669. It uses the method by Downey briefly mentioned at the end of the article.
- toxik 3y ago
  Hm, this method seems odd to me. Unbiased yes but then rounded to nearest float. Consider each floating point value a bin. Then this will fill bins where the coverage is sparse faster. Eg say the possible floats are 1, 10, 100 and so on - you generate from 0 to 10000. There will be far more in bin 1000 than 100.
leiroigh 3y ago
This problem is much more acute in Float32. Another sample implementation (maybe easier to read?) is in https://discourse.julialang.org/t/output-distribution-of-rand/12192 https://discourse.julialang.org/t/output-distribution-of-ran... As far as I remember, the main reasons for not rocking the boat on that in julia were: Everybody generates rand floats wrong, so doing it right breaks expectations and compat; and there is a perf price of maybe 0.5-1.5 cycles/number to pay; and this is nontrivial to SIMD (basically because we don't have simd TZCNT / LZCNT / access to exponential distribution)
- wodenokoto 3y ago
  What kind of things are relying on random floats being generated wrong?
  - leiroigh 3y ago
    Comparability between implementations: Say GPU vs CPU, or between languages, or to pseudocode in old papers. Typical example where the badness of the floats bites you is if you do something like log(rand()), or 1/x, or more generally you map your uniform (0,1] interval via a cdf to generate a different distribution, and your mapping has a singularity at zero (this is like super standard -- e.g. how you generate exponentially distributed numbers. Or if you generate random vectors in an N-d ball, you're using a singular cdf to compute the magnitude your vector, multiplied with a normalized normal distributed N-d variable to generate the angular component). Then the bad random floats (i.e. the non-smooth distribution close to zero) introduce real and visible bias after transformation. Afaiu the problem is well-known and serious people correct for it somewhere. If you fix the underlying problem without revisiting the now-obsolete fixes, then your results are biased again. I'm not arguing for keeping the bad random float generation everywhere. I think it should be fixed, not just in julia but everywhere. I'm just saying that it's not a no-brainer and there is discussion and compat and communication and code audit involved in doing such a momentuous change. Also, I'm not really up to putting in the work to champion that at the moment (arguing on the internet doesn't count ;)).
- FabHK 3y ago
  Nice. I am not surprised that the approach presented here was discussed in Julia circles already some 5 years ago. So it was about half as fast, if I interpret the thread correctly.
ykonstant 3y ago
Lovely article and a very nice C implementation.
GerryB 3y ago
One could argue that the most important criteria to optimize for is insuring that the mean expected value is exactly 0.5. That is a more practical criteria than trying to allow super low values like 2^-70 which are indeed expressible as floating point, but just less important than maintaining the simple behavior of a proper average value.
- jksk61 3y ago
  one could argue that the most important criteria is that if you take any interval between 0 and 1 you get the same density... at least that would be a uniform distribution, since even a normal, exponential or whatever could have mean 0.5 and simulate that instead.
- eigenket 3y ago
  Optimizing for that criteria you obtain a "random number generator" which deterministically returns 0.5 every time.
- lifthrasiir 3y ago
  The proposed generator has, assuming the true uniform integer generator, the mean expected value of 0.5 + 2^-152 [1]. If you can ever distinguish this difference you are more likely to see a biased RNG output much earlier. [1] Start with the first `rand_between_zero_and_one` snippet. `x = ((x + 1) >> 1) + (e << 52)` can be rewritten as `d = (1.0 + ((x + 1) >> 1) * 2^-52) * 2^(e - 1023)` (since it always generates a normal number). `E[(x + 1) >> 1]` exactly equals to 2^-51, and `E[2^e] = 2^1022 * 2^-1 + 2^(1022-1) * 2^-2 + ... + 2^(1022-74) * 2^-75 + 2^(1022-75) * 2^-75 = 2^1023 (1/4^1 + 1/4^2 + ... + 1/4^75 + 0.5/4^75) = 2^1023 (1/3 + 1/(6*4^75))`. So `E[d] = (1 + 2^-51 * 2^52) * (1/3 + 1/(6*4^75)) = 1/2 + 1/4^76`.
Sharlin 3y ago
Often the naive method is all right as you aren’t really concerned by the fact that numbers close to zero don’t have the maximum precision possible – particularly when your use case has a certain "native scale" and what amounts to a fixed point approach works fine. And especially if the results are going to be discretized at some point anyway (eg. computer graphics).
camel-cdr 3y ago
The holy grail of random float generation is a function that takes two bounds and can generate ALL representable floating point numbers within the given bounds with the proper distribution, here uniform (and does so efficiently). I've written such an implementation with explanations [0], but it doesn't handle subnormals correctly. I've yet to find the time to revisit the problem. [0] https://github.com/camel-cdr/cauldron/blob/main/cauldron/random.h#L1462 https://github.com/camel-cdr/cauldron/blob/main/cauldron/ran...
- NeutralForest 3y ago
  what does subnormal mean?
  - camel-cdr 3y ago
    Usually floating point numbers are defined as `1.fraction * 2^(exponent-bias)`, but if the exponent is zero `0.fraction * 2^(1-bias)` is used, to allow representing very small numbers. IIRC I also explained it here: https://m.youtube.com/watch?v=VHJUlRiRDCY&t=2774s https://m.youtube.com/watch?v=VHJUlRiRDCY&t=2774s
    - NeutralForest 3y ago
      thanks!
    - Aardwolf 3y ago
      And the reason subnormals exist at all is that normally there's an implicit 1-bit (without the implicit 1-bit, there would be multiple representations of the same value, e.g. all values with mantissa 0 would represent 0 regardless of exponent), but you could never represent the value 0 itself if the implicit 1-bit is always there. So for the lowest exponent value, the implicit bit is set to 0, to allow representing 0, but also those other very tiny values since there's still an entire mantissa to give meaning to in addition to value 0.
      - magicalhippo 3y ago
jkafjanvnfaf 3y ago
After thinking about the problem a bit, I've come up with a completely different approach. I'm sure this is known, but I can't find anything so I'd appreciate if someone can point me to a citation. Generate two independent doubles a and b the "usual way", meaning from [0, 1) with equal spacings of 2^-53. Then calculate 2^-53*b + a. Assuming for simplicity that the rounding mode is towards zero, this gives a new random double in [0, 1). I've only done the math in my head so I'm not 100% sure on all the numbers, but I believe it scores pretty well at the criteria in the blog post; tentative results are Probability of zero: 2^-106 Smallest nonzero: 2^-106 Largest non-seeable: approx. 2^-54 Distinct values: 2^58.7 The usual rounding mode (round to nearest, ties to even) will generate numbers in [0, 1] with very slight biases that need closer analysis. One could also go further and do e.g. 2^-106*c + 2^-53*b + a, but something like that should only be needed for for single-precision floats. The obvious downside is that you always need two random integers and not just in the slow path. However, SIMD becomes easier because you don't need AVX512 for lzcnt.
10000truths 3y ago
This is exactly how random floats are computed in the Zig standard library: https://github.com/ziglang/zig/pull/10428/ https://github.com/ziglang/zig/pull/10428/
- leiroigh 3y ago
  This is awesome and the first precedent I've ever seen for a standard library doing the right thing on rand floats. Big kudos to the zig people and thanks for brightening my day!
- corsix 3y ago
  Oh, cute. It looks like they ever-so-slightly overweight the probability of values whose mantissa is entirely zeroes though. For example, the probability of hitting exactly 0.5 should be 2^-54 + 2^-55, whereas zig looks to give it 2^-54 + 2^-54.
kwillets 3y ago
The final methods seem equivalent to reading a non-float binary fractional number .xyz... from a random bitstream and stopping at the point where further digits will be truncated in the float representation. eg if x = 1 read 53 more mantissa bits and stop, if x = 0 decrement the exponent and repeat. There's no bias in such a process.

116 points by ingve 3y ago

FWIW, there is a recipe for this in the Python docs[0]. It was developed in collaboration with Alan Downey and Tim Peters: from random import Random from math import ldexp class FullRandom(Random): def random(self): mantissa = 0x10_0000_0000_0000 | self.getrandbits(52) exponent = -53 x = 0 while not x: x = self.getrandbits(32) exponent += x.bit_length() - 32 return ldexp(mantissa, exponent) [0] https://docs.python.org/3/library/random.html#recipes https://docs.python.org/3/library/random.html#recipes

FabHK 3y ago

Very, very nice. Julia currently uses the same trick as Lua (generate uniform Float64 bit pattern between 1 and 2 by masking a random UInt64 and reinterpreting, then subtract 1). I wonder what the speed is. Having the “other” half open unit interval in addition to the traditional one might come in handy for several applications.

adgjlsfhk1 3y ago
For rngs that natively generate 64 bits (some only give you 52 bits at a time), the speed is really good.

orlp 3y ago

Readers may also be interested in this old comment of mine: https://news.ycombinator.com/item?id=30895669 https://news.ycombinator.com/item?id=30895669. It uses the method by Downey briefly mentioned at the end of the article.

toxik 3y ago
Hm, this method seems odd to me. Unbiased yes but then rounded to nearest float. Consider each floating point value a bin. Then this will fill bins where the coverage is sparse faster. Eg say the possible floats are 1, 10, 100 and so on - you generate from 0 to 10000. There will be far more in bin 1000 than 100.

leiroigh 3y ago

This problem is much more acute in Float32. Another sample implementation (maybe easier to read?) is in https://discourse.julialang.org/t/output-distribution-of-rand/12192 https://discourse.julialang.org/t/output-distribution-of-ran... As far as I remember, the main reasons for not rocking the boat on that in julia were: Everybody generates rand floats wrong, so doing it right breaks expectations and compat; and there is a perf price of maybe 0.5-1.5 cycles/number to pay; and this is nontrivial to SIMD (basically because we don't have simd TZCNT / LZCNT / access to exponential distribution)

wodenokoto 3y ago
What kind of things are relying on random floats being generated wrong?
- leiroigh 3y ago
  Comparability between implementations: Say GPU vs CPU, or between languages, or to pseudocode in old papers. Typical example where the badness of the floats bites you is if you do something like log(rand()), or 1/x, or more generally you map your uniform (0,1] interval via a cdf to generate a different distribution, and your mapping has a singularity at zero (this is like super standard -- e.g. how you generate exponentially distributed numbers. Or if you generate random vectors in an N-d ball, you're using a singular cdf to compute the magnitude your vector, multiplied with a normalized normal distributed N-d variable to generate the angular component). Then the bad random floats (i.e. the non-smooth distribution close to zero) introduce real and visible bias after transformation. Afaiu the problem is well-known and serious people correct for it somewhere. If you fix the underlying problem without revisiting the now-obsolete fixes, then your results are biased again. I'm not arguing for keeping the bad random float generation everywhere. I think it should be fixed, not just in julia but everywhere. I'm just saying that it's not a no-brainer and there is discussion and compat and communication and code audit involved in doing such a momentuous change. Also, I'm not really up to putting in the work to champion that at the moment (arguing on the internet doesn't count ;)).
FabHK 3y ago
Nice. I am not surprised that the approach presented here was discussed in Julia circles already some 5 years ago. So it was about half as fast, if I interpret the thread correctly.

ykonstant 3y ago

Lovely article and a very nice C implementation.

GerryB 3y ago

One could argue that the most important criteria to optimize for is insuring that the mean expected value is exactly 0.5. That is a more practical criteria than trying to allow super low values like 2^-70 which are indeed expressible as floating point, but just less important than maintaining the simple behavior of a proper average value.

jksk61 3y ago
one could argue that the most important criteria is that if you take any interval between 0 and 1 you get the same density... at least that would be a uniform distribution, since even a normal, exponential or whatever could have mean 0.5 and simulate that instead.
eigenket 3y ago
Optimizing for that criteria you obtain a "random number generator" which deterministically returns 0.5 every time.
lifthrasiir 3y ago
The proposed generator has, assuming the true uniform integer generator, the mean expected value of 0.5 + 2^-152 [1]. If you can ever distinguish this difference you are more likely to see a biased RNG output much earlier. [1] Start with the first `rand_between_zero_and_one` snippet. `x = ((x + 1) >> 1) + (e << 52)` can be rewritten as `d = (1.0 + ((x + 1) >> 1) * 2^-52) * 2^(e - 1023)` (since it always generates a normal number). `E[(x + 1) >> 1]` exactly equals to 2^-51, and `E[2^e] = 2^1022 * 2^-1 + 2^(1022-1) * 2^-2 + ... + 2^(1022-74) * 2^-75 + 2^(1022-75) * 2^-75 = 2^1023 (1/4^1 + 1/4^2 + ... + 1/4^75 + 0.5/4^75) = 2^1023 (1/3 + 1/(6*4^75))`. So `E[d] = (1 + 2^-51 * 2^52) * (1/3 + 1/(6*4^75)) = 1/2 + 1/4^76`.

Sharlin 3y ago

Often the naive method is all right as you aren’t really concerned by the fact that numbers close to zero don’t have the maximum precision possible – particularly when your use case has a certain "native scale" and what amounts to a fixed point approach works fine. And especially if the results are going to be discretized at some point anyway (eg. computer graphics).

camel-cdr 3y ago

The holy grail of random float generation is a function that takes two bounds and can generate ALL representable floating point numbers within the given bounds with the proper distribution, here uniform (and does so efficiently). I've written such an implementation with explanations [0], but it doesn't handle subnormals correctly. I've yet to find the time to revisit the problem. [0] https://github.com/camel-cdr/cauldron/blob/main/cauldron/random.h#L1462 https://github.com/camel-cdr/cauldron/blob/main/cauldron/ran...

NeutralForest 3y ago
what does subnormal mean?
- camel-cdr 3y ago
  Usually floating point numbers are defined as `1.fraction * 2^(exponent-bias)`, but if the exponent is zero `0.fraction * 2^(1-bias)` is used, to allow representing very small numbers. IIRC I also explained it here: https://m.youtube.com/watch?v=VHJUlRiRDCY&t=2774s https://m.youtube.com/watch?v=VHJUlRiRDCY&t=2774s
  - NeutralForest 3y ago
    thanks!
  - Aardwolf 3y ago
    And the reason subnormals exist at all is that normally there's an implicit 1-bit (without the implicit 1-bit, there would be multiple representations of the same value, e.g. all values with mantissa 0 would represent 0 regardless of exponent), but you could never represent the value 0 itself if the implicit 1-bit is always there. So for the lowest exponent value, the implicit bit is set to 0, to allow representing 0, but also those other very tiny values since there's still an entire mantissa to give meaning to in addition to value 0.
    - magicalhippo 3y ago

jkafjanvnfaf 3y ago

After thinking about the problem a bit, I've come up with a completely different approach. I'm sure this is known, but I can't find anything so I'd appreciate if someone can point me to a citation. Generate two independent doubles a and b the "usual way", meaning from [0, 1) with equal spacings of 2^-53. Then calculate 2^-53*b + a. Assuming for simplicity that the rounding mode is towards zero, this gives a new random double in [0, 1). I've only done the math in my head so I'm not 100% sure on all the numbers, but I believe it scores pretty well at the criteria in the blog post; tentative results are Probability of zero: 2^-106 Smallest nonzero: 2^-106 Largest non-seeable: approx. 2^-54 Distinct values: 2^58.7 The usual rounding mode (round to nearest, ties to even) will generate numbers in [0, 1] with very slight biases that need closer analysis. One could also go further and do e.g. 2^-106*c + 2^-53*b + a, but something like that should only be needed for for single-precision floats. The obvious downside is that you always need two random integers and not just in the slow path. However, SIMD becomes easier because you don't need AVX512 for lzcnt.

10000truths 3y ago

This is exactly how random floats are computed in the Zig standard library: https://github.com/ziglang/zig/pull/10428/ https://github.com/ziglang/zig/pull/10428/

leiroigh 3y ago
This is awesome and the first precedent I've ever seen for a standard library doing the right thing on rand floats. Big kudos to the zig people and thanks for brightening my day!
corsix 3y ago
Oh, cute. It looks like they ever-so-slightly overweight the probability of values whose mantissa is entirely zeroes though. For example, the probability of hitting exactly 0.5 should be 2^-54 + 2^-55, whereas zig looks to give it 2^-54 + 2^-54.

kwillets 3y ago

The final methods seem equivalent to reading a non-float binary fractional number .xyz... from a random bitstream and stopping at the point where further digits will be truncated in the float representation. eg if x = 1 read 53 more mantissa bits and stop, if x = 0 decrement the exponent and repeat. There's no bias in such a process.