4 ms·

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

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
      Subnormals or denormalized floats can lead to some unexpected results. IIRC dividing a subnormal by a subnormal returns infinity. I experienced this where a colleague had tried to reduce the number of divisions by rewriting (a/x)*(b/y)*(c/z) to (a*b*c)/(x*y*z). The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1.
      - bee_rider 3y ago
        Yeah, I think a good general guideline is that if you see a bunch of divisions like that, the author did it intentionally (to avoid overflow or subnormals). Anyway, I wonder if the optimization even helps at all. Either way the critical path is a division and two multiplications. That seems like the sort of thing that, between SIMD an OoO execution, a modern CPU ought to be able to work out, haha. No promises but I’d want to benchmark it before caring.
      - Aardwolf 3y ago
        Maybe depends on the platform? I tried 5e-323 / 6e-323 in javascript console and it returns 0.8333333333333334 (and e-323 should be subnormal for double precision)
        leiroigh 3y ago
        Depends on control registers like e.g. MXCSR. It's an utter mess. Consider e.g. https://news.ycombinator.com/item?id=32738206 https://news.ycombinator.com/item?id=32738206
      - thethirdone 3y ago
        > IIRC dividing a subnormal by a subnormal returns infinity. This should not be true in a conformant implementation of IEEE 754-2008. You can get infinity by dividing by a subnormal (due to overflow), but that should only be possible with sufficiently large numerators. > The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1. My guess is that this problem was actually the numerator being subnormal and the denominator underflowing all the way to zero and therefore the result going to infinity.
leiroigh 3y ago
I don't think that would be a good idea -- sticking to positive normals and truncating below is enough for float32 and float64. I mean, consider that 1.1754944f-38 is a normal 32 bit float. The probability that your "perfect" code would ever get triggered is such that humankind has not built enough compute to expect to get that unlucky yet. (you may ask: why then care about that possibility if it's dead code anyway? Because of security vulns in context of flawed RNGs. Bad underlying RNG leading to bad distributions is expected; bad underlying RNG leading to RCE is not OK. So don't ever output zero or subnormals!) For bfloat16 the same holds; for ieee float16, well, I have no clue what we want.
- camel-cdr 3y ago
  I don't think such a function would be a good default, but I think there may be applications, probably some complex very large simulations, where this might matter.
  - leiroigh 3y ago
    You think there might be applications like complex very large simulations where an event with probability 1:10^38 matters? You are aware that the current age of the universe is < 10^27 nanoseconds, just for comparison?
    - camel-cdr 3y ago
      But what if you want a random number between two very small floats?
      - leiroigh 3y ago
        Fair enough, if a user asks for a random float between [0, 1e-38f] then subnormals are expected. I was just thinking about the (0,1) case, under the mistaken assumption that one could map it to (a,b) via multiplication/addition, but you're right -- if you want (a,b) perfectly () then it's not obvious to me. () up to inaccuracies of cosmologically negligible scale
ot 3y ago
How much entropy do you need for that?
- hansvm 3y ago
  Average entropy is low because the events are rare. Worst-case though, you need a number of bits linear in the number of states the exponent can have, plus the number of bits of precision you need. Subnormality doesn't make that much worse because they all have the same exponent (the thing that consumes a lot of worst-case bits), and because their relative frequencies are precisely what you would obtain by filling the significand with uniformly random bits. It's been awhile since I've thought about this problem, so I'll reserve the right to be off by a small constant additive or multiplicative factor. You need around 24+256 bits for 32-bit IEEE-754 floats (less by roughly half if you're sticking in a range like 0 to 1), and 53+2048 for 64-bit. The average number of bits required is something closer to 25 and 54, respectively.
toast0 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). Uniform as in all representable numbers have equal probability? Or uniform as in approximating a uniform distribution over real numbers?
- camel-cdr 3y ago
  Uniform over the real numbers.
  - StefanKarpinski 3y ago
    It’s a little unclear what you mean by that without further explanation. Do you mean that conceptually one selects a real number at random and then rounds that real number to the closest representable float? (And if so, which rounding mode?)
    - zamadatix 3y ago
      Floats aren't uniformly distributed over the reals (e.g. larger floats are sparser) so the above approach simply tries to account for that when generating the number by using a precisely counteracting non-uniform selection method which will result in even distribution in relation to where in the reals the picks are. It does so directly.

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
    Subnormals or denormalized floats can lead to some unexpected results. IIRC dividing a subnormal by a subnormal returns infinity. I experienced this where a colleague had tried to reduce the number of divisions by rewriting (a/x)*(b/y)*(c/z) to (a*b*c)/(x*y*z). The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1.
    - bee_rider 3y ago
      Yeah, I think a good general guideline is that if you see a bunch of divisions like that, the author did it intentionally (to avoid overflow or subnormals). Anyway, I wonder if the optimization even helps at all. Either way the critical path is a division and two multiplications. That seems like the sort of thing that, between SIMD an OoO execution, a modern CPU ought to be able to work out, haha. No promises but I’d want to benchmark it before caring.
    - Aardwolf 3y ago
      Maybe depends on the platform? I tried 5e-323 / 6e-323 in javascript console and it returns 0.8333333333333334 (and e-323 should be subnormal for double precision)
      - leiroigh 3y ago
        Depends on control registers like e.g. MXCSR. It's an utter mess. Consider e.g. https://news.ycombinator.com/item?id=32738206 https://news.ycombinator.com/item?id=32738206
    - thethirdone 3y ago
      > IIRC dividing a subnormal by a subnormal returns infinity. This should not be true in a conformant implementation of IEEE 754-2008. You can get infinity by dividing by a subnormal (due to overflow), but that should only be possible with sufficiently large numerators. > The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1. My guess is that this problem was actually the numerator being subnormal and the denominator underflowing all the way to zero and therefore the result going to infinity.

leiroigh 3y ago

I don't think that would be a good idea -- sticking to positive normals and truncating below is enough for float32 and float64. I mean, consider that 1.1754944f-38 is a normal 32 bit float. The probability that your "perfect" code would ever get triggered is such that humankind has not built enough compute to expect to get that unlucky yet. (you may ask: why then care about that possibility if it's dead code anyway? Because of security vulns in context of flawed RNGs. Bad underlying RNG leading to bad distributions is expected; bad underlying RNG leading to RCE is not OK. So don't ever output zero or subnormals!) For bfloat16 the same holds; for ieee float16, well, I have no clue what we want.

camel-cdr 3y ago
I don't think such a function would be a good default, but I think there may be applications, probably some complex very large simulations, where this might matter.
- leiroigh 3y ago
  You think there might be applications like complex very large simulations where an event with probability 1:10^38 matters? You are aware that the current age of the universe is < 10^27 nanoseconds, just for comparison?
  - camel-cdr 3y ago
    But what if you want a random number between two very small floats?
    - leiroigh 3y ago
      Fair enough, if a user asks for a random float between [0, 1e-38f] then subnormals are expected. I was just thinking about the (0,1) case, under the mistaken assumption that one could map it to (a,b) via multiplication/addition, but you're right -- if you want (a,b) perfectly () then it's not obvious to me. () up to inaccuracies of cosmologically negligible scale

ot 3y ago

How much entropy do you need for that?

hansvm 3y ago
Average entropy is low because the events are rare. Worst-case though, you need a number of bits linear in the number of states the exponent can have, plus the number of bits of precision you need. Subnormality doesn't make that much worse because they all have the same exponent (the thing that consumes a lot of worst-case bits), and because their relative frequencies are precisely what you would obtain by filling the significand with uniformly random bits. It's been awhile since I've thought about this problem, so I'll reserve the right to be off by a small constant additive or multiplicative factor. You need around 24+256 bits for 32-bit IEEE-754 floats (less by roughly half if you're sticking in a range like 0 to 1), and 53+2048 for 64-bit. The average number of bits required is something closer to 25 and 54, respectively.

toast0 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). Uniform as in all representable numbers have equal probability? Or uniform as in approximating a uniform distribution over real numbers?

camel-cdr 3y ago
Uniform over the real numbers.
- StefanKarpinski 3y ago
  It’s a little unclear what you mean by that without further explanation. Do you mean that conceptually one selects a real number at random and then rounds that real number to the closest representable float? (And if so, which rounding mode?)
  - zamadatix 3y ago
    Floats aren't uniformly distributed over the reals (e.g. larger floats are sparser) so the above approach simply tries to account for that when generating the number by using a precisely counteracting non-uniform selection method which will result in even distribution in relation to where in the reals the picks are. It does so directly.