Hacker News Reader

Comment by j_not_j | Hacker News Reader

[flagged]

> If it is not [0,1) then it's not useful.

I can understand [0, 1) being useful in some use cases but saying it's entirely useless is a bit dramatic, don't you think? I've certainly had uses for [0, 1].

weirdly enough, the provided implementation of rand_float() actually generates numbers in the interval [0, 1).

I haven't read the implementation of _generate_canonical_generic(), but the following special case in generate_canonical() applies here:

  if constexpr (prng_is_bit_uniform && sizeof(T) == 4 && sizeof(generated_type) == 8) {
    return (static_cast<std::uint32_t>(gen()) >> exponent_bits_32) * mantissa_hex_32;
  }

which boils down to:

  return (static_cast<std::uint32_t>(gen()) >> 8) * 0x1.p-24;

that is, a number in the interval [0, 2^24) divided by 2^24.

Thank you for noticing! Turns out [0, 1] is an artifact of the old documentation carried back from the time where generic GCC / clang approach used to produce occasional 1's due to some issues in N4958 specification. This is fixed in a new commit.

For floats there are essentially 3 approaches that are selected based on the provided range & PRNG:

  1) Shift + multiply (like `(rng >> 11) * 0x1.0p-53`)
  2) "Low-high" from the paper by J. Doornik
  https://www.doornik.com/research/randomdouble.pdf
  3) Canonical GCC implementation

All of these generate values in [0, 1) range which is now reflected in the docs properly. For most actual cases 1st method is the one selected.

In this case I would suggest using the high bits of the RNG output when generating a float, since some generators have better entropy around the high bits.

So when you're generating floats with a 64-bit generator, instead of masking out the high bits with the static_cast, you may want to use the following:

  return (gen() >> 40) * 0x1.0p-24;

This article shows how to generate perfectly-distributed random floats: https://specbranch.com/posts/fp-rand/

That page ends with "This is the first efficient algorithm for generating random uniform floating-point numbers that can access the entire range of floating point outputs with correct probabilities". I am very skeptical.

The same basic idea is described here [1] (though without any pretty pictures, but with a link to working code), and I know others have implemented the same insight [2].

This seems like the kind of thing that is reinvented every time some who cares about random numbers learns how IEEE 754 works.

[1] http://mumble.net/~campbell/2014/04/28/uniform-random-float

[2] https://github.com/camel-cdr/cauldron/blob/a673553dd7925b0f1...

Downey published a (big-endian, SPARC) method in 2007:

https://allendowney.com/research/rand/

And Lemire in 2017 asked the question

https://lemire.me/blog/2017/02/28/how-many-floating-point-nu...

The basic approach: collect top 60 bits of a 64-bit PRNG. (Assume the LSBs are corrupt or zero or nonrandom). Set exponent zero. If the top mantissa bit is zero, shift left, subtract one from exponent. Repeat. When you run out of bits, collect another PRNG from your generator and resume until you have 56 bits of mantissa. When done, your floating-point PRNG is 2^exponent * mantissa.

My short explanation: Suppose you want a random distance. Multiplying a PR integer is the same as selecting a random tile in the distance, and measuring the distance to the (far) edge of the tile. This is the multiply method.

But if you want a real-valued distance even for very near distances, you need to scale your random number and ensure you have random bits throughout the mantissa. So reduce the exponent for every leading zero in your PR integer and shift. Add more bits if needed.

Test: the histogram of exponents for a large-enough set of samples is linear, give or take. Very small floating numbers (less than say 1/32768) are 1/16 as likely as numbers in [0.5,1).

Now that I've actually looked at the "utl::random" code in the OP, I see that its UniformRealDistribution is a wrapper around std::generate_canonical, so the juicy bits about turning a random into into a random float are not exposed here at all. But the utl::random code does include an pointer* to an informative C++ working group note.

* https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p09...

Thanks for the additional pointers!

I like the idea of looking at the histogram of exponents: incrementing (by 1) the exponent doubles the width of the interval so should double the number of hits. Conversely the histogram of the significand (or any subset of its bits) should be flat.

Beebe publishes a lengthy bibliography.

https://www.netlib.org/tex/bib/prng.pdf

758 pages and counting....

Lengthy indeed. Do you have advice for how to find amongst these the first paper on correctly sampling [0,1) with real-world floating point?

Reply on news.ycombinator.com