# Why different sums for these two functions for large sets of data?

Discussion in 'C++' started by Eric Lilja, May 24, 2005.

1. ### Eric LiljaGuest

Hello, consider the following two functions:

/* function foo() */
void foo()
{
float y = 0.0f;
float sum = 0.0f;

for(int i = 0; i < num; ++i)
{
for(int j = 0; j < num; ++j)
{
const float& x = myDataX;

y = myDataY[j];

sum += x + y;
}
}

std::cout << "sum = " << sum << std::endl;
}

/* function bar() */
void bar()
{
float sum = 0.0f;
float y_sum = 0.0f;

for(int y = 0; y < num; ++y)
y_sum += myDataY[y];

for(int x = 0; x < num; ++x)
{
float x_times_num = myDataX[x] * num;

float value = (x_times_num + y_sum);
sum += value;
}

std::cout << "sum = " << sum << std::endl;
}

If myDataX and myDataY are small arrays, both functions produce the correct
sum and bar() seems to be a bit faster (I removed code for timing when
posting). But for large sets of data, bar() produces a too large sum but is
very much faster. Why? Where's the bug in bar()? I have a feeling it has to
do with the fact that I'm using floats and floats don't always behave as one
might expect unless you understand in detail how they're represented
internally. I want to fix bar() so it produces the correct sum for large
data sets, and, hopefully still doing it faster than foo(). I am compiling
with debug info and will all optimizations disabled.

If you see any bugs, have suggestions on how to improve the speed of my
functions, please reply. I need to calculate the sum calculated by foo() and
bar() as fast as possible.

Example data set that works for both:
const int num = 4;
myDataX = new float[num];
myDataY = new float[num];
myDataX = 1;
myDataX = 2;
myDataX = 3;
myDataX = 4;

myDataY = 5;
myDataY = 6;
myDataY = 7;
myDataY = 8;

Example data set (large) where bar() produces an incorrect sum:
const int num = 20000;
myDataX = new float[num];
myDataY = new float[num];

std::fill(myDataX, &myDataX[num-1], 0.1f);
std::fill(myDataY, &myDataY[num-1], 0.1f);

These variables are static globals in my test program.

Thanks for any replies

/ E

Eric Lilja, May 24, 2005

2. ### HowardGuest

Is that correct? I get totally bogus results when trying to run this code.
Try setting num to 4 for this case, and see what I mean. What do you get?
I get large negative values from both. I think the fill statement isn't
right. (Although I don't know how to fix it, since I've never used
std::fill with an array.)

A few suggestions:

1) why not use a vector, instead of an array?
2) why not use double instead of float, (since it has better precision,
which will help in any large summing of float values)?
3) isn't that formula equal to just: num times (the sum of the Xs plus the
sum of the Ys)?

-Howard

Howard, May 24, 2005

3. ### Donovan RebbechiGuest

Are you sure you've got this the right way around ? Your code looks very
inefficient. If you're trying to compute the sum

Sum(i=0..num-1,j=0..num-1) (x_i + y_j)

then you're doing this very inefficiently. Note that this sum exands as

Sum(i,j) x_i + Sum(i,j) y_j

Now because x_i is independent of j, and x_j is independent of x_i, you can
collapse a subscript of each sum, and get:

size * (Sum(i) x_i +sum(j) x_j)

Sum

The resulting code looks like this:

float bar ( float*x, float* y, int size )
{
double sumx = 0, sumy = 0;
for (int i = 0; i < size; ++i)
{
sumx += *x++;
sumy += *y++;
}
return size*(sumx + sumy);
}

or alternatively, as an STL one liner:
float bar ( float*x, float* y, int size )
{
return size * (std::accumulate(x,x+size,0.f) +
std::accumulate(y,y+size,0.f)
}

This should be (and is, according to my test with your data) more accurate
and one multiplication, compared to yours : 2*size*size floating point
additions. When I ran the test (with my version of bar), I got these results:
foo: 4.1943e+06
bar: 8e+07

even with only 200 elements in the array, I get:
foo: 8003.11
bar: 8000

Cheers,

Donovan Rebbechi, May 24, 2005
4. ### Eric LiljaGuest

[my OP snipped]
I tried that version and compared it with the original version that I was
given (that I've been trying to improve upon and that I didn't show in my
OP). Here's a complete program with output for the two varaints. As you can
see the results are equal for the small data set but way off for the large
data set:

#include <algorithm>
#include <cstddef>
#include <cstring>
#include <iomanip>
#include <iostream>

static double
un_optimized(const float *myDataX, const float *myDataY,
std::size_t num)
{
unsigned int i, j;
float sum = 0;

for( i = 0; i < num; i++ )
{
for( j = 0; j < num; j++ )
{
float x = myDataX;
float y = myDataY[j];
float value = x+y;
sum+=value;
}
}

return sum;
}

static double
optimized(const float *myDataX, const float *myDataY, std::size_t num)
{
double sumx = 0, sumy = 0;

for(unsigned int i = 0; i < num; ++i)
{
sumx += *myDataX++;
sumy += *myDataY++;
}

return (num*(sumx + sumy));
}

int
main()
{
std::size_t num = 20000;
float *myDataX = new float[num];
float *myDataY = new float[num];

std::fill(myDataX, &myDataX[num], 0.1f);
std::fill(myDataY, &myDataY[num], 0.1f);
//std::memset(myDataX, 1, num);
//std::memset(myDataY, 1, num);

std::cout << std::fixed << std::showpoint << std::setprecision(8);

std::cout << "sum for un_optimized(): "
<< un_optimized(myDataX, myDataY, num) << std::endl;
std::cout << "sum for optmized(): "
<< optimized(myDataX, myDataY, num) << std::endl;

delete [] myDataX;
delete [] myDataY;

num = 4;
float myDataX2[] = {1,2,3,4};
float myDataY2[] = {5,6,7,8};

std::cout << "sum for un_optimized(): "
<< un_optimized(myDataX2, myDataY2, num) << std::endl;
std::cout << "sum for optmized(): "
<< optimized(myDataX2, myDataY2, num) << std::endl;

return 0;
}

Output:
sum for un_optimized(): 4194304.00000000
sum for optmized(): 80000001.19209290
sum for un_optimized(): 144.00000000
sum for optmized(): 144.00000000

/ Eric

Eric Lilja, May 24, 2005
5. ### Shezan BaigGuest

static double optimized(const float *myDataX,
const float *myDataY,
std::size_t num)
{
float x = *myDataX, y = *myDataY;
return num * num * x * y;
}

(Since you are relying on the fact that all elements in each array
contain the same value, you don't even need the loop).

Hope this helps,
-shez-

Shezan Baig, May 24, 2005
6. ### Larry I SmithGuest

That should be:

std::fill(myDataX, myDataX + num, 0.1f);
std::fill(myDataY, myDataY + num, 0.1f);

The 2nd arg needs to be an 'end()' iterator,
i.e. 'one past the end'.

Larry

Larry I Smith, May 24, 2005
7. ### Larry I SmithGuest

Eric Lilja wrote:
[snip]

Oops, the above MUST be 'double sum = 0;'
Otherwise it will overflow for large values of 'num'.

[snip]

Change 'sum' to 'double' in un_optimized() (see my comments above)
then this program outputs:

sum for un_optimized(): 80000001.19209290
sum for optmized(): 80000001.19209290
sum for un_optimized(): 144.00000000
sum for optmized(): 144.00000000

Regards,
Larry

Larry I Smith, May 24, 2005
8. ### Niels DybdahlGuest

floats usually have a precision of 6-7 digits. When your result becomes
larger than that compared to the values that you are adding, you loose the
values. This is a bigger problem in foo() because you only add x+y in each
loop. In bar() you add num*x+num*y in each loop, which is a larger value.
If you use double instead you will have around 15 digits precision and that
will solve the problem in both functions.

Niels Dybdahl

Niels Dybdahl, May 25, 2005
9. ### Karl Heinz BucheggerGuest

Until you know what you do and are willing and have the knowledge
to fight that beast:
* avoid 'float'

Karl Heinz Buchegger, May 25, 2005
10. ### Old WolfGuest

This is a bizarre way to write that algorithm, consider:

for(int i = 0; i < num; ++i)
for(int j = 0; j < num; ++j)
sum += myDataX + myDataY[j];
What makes you so sure that foo() is right and bar() is wrong?

Every floating point operation incurs a small precision error,
and foo does many many more operations than bar, so my guess
would be that bar's result is more correct.
1. Use the simpler version everyone's already suggested
2. Use 'double' for your intermediate results. That will probably
be faster than using 'float' , as an added bonus.
3. In your compiler, turn on 'fast floating point' option. This
will exponentially increase the speed of your program, but
it may violate ISO C requirements for what to do when precision
is lost (ie. you may get a slightly different answer, for
better or worse).

Why is 'double' maybe faster ? Many CPUs have builtin instructions
for arithmetic in 'double', and when you store the result in a
float, it has to perform a truncation operation (and if it has to
conform to ISO C, an even slower truncation operation).
Those fills are wrong, you are not filling the last member.
As well as the option suggested by Larry, consider this one:

std::fill_n(myDataX, num, 0.1f);

An even better option would be to use RAII:

std::vector<float> myDataX(num, 0.1f);

The correct answer to this calculation is 20000 *
(20000 * (0.1 + 0.1)) ie 80000000.
So it looks like bar() is indeed giving the more accurate