C++ 11 std::async for Fast Float Format Finding

After a recent post on float precision there was some debate about round-tripping of floats. My claim was that if you print a float with printf(“%1.8e”, f); and then scan it back in then you are guaranteed to get back exactly the same float, but there was some skepticism.

Bold claims need careful proof, and this post shows how C++ 11 can be leveraged to verify this statement faster and more easily than with old-style C++.

The best way to prove my claim is to write some test code. The code to see whether you can round-trip floats to decimal and back without ‘drift’ is pretty simple. You just have to loop through all two billion positive floats (you can use this code as a starting point) and for each float, sprintf it to ASCII and then sscanf it back to a float. If you always get back exactly the same number then the hypothesis is proven.

However this is time-consuming. On my laptop it took about fifty minutes to run. Normally multi-threading it would be too much work for a throwaway task like this, but with a fresh install of VS 11 beta on my laptop I decided to give std::async a try. It worked really well and gave me a 3.6x speedup.

In the code below you can see calls to futures.push_back where I queue up asynchronous calls to FindDriftingFloats(start, end). A bit later I loop through those futures, retrieving the results. This code will automatically parallelize across up to 255 cores, and increasing the level of parallelism is just a matter of doing smaller batches by incrementing start.i by a smaller amount.

The cool thing is that, for the first time, it is possible to write portable multithreaded C++ code.

Here’s the code – focus on the futures vector.

#include <future>
#include <vector>
#include <stdio.h>

/* See
Tricks With the Floating-Point Format

for the potential portability problems with the union and bit-fields below.
*/
union Float_t
{
    Float_t(float num = 0.0f) : f(num) {}
    // Portable extraction of components.
    bool Negative() const { return (i >> 31) != 0; }
    int32_t RawMantissa() const { return i & ((1 << 23) - 1); }
    int32_t RawExponent() const { return (i >> 23) & 0xFF; }

    int32_t i;
    float f;
#ifdef _DEBUG
    struct
    {   // Bitfields for exploration. Do not use in production code.
        uint32_t mantissa : 23;
        uint32_t exponent : 8;
        uint32_t sign : 1;
    } parts;
#endif
};

// Print each float in the range and verify that it round-trips.
// This function is designed to be called on multiple threads.
// Therefore no global variables should be modified.
int FindDriftingFloats(Float_t start, Float_t end)
{
    int failures = 0;
    // These print statements may occur out of order due to this
    // code being run on multiple threads.
    printf("Scan %1.8ef to %1.8ef (%x to %x)\n",
                start.f, end.f, start.i, end.i);

    while (start.f < end.f)
    {
        char buffer[30];
        sprintf_s(buffer, "%1.8e", start.f);
        float fRead;
        int count = sscanf_s(buffer, "%f", &fRead);
        if (count != 1 || fRead != start.f)
        {
            printf("Conversion failure on %1.8e (%08x)!\n",
                        start.f, start.i);
            ++failures;
        }

        start.i += 1;
    }

    return failures;
}

int main(int argc, char* argv[])
{
    auto startTime = clock();
    std::vector<std::future<int>> futures;

    // This represents infinity. We stop here, without processing
    // this number.
    const uint32_t maxRep = 255 << 23;
    Float_t start(0);
    while (start.i < maxRep)
    {
        Float_t end;
        // Increment by any amount. 1<<23 means that each exponent
        // is processed as one batch, but any positive integer will do.
        end.i = start.i + (1 << 23);
        if (end.i >= maxRep)
            end.i = maxRep;

        // Call FindDriftingFloats(start, end) on another thread.
        futures.push_back(std::async(std::launch::async,
                FindDriftingFloats, start, end));
        // Advance to the next range.
        start = end;
    }

    // Record the results. The results may not complete in order, but we
    // wait for them in order.
    int totalFailures = 0;
    for (size_t i = 0; i < futures.size(); ++i)
    {
        int failures = futures[i].get();
        totalFailures += failures;
        // Cast size_t to a type that can be portably printed. unsigned
        // long long would be better if the values could actually be large.
        printf("%d failures found in range %u.\n", failures, (unsigned)i);
    }

    printf("%d floats failed to round-trip to decimal.\n", totalFailures);
    auto elapsedTime = clock() - startTime;
    printf("Took %1.3fs\n", elapsedTime / double(CLOCKS_PER_SEC));

    return 0;
}

Typical output is shown here:

Scan 0.00000000e+000f to 1.17549435e-038f (0 to 800000)
Scan 4.70197740e-038f to 9.40395481e-038f (1800000 to 2000000)
Scan 1.88079096e-037f to 3.76158192e-037f (2800000 to 3000000)
Scan 3.76158192e-037f to 7.52316385e-037f (3000000 to 3800000)
Scan 1.17549435e-038f to 2.35098870e-038f (800000 to 1000000)
Scan 2.35098870e-038f to 4.70197740e-038f (1000000 to 1800000)
Scan 9.40395481e-038f to 1.88079096e-037f (2000000 to 2800000)
Scan 7.52316385e-037f to 1.50463277e-036f (3800000 to 4000000)
Scan 1.50463277e-036f to 3.00926554e-036f (4000000 to 4800000)
Scan 3.00926554e-036f to 6.01853108e-036f (4800000 to 5000000)
Scan 6.01853108e-036f to 1.20370622e-035f (5000000 to 5800000)
Scan 1.20370622e-035f to 2.40741243e-035f (5800000 to 6000000)
Scan 2.40741243e-035f to 4.81482486e-035f (6000000 to 6800000)
Scan 4.81482486e-035f to 9.62964972e-035f (6800000 to 7000000)
Scan 9.62964972e-035f to 1.92592994e-034f (7000000 to 7800000)
0 failures found in range 0.
Scan 1.92592994e-034f to 3.85185989e-034f (7800000 to 8000000)
0 failures found in range 1.
0 failures found in range 2.
0 failures found in range 3.

…

Scan 8.50705917e+037f to 1.70141183e+038f (7e800000 to 7f000000)
0 failures found in range 246.
Scan 1.70141183e+038f to 1.#INF0000e+000f (7f000000 to 7f800000)
0 failures found in range 247.
0 failures found in range 248.
0 failures found in range 249.
0 failures found in range 250.
0 failures found in range 251.
0 failures found in range 252.
0 failures found in range 253.
0 failures found in range 254.
0 floats failed to round-trip to decimal.
Took 830.424s

Note that the “Scan …” printouts don’t always appear in order. That is because they are being done on separate threads and their precise timing is at the whim of the Windows scheduler.

Note also that fifteen jobs were started before the first result was retrieved. That may indicate that the first job took a particularly long time (results are retrieved in order), or it may mean that the main thread was starved for CPU time and was briefly unable to retrieve results.

My laptop has four hyperthreaded cores. The fact that I only got a speedup of 3.6x is interesting, and probably has to do with the single-threaded CPU frequency being higher due to TurboBoost, together with lackluster improvements from hyperthreading. It happens. I did verify that all eight hardware threads were pegged the entire time.

That’s it. Multi-threading done easily enough that even on run-once jobs like this I can multithread my tests in less time than the 36 minutes that multi-threading will save me. Win-win.

The total time to run to run the multithreaded test was under 14 minutes. On a desktop machine with more and faster cores it could easily be done in half that time or less – not bad for running a test on two billion floats.