Hacker News Reader

Comment by camel-cdr | Hacker News Reader

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/ran...

View on HN

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.

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.

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?

But what if you want a random number between two very small floats?

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

> 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?

Uniform over the real numbers.

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?)

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.

How much entropy do you need for that?

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.

what does subnormal mean?

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

thanks!

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.

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.

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.

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)

Depends on control registers like e.g. MXCSR. It's an utter mess. Consider e.g. https://news.ycombinator.com/item?id=32738206

> 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.

It's been over a decade, so I struggle to recall the details.

> 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.

Yes thinking about you're probably correct.

Reply on news.ycombinator.com