Proper way to generate a random float given a binary random number generator?

Question

Let's say we have a binary random number generator, int r(); that will return a zero or a one both with propability 0.5.

I looked at Boost.Random, and they generate, say, 32 bits and do something like this (pseudocode):

x = double(rand_int32());
return min + x / (2^32) * (max - min);

I have some serious doubts about this. A double has 53 bits of mantissa, and 32 bits can never properly generate a fully random mantissa, among other things such as rounding errors, etc.

What would be a fast way to create a uniformly distributed float or double in the half-open range [min, max), assuming IEEE754? The emphasis here lies on correctness of distribution, not speed.

To properly define correct, the correct distribution would be equal to the one that we would get if we would take an infinitely precise uniformly distributed random number generator and for each number we would round to the nearest IEEE754 representation, if that representation would still be within [min, max), otherwise the number would not count for the distribution.

P.S.: I would be interested in correct solutions for open ranges as well.

Magenta

I assume "take an infinitely precise uniform random number generator on the range [min, max)", because without a range restrictions, there are no uniform random number generators. :) – Garmon

@Yakk That would create a malformed distrbution on the very ends due to the rounding, [min - n*ULP, max + n*ULP) would be fine though for some value of n that I'm too lazy to think about but is probably 1. – Magenta

@nightcracker To make sure I'm understanding correctly: You want a dsitrbution where some floating point bit patterns are much more common than others (due to being farther away from 0), NOT a distribution where each floating point bit pattern is equally probable? – Derwon

@MarkB yes, nightcracker does. – Garmon

@nightcracker ok, uniform on [std::numeric_limits<double>::Min(), std::numeric_limits<double>::Max()), rounded, then discarded if outside of [min, max). ;) – Garmon

@Mark B: No. I assume he does not want this. He is confused by the fact that the sampling interval is larger in some spots. But his is not a problem. For equidistribution on (b-a) you have P(x \in (a,b)) = b-a. This holds up to an error of the sampling intervall if you use equi distributed bit patterns. – Glacialist

@MarkB Exactly, I want the distribution of IEEE754 numbers to be such that, when converted back into "true numbers" it closely resembles a uniform distribution. This means that some bit patterns are much more common than others because numbers close to 0 are much more dense. – Magenta

@nightcracker: Ah. Sorry. Maybe I was confusred. For denormalized numbers the method I had in mind works, because they are just a series of {0,1} * 1/2^k. Maybe it does not work if min and max are different form 0,1. – Glacialist

@ChristianFries Yes, min and max can be any range. – Magenta

But then min + x * (max-min) where x is a denormalized number generated from equip-distributed bit pattern is correct - isn't it? – Glacialist

@ChristianFries That is my problem/question - I don't exactly know. Worse than that, I don't really trust processor floating point units enough when it comes to multiplying really small numbers (x) with normal numbers when precision is paramount. – Magenta

Answer 1

AFAIK, the correct (and probably also fastest) way is to first create a 64 bit unsigned integer where the 52 fraction bits are random bits, and the exponent is 1023, which if type punned into a (IEEE 754) double will be a uniformly distributed random value in the range [1.0, 2.0). So the last step is to subtract 1.0 from that, resulting in a uniformly distributed random double value in the range [0.0, 1.0).

In pseudo code:

rndDouble = bitCastUInt64ToDouble(1023 << 52 | rndUInt64 & 0xfffffffffffff) - 1.0

This method is mentioned here: http://xoroshiro.di.unimi.it (See "Generating uniform doubles in the unit interval")

EDIT: The recommended method has since changed to: (x >> 11) * (1. / (UINT64_C(1) << 53))

See above link for details.

Vimen

Just wanted to add that if you want to generate a single precision 32bit float (what "float" means in most languages), the formula becomes: rndFloat = bitCastUInt32ToDouble(127 << 23 | rndUInt32 >> 9) - 1.0 (here, rndUInt32 >> 9 is the same as rndUInt32 & 0x1ffffff) – Aurist

Re: > The recommended change has since changed ... I'm not convinced that "you will be generating half the values you could actually generate" due to "all doubles generated will have the lowest significand bit set to zero". I mean, we are already throwing away bits to make 64 into 52; who cares? I am seeing the issue, BTW. Looking at hex dumps of successive generated values, I see that the least significant byte in the generate double is never an odd number. – Counterproposal

Can't we fix this, I wonder, like this: do the subtraction of 1.0 against the original union member. Then access it as an integer again, to put a random bit into the LSB position. Then extract the final double from the union. We can use any of the bits that we threw away from the exponent field; or just the LSB of the original integer value. – Counterproposal

It seems every so slightly non-uniform to sneak the bit back into the value this way, because after the subtraction, the values are no longer uniformly distributed within the double representation. Whenever we add a 1 bit, we are effectively adding a small epsilon to the value. That epsilon is not a constant over the range, though! Its value depends on whatever the exponent is. – Counterproposal

Answer 2

Here is a correct approach with no attempt at efficiency.

We start with a bignum class, and then a rational wrapper of said bignums.

We produce a range "sufficiently bigger than" our [min, max) range, so that rounding of our smaller_min and bigger_max produces floating point values outside that range, in our rational built on the bignum.

Now we subdivide the range into two parts perfectly down the middle (which we can do, as we have a rational bignum system). We pick one of the two parts at random.

If, after rounding, the top and bottom of the picked range would be (A) outside of [min, max) (on the same side, mind you!) you reject and restart from the beginning.

If (B) the top and bottom of your range rounds to the same double (or float if you are returning a float), you are done, and you return this value.

Otherwise (C) you recurse on this new, smaller range (subdivide, pick randomly, test).

There are no guarantees that this procedure halts, because you can either constantly drill down to the "edge" between two rounding doubles, or you could constantly pick values outside of the [min, max) range. The probability of this happening is (never halting), however, zero (assuming a good random number generator, and a [min, max) of non-zero size).

This also works for (min, max), or even picking a number in the rounded sufficiently fat Cantor set. So long as the measure of the valid range of reals that round to the correct floating point values is non zero, and the range has a compact support, this procedure can be run and has a probability of 100% of terminating, but no hard upper bound on the time it takes can be made.

Garmon

Sounds solid. Is there any rounding algorithm you would recommend for the rational -> IEEE754 conversion? – Magenta

@nightcracker Rounding down would work (give a valid distribution): it really depends on what "region" you want each IEEE754 to represent. I guess we could look at the IEEE754 spec to figure out how it rounds when you do something like multiplication? That would give us a clue asto what region they are "supposed" to represent. – Garmon

IEEE754 arithmetic is always "arithmetic at infinite precision followed by rounding", i.e. the IEEE754 number representing ab is the one being closest to ab (put differently, you can think of a*b calculated at infinite precision, then rounded to the closes IEEE754 number). Hence I believe that min + (max-min) * x where x is a denomalized number is a feasible approach. – Glacialist

@ChristianFries "by rounding" -- what happens to the midway point? Up, down, even, odd? I assume linear rounding, not exponential? – Garmon

@nightcracker I have come up with an alternative approach: you start generating starting with a sign bit, then the MSB of the largest (in magnitude) element of the range. You stop generating bits when the value lies completely within the rounding range of a particular floating point value (at which point you are done), or you the value generated lies completely outside the range you are targeting (in which case, you discard and repeat). This ends up being logically equivalent to the above method with a particular range (bits are equivalent to which half subrange), but less abstract. – Garmon

@yakk: I see. So what he is looking for is just a "random rounding" which is in accordance of the distribution. For the denormalized numbers this can be achieved by just generating one bit more and use it as a rounding indicator. However: I doubt if this will solve a real problem. Who needs equidistributed random numbers between min and max? Usually they are input to some transformation which generates another distribution (e.g. normal). Then you face the same rounding problem again distorting the distribution. Then the distribution is just correct up to "rounding errors". – Glacialist

@Yakk Would you mind providing pseudocode of your alternative approach in the answer? Then I would accept it. Right now I have some troubles with converting it to code. – Magenta

@nightcracker sure. What bignum library did you try to work with? – Garmon

@Yakk I was under the impression your alternative approach didn't require a bignum library and you could just directly generate bits inside the float. – Magenta

@nightcracker oh, that approach. Hmm, easiest way would be to generate an unsigned sequence of bits, with a fixed max precision (discounting leading 0s), and a min value. Reducimg that to float generation will be feasible. On my phone, will try to remember on Monday if this comment is not enough. – Garmon

Answer 3

The problem here is that in IEEE754 the doubles which may be represented are not equi-distributed. That is, if we have a generator generating real numbers, say in (0,1) and then map to IEEE754 representable numbers, the result will not be equi-distributed.

Thus, we have to define "equi-distribution". That said, assuming that each IEEE754 number is just a representative for the probability of lying in the interval defined by the IEEE754 rounding, the procedure of first generating equi-distributed "numbers" and the round to IEEE754 will generate (by definition) an "equi-distribution" of IEEE754 numbers.

Hence, I believe that the above formula will become arbitrary close to such a distribution if we just choose the accuracy high enough. If we restrict the problem to finding a number in [0,1) this means to restricting to the set of denomalized IEEE 754 numbers, which are one-to-one to a 53 bit integer. Thus it should be fast and correct to generate just the mantissa by a 53 bit binary random number generator.

IEEE 754 arithmetic is always "arithmetic at infinite precision followed by rounding", i.e. the IEEE754 number representing ab is the one being closest to ab (put differently, you can think of a*b calculated at infinite precision, then rounded to the closes IEEE754 number). Hence I believe that min + (max-min) * x, where x is a denomalized number, is a feasible approach.

(Note: As clear from my comment, I was first not aware that you where pointing to the case with min and max different from 0,1. The denormalized numbers have the property that they are evenly spaced. Hence you get the equi distribution by mapping the 53 bits to the mantissa. Next you can use the floating point arithmetic, due fact that it is correct up to machine precistion. If you use the reverse mapping you will recover the equi-distribution.

See this question for another aspect of this problem: Scaling Int uniform random range into Double one

Glacialist

This problem is exactly the reason I asked this question. There are many more IEEE754 numbers close to 0 than there are ones close to 1, for example. – Magenta

Yes. But what I wanted to say is: an IEEE 754 number is a representative of its rounding interval (in the sense of an equivalence relations) and the random number represents the event of "having a real number falling in that interval". Endowed with this notion you are sampling an equi-distribution. Its just that your sampling intervals are not equally-spaced - but this does not harm the fact that you have sampled an equi-distribution. – Glacialist

Answer 4

std::uniform_real_distribution.

There's a really good talk by S.T.L. from this year’s Going Native conference that explains why you should use the standard distributions whenever possible. In short, hand-rolled code tends to be of laughably poor quality (think std::rand() % 100), or have more subtle uniformity flaws, such as in (std::rand() * 1.0 / RAND_MAX) * 99, which is the example given in the talk and is a special case of the code posted in the question.

EDIT: I took a look at libstdc++’s implementation of std::uniform_real_distribution, and this is what I found:

The implementation produces a number in the range [dist_min, dist_max) by using a simple linear transformation from some number produced in the range [0, 1). It generates this source number using std::generate_canonical, the implementation of which my be found here (at the end of the file). std::generate_canonical determines the number of times (denoted as k) the range of the distribution, expressed as an integer and denoted here as r*, will fit in the mantissa of the target type. What it then does is essentially to generate one number in [0, r) for each r-sized segment of the mantissa and, using arithmetic, populate each segment accordingly. The formula for the resulting value may be expressed as

Σ(i=0, k-1, X/(r^i))

where X is a stochastic variable in [0, r). Each division by the range is equivalent to a shift by the number of bits used to represent it (i.e., log2(r)), and so fills the corresponding mantissa segment. This way, the whole of the precision of the target type is used, and since the range of the result is [0, 1), the exponent remains 0** (modulo bias) and you don’t get the uniformity issues you have when you start messing with the exponent.

I would not trust implicity that this method is cryptographically secure (and I have suspicions about possible off-by-one errors in the calculation of the size of r), but I imagine it is significantly more reliable in terms of uniformity than the Boost implementation you posted, and definitely better than fiddling about with std::rand.

It may be worth noting that the Boost code is in fact a degenerate case of this algorithm where k = 1, meaning that it is equivalent if the input range requires at least 23 bits to represent its size (IEE 754 single-precision) or at least 52 bits (double-precision). This means a minimum range of ~8.4 million or ~4.5e15, respectively. In light of this information, I don’t think that if you’re using a binary generator, the Boost implementation is quite going to cut it.

After a brief look at libc++’s implementation, it looks like they are using what is the same algorithm, implemented slightly differently.

(*) r is actually the range of the input plus one. This allows using the max value of the urng as valid input.

(**) Strictly speaking, the encoded exponent is not 0, as IEEE 754 encodes an implicit leading 1 before the radix of the significand. Conceptually, however, this is irrelevant to this algorithm.

Organza

-1 With all due respect, this is not helpful at all. My question explicitly states that I don't trust the Boost.Random implementation, and std::uniform_real_distribution is copied pretty much from Boost.Random. And std::uniform_real_distribution is just that - a name. I'm looking for an actual algorithm. – Magenta

@nightcracker The interface to the standard random library is definitely heavily influenced by Boost, but that doesn’t necessarily mean that the implementation is. In the case of libstdc++, it appears to use std::generate_canonical to map the generator’s output to [0, 1) and does a linear transformation from that to [min, max). The implementation of std::generate_canonical is rather interesting, and I’ll see if I can update with an explanation. – Organza