Taming floating-point sums
orlp.net
orlp.net
Rather than striving for exact values for every computation, we make it so that the expected value is exact. For example, suppose our number system was restricted to the integers and we have incoming value 1.1, with stochastic rounding we'd round down to 1 90% of the time and 2 10% of the time giving an expected value of 1.1!
Further reading:
It is used all the time however in audio and image processing, where it is called dithering.
Also, sometimes you’re not limited by processing speed, but by the destination data structure (eg quantized to integers).
Your code that runs every frame to update the graphics logically wants to do something like this:
X += Vx
Y += Vy
To allow for velocities that aren't an integer number of pixels per frame and that are slower than 1 pixel per frame you'd want to actually store X in a fixed point format, say (Xi, Xf) where Xi is the integer part of the X position, and Xf is the fractional part times 256, so X = Xi + Xf/256. Similarly for Y. The object only actually moves on the screen when Xi or Yi changes.
Similarly velocity would also be in that format: Vx = Vxi + Vxf/256, and similar for Vy.
With that, the position update in your loop would be something like this:
Xf += Vxf
if that wrapped
Xi += 1
Xi += Vxi
For each object you end up needing 8 bytes (1 byte for each of Xi, Xf, Yi, Yf, Vxi, Vxf, Vyi, Vyf). That doesn't sound like much but the Intellivision only had something like 240 bytes in the console available. (If you couldn't get your RAM requirements down to that it was possible to have extra RAM in the cartridge but that would raise the cost).
So someone figured out that you didn't actually need to store Xf and Yf. Just generate then at random as needed! The loop then becomes something like this:
rb = random_unsigned_byte()
if Vxf + rb wraps
Xi += 1
Xi += Vxi
That turns out to work quite reasonably. Essentially it is interpreting Vxf as meaning that the object has a Vxf/256 chance of crossing a pixel boundary on a given frame. Thinking of it that way then suggests getting rid of the addition of the random byte with a wrap check and just doing a compare instead:
if Vxf > random_unsigned_byte()
Xi += 1
Xi += Vxi
Net result: we've cut the RAM for storing position and velocity from 8 bytes per object to 6, at the cost of needing to generate 2 random bytes per object per frame.
But you are right--a game might not need 16 bits for velocity. Suppose a game needed a maximum of just under 4 pixels/frame X velocity and needed a minimum of 1/64 pixel/frame.
It could have stored velocity in one byte like this FFFFFFII where the Fs are the bits of the fractional part of the velocity in units of 1/64 pixels/frame and the Is are the integer part in pixels/frame. By putting the fractional part in the upper bits and the integer part in the lower bits no shifting will be required--just AND and OR.
Then updating X would be something like this:
X += Vx % 0x0003
if Vx > (random_unsigned_byte() | 0x0003)
X += 1
I don't think anyone would notice. For two objects in a race all the way across the screen with the same Vx, with one including the OR and one not, my back of the envelope calculation says that the speed boost would only get that one to the finish line less than half a millisecond earlier.
[1] http://www.bitsavers.org/components/gi/CP1600/CP-1600_Microp...
[2] https://en.wikipedia.org/wiki/General_Instrument_CP1600#Inst...
Oops! That should be & not %.
[1]: https://www.allaboutcircuits.com/technical-articles/what-is-...
[2]: https://www.analog.com/en/resources/analog-dialogue/articles...
See my paper at https://arxiv.org/abs/1505.05571 (repository at https://gitlab.com/radfordneal/xsum), which is used in the Julia xsum package (https://docs.juliahub.com/General/Xsum/stable/). There's also recent work by Marko Lange at https://web.archive.org/web/20220616031105id_/https://dl.acm...
I think these methods are probably faster than the exact methods you benchmark.
A new copy costs $112 on Amazon: https://www.amazon.com/Handbook-Floating-Point-Arithmetic-Je...
That said, I should have included it, because it is on the Pareto frontier. On my machine it is ~28.9 GB/s, with 0 error (note that this doesn't mean it always produces a correctly-rounded sum, just in this benchmark test input).
I'll add it to the article tomorrow or the day after.
I'd be curious to see the performance of hypothetical block_pairwise_addfast and block_kahan_addfast.
This is exactly what you are loosing when using fadd_fast: fine control over the errors terms of floating point operations that do matter in a lot of cases.
An other thing you are loosing is reproducibility: depending on the machine, compiler version, etc, your program will compute differently and for example may switch from linear to quadratic errors terms when you recompile. It could be the difference between a numerical algorithm converging or not, this kind of things.
The sum part can be done with AVX-512 [1] if N=8, and N=8 is probably enough for a lot of real world summations. The conversion back to floating point at the end is only done once so not very costly.
The conversion to fixed point is probably the worst part. Is there a faster way to do it than a loop that extracts one 64-bit integer per iteration? If not we could convert several input numbers at a time with SIMD. But it would be nifty to be able to convert one double to N integers with SIMD instead of a loop.
[1] http://www.numberworld.org/y-cruncher/internals/addition.htm...
I wonder what methods give good results for this sort of input though.
Kind of unfortunate that that syntax is so trivial to misread, but it is what it is.
The code from the blog post is all on Github, feel free to try and add/benchmark it yourself: https://github.com/orlp/sum-bench/.
a/b + c/d -> ( ad + cb ) / bd
usually in computers there is someone that has already "tried it" (example - trinary computers) but i haven't found stories of rational CPU experiments. surely some mad scientist tried it at some point. hoping HN commenters can point me to some history.