There are Only Four Billion Floats–So Test Them All! | Random ASCII – tech blog of Bruce Dawson

Random ASCII – tech blog of Bruce Dawson

Forecast for randomascii: programming, tech topics, with a chance of unicycling

Skip to content

Home<br>About

&larr; Bugs I Got Other Companies to Fix in 2013

Process Tree from an Xperf Trace &rarr;

There are Only Four Billion Floats–So Test Them All!

Posted on January 27, 2014 by brucedawson

A few months ago I saw a blog post touting fancy new SSE3 functions for implementing vector floor, ceil, and round functions. There was the inevitable proud proclaiming of impressive performance and correctness. However the ceil function gave the wrong answer for many numbers it was supposed to handle, including odd-ball numbers like ‘one’.

The floor and round functions were similarly flawed. The reddit discussion of these problems then discussed two other sets of vector math functions. Both of them were similarly buggy.

Fixed versions of some of these functions were produced, and they are greatly improved, but some of them still have bugs.

Floating-point math is hard, but testing these functions is trivial, and fast. Just do it.

The functions ceil, floor, and round are particularly easy to test because there are presumed-good CRT (C RunTime) functions that you can check them against. And, you can test every float bit-pattern (all four billion!) in about ninety seconds. It’s actually very easy. Just iterate through all four-billion (technically 2^32) bit patterns, call your test function, call your reference function, and make sure the results match. Properly comparing NaN and zero results takes a bit of care but it’s still not too bad.

Aside: floating-point math has a reputation for producing results that are unpredictably wrong. This reputation is then used to justify sloppiness, which then justifies the reputation. In fact IEEE floating-point math is designed to, whenever practical, give the best possible answer (correctly rounded), and functions that extend floating-point math should follow this pattern, and only deviate from it when it is clear that correctness is too expensive.

Later on I’ll show the implementation for my ExhaustiveTest function but for now here is the function declaration:

typedef float(*Transform)(float);

// Pass in a range of float representations to compare against.<br>// start and stop are inclusive. Pass in 0, 0xFFFFFFFF to scan all<br>// floats. The floats are iterated through by incrementing<br>// their integer representation.<br>void ExhaustiveTest(uint32_t start, uint32_t stop, Transform TestFunc,<br>Transform RefFunc, const char* desc)<br>Typical test code that uses ExhaustiveTest is shown below. In this case I am testing the original SSE 2 _mm_ceil_ps2 function that started the discussion, with a wrapper to translate between float and __m128. The function didn’t claim to handle floats outside of the range of 32-bit integers so I restricted the test range to just those numbers:

float old_mm_ceil_ps2(float f)<br>__m128 input = { f, 0, 0, 0 };<br>__m128 result = old_mm_ceil_ps2(input);<br>return result.m128_f32[0];

int main()<br>// This is the biggest number that can be represented in<br>// both float and int32_t. It’s 2^31-128.<br>Float_t maxfloatasint(2147483520.0f);<br>const uint32_t signBit = 0×80000000;<br>ExhaustiveTest(0, (uint32_t)maxfloatasint.i, old_mm_ceil_ps2, ceil,<br>"old _mm_ceil_ps2");<br>ExhaustiveTest(signBit, signBit | maxfloatasint.i, old_mm_ceil_ps2, ceil,<br>"old _mm_ceil_ps2");<br>Note that this code uses the Float_t type to get the integer representation of a particular float. I described Float_t years ago in Tricks With the Floating-Point Format .

How did the original functions do?

_mm_ceil_ps2 claimed to handle all numbers in the range of 32-bit integers, which is already ignoring about 38% of floating-point numbers. Even in that limited range it had 872,415,233 errors – that’s a 33% failure rate over the 2,650,800,128 floats it tried to handle. _mm_ceil_ps2 got the wrong answer for all numbers between 0.0 and FLT_EPSILON * 0.25, all odd numbers below 8,388,608, and a few other numbers. A fixed version was quickly produced after the errors were pointed out.

Another set of vector math functions that was discussed was DirectXMath. The 3.03 version of DirectXMath’s XMVectorCeiling claimed to handle all floats. However it failed on lots of tiny numbers, and on most odd numbers. In total there were 880,803,839 errors out of the 4,294,967,296 numbers (all floats) that it tried to handle. The one redeeming point for XMVectorCeiling is that these bugs have been known and fixed for a while, but you need the latest Windows SDK (comes with VS 2013) in order to get the fixed 3.06 version. And even the 3.06 version doesn’t entirely fix XMVectorRound.

The LiraNuna / glsl-sse2 family of functions were the final set of math functions that were mentioned. The LiraNuna ceil function claimed to handle all floats but it gave the wrong answer on 864,026,625 numbers. That’s better...