How to generate a double-precision floating-point number in [0, 1] uniformly

If you can live with [0, 1), there are simpler ways of doing this. Generate a number [1,2) by filling the 52-bit fractional component with a random integer. Then subtract 1. And you're done.

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/dSF...

I appreciate the effort to write this up, but I'm not seeing that this has applicability to any real problem. If you put 53 random bits in the mantissa with a fixed exponent, you get a smooth distribution from from [0, 1 - 2^-53] in increments of 2^-53. He seems to be concerned that you will never get a number between(0, 2^-53). True, but you'll also never get any other number between (n * 2^-53,(n+1) * 2^-53) for all the rest of the range.

[wrong statement about denormals deleted]

His talk about "dividing" to produce the desired range seems misguided. Floating point division is ~10x the cost of any other floating point operations. If you are trying to produce precise floating point numbers, you almost certainly want to be working at the level of bits rather than hoping that division will do the work for you.

First off this seems like unordered notes and I spent 15 minutes trying to figure out what the actual requirements were, because it was saying things it needed and then talking about solutions that didn't have those features.

Anyway, the main goal seems to be even density of chosen numbers, while also not skipping floats. In other words the probability of any specific float being chosen is proportional to its width.

This is easy to do, and needs no rounding. At least for [0,1)

Observe that floats [.5,1) are half the probability space. Use one random bit to choose between them and smaller values.

Then floats [.25,5) are half the space. Use another random bit.

Once you've chosen the exponent, fill in the mantissa with 52 random bits. All floats in the range are the same distance apart, so random bits work fine.

  mantissa := random(2^52)
  exponent := -1
  while random(2) == 1:
    exponent -= 1
  return (1 + mantissa/2^52) * 2^exponent

If you really want [0,1], then you must first define how wide of a range of real numbers you want to assign to "1". Rounding toward zero makes more sense. Alternatively, note how nonzero double precision floats go down below 2^-1000. You will never ever actually generate a zero, no matter what you do. Why should "1" be more likely than "0"?

But if you insist on rounding to nearest then you can put in a special clause where a mantissa of 0 has a 50% chance of being treated as 2^52.

I like the idea of an random number generator that exploits the full precision of the floats. However, the issue with the naive method is loss of accuracy, not bias. When the article states

If k is 64, a natural choice for the source of bits, then because the set {0, 1, ..., 2^53 - 1} is represented exactly, whereas many subsets of {2^53, 2^53 + 1, ..., 2^64 - 1} are rounded to common floating-point numbers, outputs in [0, 2^-11) are underrepresented.

this is very misleading, because the points that are "underrepresented" are closer together. The discrete probability density created by the naive method is a perfectly good approximation to the uniform distribution.

Furthermore, although as the article shows, we can do better, this can, and must, come at a cost. Since the small numbers occur with much lower probability, the only way to generate them (without importance sampling) is by generating much more randomness (1088 bits in the article). Generating random bits is expensive, and this uses 17 times more randomness than the naive method.

EDIT: it would be possible to avoid generating all the extra bits in most cases, and so an efficient average time algorithm should be possible, but the worst case still requires 1088 bits to be generated.

For [0, 1), would the naive approach of generating a string of, say, 16 digits, sticking a "0." in front of it and calling strod / atof yield similarly uniform results?

Wouldn't the best (and coincidentally simplest) way be to synthesize the double-precision word by combining a random integer exponent, a random integer mantissa, and a random sign bit? You'd need to account for the logarithmic distribution of the exponent, but that seems like the biggest complication.

I don't understand the debate about rounding. What does it mean to round a random number? If the rounding process is biased, the number will be less random as a result, so you shouldn't have rounded it. If the process isn't biased, the number will be no more or less randomly-distributed than it was before, so there was no point in rounding it.

I'm sure that there are valid reasons to want to do this, but if your find yourself needing to do this, I think you should at least consider whether there might be some other way to get what you want.

Using a computer algorithm, http://en.wikipedia.org/wiki/Mersenne_twister is a good choice.

"Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin." -- Von Neumann.

Using a computer, maybe look for small current variations around a resistor or similar device? http://en.wikipedia.org/wiki/Hardware_random_number_generato...