Add efficient method for precise integer to floating point conversion

[joshbainbridge]

Conversion from an integer value across the full range of representable values, to floating point values within the range of [0, 1) is a key part to QMC algorithms, as most calculations are done using integer arithmetic, but the resulting output often needs to be floating point.

The current method employed uses a single division, and relies heavily on the guarantees of IEEE 754. This states that the resulting floating point number in the [0, 1) range will be highly accurate. However, the operation will both round up, and down, and so results in uneven probabilities across the final distribution of values within the representable range.

Rounding up also means a value of 1 itself can be the result. Due to this, we add a max operation to clamp all values to the last representable number before 1, which also adds bias.

A fully accurate mapping can be implemented, but this is computationally much more expensive than a simple division, and often cost more than the creation of the QMC point. Such an implementation was tested in this commit: a97ad21

Add WIP implementation for float mapping using clz

Mapping 32 bit unsigned integers to floating point representation could
be improved as the current method has an undesirable rounding behaviour.

This method using a clz instruction based on the algorithm in Section
2.1 of 'Quasi-Monte Carlo Algorithms (not only) for Graphics Software'
directly maps the unsigned integer bits to the floating point exponent
and mantissa with minimal loss of information.

Testing of this method did show a considerable computational cost when
drawing samples. With some sampler implementations reaching almost 2x
time increase when benchmarked.

Implementing a method with the same output as the above commit, but with at least the same computational cost of the current approach using division would be the best outcome.

[toxieainc]

Hi! Author of the mentioned paper here.. I spent a lot of time thinking back then if this could be further optimized, but the approach in the paper is already so minimalistic that the only thing that might be possible, is to avoid the special case handling (0 & 1) by some trickery, but also unclear how exactly that could happen.

Note though that also the current implementation that maps the full 32bits via a single multiplication to [0..1) is debatable. I would recommend throwing away the lowest 8bits instead, as then all remaining bits can perfectly fit into a fp32. This way you'll get a perfectly even distribution.

[fpsunflower]

+1 to what @toxieainc said.

I would also advocate for float(bits >> 8) * 0x1p-24f as the safest way to remap from ints to float. In fact, this agrees with the "fancy" method above for ~2/3 of inputs due to how floats are spaced on [0,1).

As I mentioned in the last meeting, even a simple operation like 1-x can immediately "destroy" any excess precision near the origin, so I'm not sure there are any major advantages to spending extra cycles (especially on modern machines where the int-to-float cast and float multiply are no longer that expensive).

It could be worthwhile to come up with a few utility functions that take a uint32 and return a float with particular distributions. For example, if you need a float on the interval [-1,1) you can gain an extra bit of precision "for free". Same for other common transforms like -log(1-x), cos(2pi*x), etc ... That being said, I've been using the version with uniform spacing for a while and haven't noticed any particular need to go beyond it.

If more precision is needed one can simply upgrade to double precision: double(bits) * 0x1p-32 will give you perfectly uniform spacing for the full 32-bits. Chances are that if you care about the difference between 2^-24 and 2^-32 spacing in your random numbers, you will probably need doubles for the difference to carry over to the thing you are using the random number for.

[joshbainbridge]

Hey @toxieainc, great paper! And I'd agree that the current approach is not ideal. At Framestore we did originally use the method of throwing away the lowest 8 bits. The uniform precision is a nice property.

We later switched to a division, similar to PBRT. I can't recall the specifics, but there was a pathological production case where the lack of precision in the lower floating point representation caused us some headaches (most likely higher level volume sampling code). I can check with Manuel Gamito who might know more.

Thanks @fpsunflower, your point that any advantage from a fancier method is unlikely to make a difference, in most practical cases, is fair. Perhaps we could revisit dropping the lower 8 bits. This could also provide more opportunity for optimization in the [0, 1) case.

There might be one other option. I know @mr-matthew-jones was looking into this, and had a promising idea. It sounded like we could get all the precision, with round down, and as few as a bit op and a multiply. I'll reach out to see if he can share details before we make a call.

However, the lack of uniform precision across the [0, 1) range would still remain, so we may still choose to go either way depending on the consensus on that point.

[fpsunflower]

That's a very elegant approach! Simple enough to write on one line, and it matches the "fancy" method exactly!

Only one or two extra operations (depending on the instruction set) compared to the uniform version, so no real argument performance wise.

I think the only discussion item left is if there is any benefit to the "biased" extra precision.

Seeing this approach also make me wonder what can be done for generating floats on other common domains.

Since we are on the topic - I'll mention this approach as well:

float u01_naive_approx(uint32_t bits) {
    return float(bits) * 2.3283063e-10f;
}

This generates a float on [0,1) by multiplying by a constant nudged slightly away from 2^-32. such that the top range still rounds away from 1.0. Its still essentially uniformly spaced. It matches the "fancy" method for ~49% of inputs.