State of the art of random double generation | Hacker News Reader

State of the art of random double generation

blahonga.org

43 pointsnkurz38 comments

This 1:2:1 distribution is the right sort of behavior for a double-precision random number generator (what the standard says notwithstanding). It's what you get if you throw a uniform random dart at the real line and round it to the nearest floating point number.

Also, a <= x <= b is really the right way for uniform floating point ranges to be, with half-probability of the edges getting selected. The API of returning [0, 1) is pretty horrid, because then when somebody uses that to generate [1, 2), by doing rand01() + 1.0, they actually create the possibility of generating 2, because (1-2^-53) + 1 rounds up to 2. (That's just like the rounding that happens when generating to [2^52, 2^52+2].)

If you want a <= x < b, getting that for practical purposes could be relatively simple: compute rand01() * (b - a) + a, and if the output is >= b, then return a. That works wonderfully if rand01() behaves like randint(2^53) / 2.0^53. It's a little suboptimal if rand01() is implemented something like randint(2^64) / 2.0^64, because it creates a weird little density spike at the bottom for contrived ranges such as [2.0^-30, 1 + 2.0^-30), but if you're getting that anal you'd also worry about the round-to-even behavior of computing rand01() + 1.0.

"It's what you get if you throw a uniform random dart at the real line and round it to the nearest floating point number."

Only when you throw the dart, repeat until you are in the target interval, and then look for the nearest representable number. If you throw the dart, look for the nearest representable number and repeat until that number is in the interval, you get something different.

Worse, neither approach gives you a uniform distribution in the general case. The hypothetical uniform dart is as likely to land between 1 and 2 as it is to land between 2 and 3, but there are more representable numbers in the first range than in the second.

I also think your rand01() + 1.0 example is flawed. It is just because floating point is hard that this kind of functionality ends up in a standard library.

> Only when you throw the dart, repeat until you are in the target interval, and then look for the nearest representable number.

Yeah, when you throw your dart at the target interval.

> If you throw the dart, look for the nearest representable number and repeat until that number is in the interval, you get something different.

I don't know why you'd want to throw a uniform dart at some other, larger interval, before culling it down. You can't throw a uniform random dart at the entire real line. I think what you're describing gives undue representation to the edges -- you wouldn't use those sampling points with those weightings to integrate a function, but you might use the trapezoid method.

> The hypothetical uniform dart is as likely to land between 1 and 2 as it is to land between 2 and 3, but there are more representable numbers in the first range than in the second.

That's a fact of life.

At first, I thought trying to sample "uniformly" from a discrete non-homogeneous set of real numbers was just mathematically nonsense, but there's actually a pretty elegant generalization:

Basically, the probability of picking any particular floating-point value x should be proportional to the distance from x to the next largest representable number (with zero probability of choosing the largest floating-point number). This has the property that for any two floating-point numbers a and b, the probability of sampling a value in the interval [a, b) is proportional to b - a, which is the same as the defining property of the usual uniform distribution on real numbers, just with some additional restrictions.

It's no better in Java. A fairly literal translation of the code [1], using the new Random::doubles method [2], gets counts like so:

  166808
  332743
  500449

If anything, this is even weirder than the C++ results. :rage1:

The doubles method makes each value "as if it's the result of calling the following method with the origin and bound':

  double nextDouble(double origin, double bound) {
    double r = nextDouble();
    r = r * (bound - origin) + origin;
    if (r >= bound) // correct for rounding
      r = Math.nextDown(bound);
    return r;
  }

Which is the naive approach that Artur has been struggling against.

[1] https://bitbucket.org/snippets/twic/qqEe

[2] https://docs.oracle.com/javase/8/docs/api/java/util/Random.h...

I guess this could be related to my question I asked some time ago on SO: http://stackoverflow.com/questions/25668600/is-1-0-a-valid-o...

A related question:

Is there any efficient way to generate a random number in a range w.r.t. the amount of entropy used? The way I know (generate a random number, if it's in the top incomplete section regenerate until it isn't, modulo the width of the range, add the starting value - there are some fencepost errors here but you get the idea) is really inefficient in the case of, for example, you want a number between zero and half the range + 1 (so, for example, if you have an 8-bit RNG and you want numbers between 0-128, inclusive.)

You end up using, on average, 8*(1/(1-127/256)) ~= 15.9 bits of entropy per call, whereas you only "need" to use log2(129) ~= 7.0 bits of entropy per call.

Take the interval [0, 1), keep cutting it in half and picking a side, until the interval's contained by a [k/n, (k+1)/n).

Consider your case where n = 129. After 8 subdivisions, you have 256 possible intervals, which means the 128 interior boundary points (numbers of the from k/129, 0 < k < 129) will lie within half of the intervals, and the other half has definitively chosen a number. Those which haven't have a 50% chance of terminating after each successive bit, so you've got on average 9 bits of entropy per call.

I don't know if that's optimal, but you can't do better than 8.

You can do better than 8, as mark-r said.

Do as I mentioned, but batch requests.

For example, if you generate 89 numbers from 0-128 inclusive at once, you'll need, on average, 624.3 bits of entropy to generate them, whereas the optimal is 624.0 bits. That's an average of 7.0 bits of entropy per number generated. (It works out to something like 1.00054 x optimal.)

So, generate a number from 0 to 128^k-1, and split it base 129, for some value of k, returning the k random numbers generated this way. For sufficiently large k you will be able to use arbitrarily close to the optimal number of bits of entropy per request.

Unfortunately, this requires you know in advance what the range of the random numbers to generate will be.

You don't need to know in advance. Generate a big r from 0 <= r < s. When calling randint(n), return r mod n and update (r, s) = (r/n, s/n), unless r >= s - (s mod n) in which case reset s and regenerate r.

Sometimes you can generate the random numbers in batches to get better efficiency.

For your example, you can combine 3 random 8-bit values into a 24 bit value. You still only get 3 outputs, but you can use all the random values up to 7 * (129^3) - you only need regenerate random numbers 10% of the time, rather than 49%. 60 bits does even better, giving you 8 outputs and requiring regeneration only 0.2% of the time.

Nice.

A pity it requires knowing sizes in advance. I wonder... Could you do something similar "afterwards"?

Sure you can, I used a short program myself to calculate the best case outcomes. The question is whether it would be worth the bother. Usually the range you request is much smaller than the range of your random number generator, so the number of regenerations required are nowhere near the worst case - and the cost of generating a random number isn't enough to get worried about.

What good is the standard implementation if the gcc guy responds with WONTFIX?

Root of the problem seems to be that the standard gives pseudo code for generate_canonical (paragraph 26.5.7.2 in my unofficial version) that implementations must follow that, when translated straightforwardly into C, gives results that are inconsistent with the textual description of the function.

See http://stackoverflow.com/a/25669510

Language-lawyers-no-oxygen-astronauts (1) strike again.

Producing quality random floating point numbers, fast, is a solved problem (maybe even since the times of John von Neumann) but the astronauts would then have to come to the Earth and not just invent the wrong code attempting to use the the latest features of the standard. The direction should be the opposite: starting from the known and proper solutions, see how they can be (if they need to be) standardized.

The standard that standardizes "something" while not knowing how the problem is actually solved in practice is actually dangerous.

1) http://www.joelonsoftware.com/articles/fog0000000018.html

In my dark-tinted view from actual experience of working on multiple standards and the associated politics, I am inclined to say most standards are like this - on multiple occasions 95% of all people discussing what needs to go into a standard didn't even read the conceptual papers associated with them, yet ended up super vocal and forced their view. Maybe I was just unlucky, yet now I know what does it mean "design by committee" from an actual experience.

Author of the blog post following the sudden surge of traffic here.

> Producing quality random floating point numbers, fast, is a solved problem

Is it? Do you have any reference I can read?

Everything I've found is obviously flawed on the level of `(rand() % 1000)/1000.0`. I don't actually care about the C++ facilities that I'm discussing in that post. I just hoped that they would be a good path to follow to find one good source that discusses the various trade-offs between mapping a set of random integers (or a stream of bits, or whatever) to floating points. I've seen serious peer reviewed papers that just randomize the bits of the exponent and mantissa and call that random or papers that do what I've tried to do (which has its flaws) - count the number of representable, uniformly spaced, floating point numbers between `from` and `to` but then generate their random numbers with `random() % count`, which is like restoring the Mona Lisa with crayons.

So I'm completely serious. I'm desperate for a good source, where can I find it?

I get the impression that it's indeed not a solved problem, and even implementations that do not have to follow standards like C++ get it wrong. I recently submitted this to Hackernews, which is a nice treatment of the problems:

http://mumble.net/~campbell/2014/04/28/uniform-random-float?

I actually got the link to his code yesterday and it's eerily close to what I wrote a few months ago, but I rejected that approach because I didn't like that any of the numbers in [0.5,1.0) were so much more likely to be returned than 0.0. This gives the highest possible entropy and all numbers representable by floating points in that range can be returned, but the trade-off is that certain values will be returned vastly much more often than others. It's this discussion of trade-offs that I'm actually missing.

I don't feel competent enough to write that paper/article/rant about the trade-offs myself, I barely understand math and randomness, but I'm getting tired of all the libraries out there that say "uniform[1] random[2] floating points[3] [1] - not actually uniform [2] - not actually random [3] - terms and conditions apply, YMMV, not valid where invalid".

What do you mean when you say you don't like that some numbers are returned more often than others? This seems completely inevitable for a uniform distribution because floating point numbers don't have a uniform density across the real line.

Only if you choose a range like [1,2] with a uniform density of floating point numbers you will see all floating point numbers in the interval - except both bounds of course - with equal probability.

You would have to not return some floating point numbers ever to see all returned floating point numbers in the interval with equal probability but that of course hardly qualifies as a uniform distribution if some perfectly valid floating point numbers will never occur.

It totally qualifies as a uniform distribution. Is it the case that generating [1, 2) generates a uniform distribution, but subtracting 1 from that does not produce a uniform distribution?

I've been working on the same problem. I've found two broad solutions, which solve subtly different definitions of the problem:

1) Generate numbers that are uniformly distributed in the sense of, if you generated a bazillion of them and divided the range into buckets of even size, each bucket would have similar numbers of numbers in them.

2) Generate numbers that meet the above definition AND so that any particular floating number is no more likely to be picked than any other.

(The discrepancy being down to there being more representable floating point numbers in some ranges than others)

For my application, I don't need to worry about #2, so I just generate an integer number uniformly across a given range (which is simple enough to do without bias), cast to float, and divide. It meets the definition of #1 unless the buckets chosen are very small.

I have seen a few people write solutions for #2 online, but haven't seen any that have been sufficiently analysed and proven correct.

#2 is extremely easy, isn't it? Use a uniform random generator to generate a bit pattern of the size of a float, reinterpret_cast to float, and repeat until you have something that is in the target range (that gets rid of NaN's as they aren't ever in any range; filter denormalized numbers, if required)

You can speed that up by forcing bits that you know have to be zero or one to those values. If you do that, I think your expected number of iterations can be proven to be fairly small.

That would violate #1, as it would generate numbers more frequently that lay in ranges of higher precision.

But you can't simultaneously have #1 and #2. #2 means that you pick a number between 1 and 2 as often as you pick one between 2 and 4, 4 and 8, etc.

You can compensate for that by never returning many numbers, but that violates #2.

2) Generate numbers that meet the above definition AND so that any particular floating number is no more likely to be picked than any other.

But you can't simultaneously have #1 and #2.

#2 = #1 AND MORE and NOT (#1 AND #2) imply that #2 is impossible independent of #1 or MORE.

You most certainly can.

If you follow the github links in that blog post, you'll see in arbitrary_range.c that that's exactly what I do in "rd_positive". Pick the highest number in the range, calculate the distance to next lower representable number, only return multiples of that (bounded by the total number of those distances within the range). For IEEE 754 doubles, unless your range can't represent that many doubles at all it makes you end up with somewhere between 52 and 53 bits of entropy per number generated.

But doesn't that mean that you decide to split he set of represent able numbers in gate range in those that your code can return and those that it will never return?

Yes, the ones that you decide to return all have the same probability. The ones that you decide to never return have Rob ability zero. That is not "so that any particular floating number is no more likely to be picked than any other."

It's referring to the numbers that could be picked at all.

I made a mistake in my comment. It should say that any particular floating number is no more likely to be picked than any other with the same exponent. (I think that's right? So basically any representable number could be chosen, but ones in 'sparser' rangers will be more likely to be picked than a single number in a 'denser' range)

As you claim on your page github page: "To be useful this would require vastly more testing and independent verificaiton from someone who actually understands set theory and ieee 754 floating point numbers. I can confidently say that I know neither, but as the industry standard on stackoverflow and the c++ libraries seems to indicate, something something something one-eyed man and the land of the blind."

I believe your confusion comes exactly from not really understanding how floating point works on the low level. You really need that before you can write the meaningful numerical code.

Moreover, the most effective generation of the range of floating point numbers has to "know" the internal representation of floats. Specifically, the papers and the implementations of random number generators from Prof Makoto Matsumoto use that knowledge, and he specifically claims that he didn't invent the approach, but that it was known at least since 1997, crediting Agner Fog.

I'm sorry, but if you mean Makoto Matsumoto who seems to be the name mentioned in the SFMT library, that library is bogus. It generates numbers in [1,2) and then maps that to [0,1) by subtracting 1. This can be trivially shown to lose information, in fact I added a little hack to my github repo two days ago showing why it's bogus when that specific library was mentioned to me (again). Unless you want to lose one bit of entropy in the first conversion every number gets. And it still doesn't get us anywhere close to arbitrary ranges (which is the hard problem).

I also just looked at stuff that Agner Fog has published, thank you for mentioning him. Same thing. Maybe my understanding of floating point is bad, but his is worse. One generator uses the classic rng()/((double)range_of_rng) approach which is the whole reason why the C++ standard generators aren't uniform (because it will only work you use a non-standard rounding mode for floating points and still doesn't get us anywhere close to arbitrary ranges which are a completely different problem). Another uses [1,2) - 1.0 which loses information. And the third one is "combines two different random number generators for improved randomness", which I guess sounds good, how does he do it? He generates two [0,1) numbers, adds them together then subtracts 1.0 when the result is bigger than 1.0 which both loses precision and compounds the rounding problems.

Literally the first thing I did when trying to figure this out months ago was to randomize the bits in the mantissa with a 0 exponent to get a [1,2) number then subtracted 1.0 from it, then it took me 10 minutes to test that this doesn't actually work.

So, maybe my understanding of floating point numbers is bad, but the authorities you cite are worse.

The fast generators as described by Matsumoto and Fog generate the numbers this way:

"If the 52-bit fraction part is uniformly randomly chosen, then the represented number is uniformly randomly distributed over [1, 2) with 52-bit precision."

Then for [0,1) 1 is subtracted.

I guess you can be confused because you don't understand the binary floating point? Specifically you can have wrong assumption about the numbers that are to be generated by the given algorithms. Let me explain you on the smaller example: let's generate the floats instead of doubles. If all possible input integer values were generated by Matsumoto's algorithm i e [0..8388607] and then from them generated the floats Xi in the range [0..1) then multiplying Xi with the constant 8388608f should get the results: 0f, 1f, 2f, ... 8388605f, 8388606f, 8388607f. And that's what is produced. Doubles have more bits but function the same: the constant which makes the tests easy is 4503599627370496.0 and the last number from the algorithm after the multiplication with that constant should be 4503599627370495, and the previous 4503599627370494 etc. The first three should of course be 0, 1 and 2.

What does your "better" algorithm produce under the given assumptions? Or do you actually have some other assumptions? If you do, it doesn't mean their algorithms are wrong, these do what they are designed to do and obviously correctly.

I suggest you you test your algorithm on floats, using absolutely all inputs and calculating all the outputs. Then you can write here what you achieved, that is, if you started from other goals, did you really reach them for all the inputs? Note that you don't need any random number generator to prove you reached your goal, just show that for given input of all ascending integers you get the desired output -- every result should also appear only once (I see your "tests" on github aren't made so).

But please don't write about C++ or use some notation from the latest C++ standards. Just write the inputs, outputs and describe your algorithms mathematically.

Got pointed back here that there was another response. Sure. If you subtract 1 from [1,2) you will get 52 bit precision (or 23 for floats as your magical constant indicates). This means that the number just before 1.0 will not be picked, then the number before that will be, the next one won't be, etc. You don't pick every other number in the range where you should have problems covering every number without introducing any irregularity to the result. So what? Map your random numbers to the set [0,1,2,3] divide by 4 and you also get a uniformly distributed number in [0,1).

With my method we'd pick the magical constant 16777216.0f and we'd also get 0, 1, 2, ... 16777216.0f is bigger than 8388608f, therefore it's better.

This [1,2)-1.0 method has 0 probability of picking the first number less than 1.0, 2^-52 probability to pick the next one, 0 probability to pick the next one, 2^-52 probability of picking the next one, etc. My algorithm has 2^-53 probability of picking the first number less than 1.0, 2^-53 probability of picking the next one, etc. Down to .5 I can pick every representable number. After that there is a choice. I chose to have the same probability of picking every smaller number (including 0.0) as the probability of picking numbers in the upper half of the range, so I have holes that I can't pick, but at least in the upper half of the range I can pick every number I can represent.

That's why I think that [1,2)-1.0 is flawed. Because it adds an arbitrary limitation to its output, without any good reasons except that it was convenient to implement and they didn't have to think hard.

> With my method we'd pick the magical constant 16777216.0f and we'd also get 0, 1, 2, ... 16777216.0f is bigger than 8388608f, therefore it's better.

The constants I've mentioned (23 bits of input and 8388608.0f) aren't "magic" but are the essence of floats when generating the uniform distribution for floats in the range 0..1.

That's what I want you to actually try in your code: to write it with the inputs of 24 bit integers 0..16777215 and to demonstrate that when you produce the float values and multiply these with the 16777216.0f you get every number from the 0..16777215.0f exactly once. Unless you show me that result (repeatable by anybody) I'll claim that your code doesn't work what it claims to have had achieved.

So once again, to be really clear: don't input any "random" bits. Input the integer numbers 0..16777215 (always only the lower 24 bits of them!) and demonstrate that your code works for floats (as explained above, every resulting float multiplied by the 16777216.0f must then be in the set of 0..16777215.0f and happen only once!). Only then I can believe you the equivalent solution works for doubles too. Before that, I don't believe you really understand what you want to do and what you're actually producing.

Do see my answer here: https://news.ycombinator.com/item?id=9370844

The textual description (of the Nov 2014 draft I'm reading) also says it returns values x, where a <= x < b, while stating that it requires that a <= b. That's a problem if a == b.

Reply on news.ycombinator.com