Add uniform function producing a random float in (0.0, 1.0)

[xavierleroy]

This is a follow-up to the discussion that started here: ocaml/RFCs#28 (comment)

This PR adds a uniform function that returns a random FP number in (0.0, 1.0), i.e. guaranteed to be neither 0.0 nor 1.0. It documents precisely the shape of the FP numbers returned.

The float function is implemented and documented in terms of uniform.

The old float_32 implementation was removed as not really useful and as producing differently-rounded results than the float_64 implementation.

Cc: @jhjourdan

[jhjourdan]

This implementation certainly works, but:

• Why masking the lowest 11 bits when we could in fact take them into account (and multiply by 0x1p-64 instead)?
• Testing n = 0 and branching can be a bit costly. In stead, the same behaviour could be obtained by adding a small number after the multiplication (0x1p-54 or 0x1p-65, depending on what you choose for the first remark). Sure, then you need to change the documentation below...

[xavierleroy]

"take [all the bits] into account (and multiply by 0x1p-64 instead)": then you get FP rounding during the int64 -> float conversion. The rounding can produce 0.0 and 1.0, and seems to introduce some (very subtle) bias in the output, as discussed in Generating Random Floating-Point Numbers by Dividing Integers: a Case Study. With the current code, we know exactly which FP values we can get, and we're sure they all have the same probability.

"Testing n = 0 and branching can be a bit costly". Really? This is one integer comparison plus a branch that is perfectly predicted (as n is almost never 0).

"adding a small number after the multiplication": the resulting rounding makes me nervous (again).

[xavierleroy]

Just for fun: I tweaked the n = 0 test so that the conditional branch is statically predicted as not taken.

[jhjourdan]

Right. Adding 0x1p-54 here would make producing 1. possible because of the round-to-even rule. And adding a slightly smaller number would introduce a -0x1p55 biais, because all values > 0.5 would be rounded down. While this is a very small biais (probably impossible to notice), I then agree that what you propose is probably the best.

My concern about branching was a remainder of memories of what's happening when writing the same code in C : then, the introduced branching would be a problem, since it would prevent some SIMD optimizations. But Ocamlopt does not support these optimizations anyway, and the branch will very likely be predicted.

[xavierleroy]

Right, I forgot about your adventures with SIMD :-)

I take it that you find the name uniform acceptable for this function ? The only alternative I could think of is float1, but it's ugly.

[jhjourdan]

I agree that uniform is not such a good name, because it does not follow the convention of the rest of the file, where the name follow, in some sense, the type of the return value, while uniform describes the distribution.

One possibility would be to not provide publicly this function, but still document the new guarantee. The downside is the additional multiplication by the constant 1., but:

• This is "just" a float multiplication, which is not very costly
• We could hope (hint, hint !) that the OCaml compiler would remove the multiplication if one of the multiplicand is 1., which would be clear after inlining. (Correct me if I'm wrong but removing multiplication by 1.0 is one of the only correct simplification of float computations.)

[xavierleroy]

I agree the additional multiplication by 1.0 is not costly, so it's OK if users continue to call [float g 1.0] instead of [uniform g]. On the other hand, I find it makes the documentation clearer. By "it" I mean "having the uniform function exported and documented separately". So let's do that.