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

> 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?

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

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".

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.

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 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.

"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.

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.

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.

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.

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.

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

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.

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."

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.