Floating point precision nuances.

Question

I found this code in NVIDIA's CUDA SDK samples.

void computeGold( float* reference, float* idata, const unsigned int len)
{
    reference[0] = 0;
    double total_sum = 0;
    unsigned int i;
    for( i = 1; i < len; ++i)
    {
        total_sum += idata[i-1];
        reference[i] = idata[i-1] + reference[i-1];
    }
    // Here it should be okay to use != because we have integer values
    // in a range where float can be exactly represented
    if (total_sum != reference[i-1])
        printf("Warning: exceeding single-precision accuracy.  Scan will be inaccurate.\n");
}
//(C) Nvidia Corp

Can somebody please tell me a case where the warning would be printed, and most importantly, why.

Answer 1

Normally, you cannot sum many floating point numbers.

Eventually the sum becomes different order of magnitude than every new added number, so precision is lost. For example, in case of float, adding a million numbers of the same order of magnitude gives the same result as ten million, because by the time it's done, every new number added doesn't change anything.

There're algorithms around this which involve a couple of multiplications for every added number (indeed, just to sum numbers properly). Yeah, floating point is tricky.

See http://floating-point-gui.de/

Answer 2

A float normally has a range of something like +/- 1e38, but a precision of only about 5 or 6 digits. This means, for example, that something like 12345678 can be stored, but it'll only be stored with ~6 digits of precision, so you'll get the warning.

Answer 3

The function is written with a certain range of input data in mind. If that input data expectation isn't met, the warning will print:

#include <stdio.h>
#define COUNT_OF(x) (sizeof(x)/sizeof(0[(x)]))

void computeGold( float* reference, float* idata, const unsigned int len)
{
    double total_sum = 0;
    unsigned int i;

    reference[0] = 0;

    for( i = 1; i < len; ++i)
    {
        total_sum += idata[i-1];
        reference[i] = idata[i-1] + reference[i-1];
    }
    // Here it should be okay to use != because we have integer values
    // in a range where float can be exactly represented
    if (total_sum != reference[i-1])
        printf("Warning: exceeding single-precision accuracy.  Scan will be inaccurate.\n");
}
//(C) Nvidia Corp


float data[] = {
    1.0,
    2.0,
    3.0,
    4.0,
    5.0
};

float data2[] = {
    123456.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    999999.0,
    123456.0
};

float ref[COUNT_OF(data2)] = {0.0};

int main()
{
    computeGold( ref, data, COUNT_OF(data));
    computeGold( ref, data2, COUNT_OF(data2));
    return 0;
}