generate FP bit patterns with the same probability distribution as you'd get from rounding genuinely random real numbers to double
Ugh, though, thinking about this, that probability distribution is actually pretty weird, because of the default round-to-nearest semantics! The value 1.0 is a viable output because things can round up to it; weird things happen at exponent boundaries, where half as large a range of reals can round up to each smaller power of 2 as can round down to it; denormals complicate things even further. If I've thought it through correctly, we can get 1.0 with probability 2−54; each float strictly between 0.5 and 1.0 with probability 2−53; 0.5 itself with probability 3 × 2−55, then the next set of nonzero ones with prob 2−54, and so on, so that each power of 2 occurs with a probability that's the arithmetic mean of the probabilities of the numbers next to it. And finally denormals are treated as an extension of the lowest exponent range, so that there's no glitch anywhere near the smallest normalised number, but 0.0 itself can occur with exactly half the probability of the denormal and exponent-1 values just above it.
(At least for this purpose we can ignore the round-to-even tie-breaking case, because it occurs with probability literally zero.)
For all that this is the Platonically-right interpretation of the normal IEEE 754 semantics, it's clearly too horrible to use! Although spigot can implement it already if you pipe random input into it, along the lines of
no subject
Date: 2023-06-24 09:35 (UTC)Ugh, though, thinking about this, that probability distribution is actually pretty weird, because of the default round-to-nearest semantics! The value 1.0 is a viable output because things can round up to it; weird things happen at exponent boundaries, where half as large a range of reals can round up to each smaller power of 2 as can round down to it; denormals complicate things even further. If I've thought it through correctly, we can get 1.0 with probability 2−54; each float strictly between 0.5 and 1.0 with probability 2−53; 0.5 itself with probability 3 × 2−55, then the next set of nonzero ones with prob 2−54, and so on, so that each power of 2 occurs with a probability that's the arithmetic mean of the probabilities of the numbers next to it. And finally denormals are treated as an extension of the lowest exponent range, so that there's no glitch anywhere near the smallest normalised number, but 0.0 itself can occur with exactly half the probability of the denormal and exponent-1 values just above it.
(At least for this purpose we can ignore the round-to-even tie-breaking case, because it occurs with probability literally zero.)
For all that this is the Platonically-right interpretation of the normal IEEE 754 semantics, it's clearly too horrible to use! Although
spigot
can implement it already if you pipe random input into it, along the lines of