# Pairwise summation

This function adds up an array pairwise. It runs fast enough and appears cache-friendly; I've tested it for correct output. It's in C++ (farther down in the file there's a function that takes a vector), but looks like valid C as well (but I haven't tried it in C). How easy is it to understand?

double pairwisesum(double *a,unsigned n)
{
unsigned i,j,b;
double sums[32],sum=0;
for (i=0;i+7<n;i+=8)
{
b=i^(i+8);
if (b==8)
sums[3]=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
else
{
sums[3]+=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
for (j=4;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (;i<n;i++)
{
b=i^(i+1);
if (b==1)
sums[0]=a[i];
else
{
sums[0]+=a[i];
for (j=1;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (i=0;i<32;i++)
if ((n>>i)&1)
sum+=sums[i];
return sum;
}


For comparison, here is the previous version, which was slower because it kept making the cache miss, and also clobbered the array:

double pairwisesum(double *a,unsigned n)
// a is clobbered.
{
unsigned i,j;
if (n)
{
for (i=1;i<n;i*=2)
for (j=0;j+i<n;j+=2*i)
a[j]+=a[j+i];
return a[0];
}
else
return 0;
}

• Just out of curiosity, what kind of data are you dealing with where summing the array this way improves precision? – apilat Feb 22 at 22:00
• *A power series, where the terms increase, then decrease. I could sort them by absolute value, then start adding the small ones, but sorting would take more time. *Matrix multiplication. *In future, computing volumes by adding lots of samples. One of the surfaces may be a road, whose alignment has an Euler spiral (computed by the above power series), which cannot be expressed in closed form, so the volume has to be computed by sampling (which will be low-discrepancy for faster convergence). *A set of floating-point checksums computed from the elevations of hundreds of millions of points. – Pierre Abbat Feb 22 at 23:37

This function adds up an array pairwise.

Not sure what that means.

It runs fast enough and appears cache-friendly; I've tested it for correct output.

Good so it works.

It's in C++ (farther down in the file there's a function that takes a vector), but looks like valid C as well (but I haven't tried it in C).

No its not.
Its C. You have just compiled it with a C++ compiler.

Being C++ is not just about writing code that can be compiled by the C++ compiler. It is a style of writing code. This is a typical C style. It just happens to compile under C++ compiler but it's not C++.

This interface:

double pairwisesum(double *a,unsigned n)


Couple of issues: We don't pass around pointers in C++ (we can but as a language we decided it was a bad idea as it does not show ownership). So we tend to use other techniques (i.e references or smart pointers).

In this case its not just a pointer: it's an array. The trouble here is that the array has decayed so we can't actually tell that it is an array any more and just have to assume the caller has not done something stupid. In C++ we like to understand the object being passed so we can use it correctly.

In C++ (up to C++17) we would normally represent this as iterators:

 template<typename T>
double pairwisesum(I begin, I end)
{
for(I loop = begin; loop != end; ++loop) {


From C++20 we would probably represent this by a range:

 tempalte<typename R>
double pairwisesum(R range)
{
for (auto const& r: range) {


Note a Range basically supports begin() and end() and is interchangeable with a View. If I got that incorrect, then sorry: we are still getting used to these concepts as they are new to C++.

How easy is it to understand?

I have read this code a couple of times and have no idea what it is trying to achieve. So I don't think it is readable or understandable without some more explanation. Some comments on what was happening would have been good.

## Code Review

Declare variables:

• As close to to the point of use as possible.

• Declare one variable per line.

unsigned i,j,b;
double sums[32],sum=0;


What are you trying to achieve? Vertical space-saving is not a good style. Make it clear as possible to use.

I know a lot of people use i/j for loop variables (especially coming from C or FORTRAN). But it is objectively a bad choice. Once you start adding comments to your code searching for i and j in the code gives you way too many false hits when trying to find all the use cases. So unless your loop is literally one line obvious and the loop variable is only used in one place don't do it. Even if all the above is true, don't do it because the code will change over time and writing it expressively now will save time later.

You can declare your variable in the loop:

  for (i=0;i+7<n;i+=8)


Write like this:

  for (unsigned loop = 0; loop + 7 < mac; loop += 8) {


There is no precedence rules you are trying to get around here. Those extra braces are just adding noise to the code:
      sums[3]=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));


Why not:

      sums[3]=a[i]+a[i+1]+a[i+2]+a[i+3]+a[i+4]+a[i+5]+a[i+6]+a[i+7];


Why not use a standard algorithm:

      auto start = a + i;
auto end   = start + 8;
sums[3]    = std::accumulate(start, end, 0);


As it turns out the parentheses here is for a reason. Because we are dealing with doubles. The parentheses reduces the amount of error that creeps into the sum by making sure we do addition in a particular order.

I would separate this additional line out into its own function (with a very explanation-based name) so it documents why you are adding this up in a particular order. This will be important for future maintainers as they may not immediately understand why all the braces are there (and do what I did and just remove them to make the code easy to read).

sums[3] = sumValuesBecauseOfSomeReasonIHaveThatReducesErrors(a + i)

inline double sumValuesBecauseOfSomeReasonIHaveThatReducesErrors(double* a) {
return (((a[0]+a[1])+(a[2]+a[3]))
+((a[4]+a[5])+(a[6]+a[7])));
}


      for (j=4;b>>(j+1);j++)


What are you testing for in b >> (j+1)? Us a named function to make the meaning clear.

      for (unsigned jLoop; insideReachValid(b, jLoop); ++jLoop) {


As a side note: prefer ++j over j++. For integers there is no difference. But if you are looping with other types it can make a minor difference and the idea is that you want to be optimal whatever the loop type. A lot of C++ code is modified by simply changing the type of the objects and you don't want to change the type then go back and change how it used as well.

Always add the braces around '{} around sub blocks. Its not always clear with indentation (especially when it does badly) what is in the sub block. Make it explicit by using the braces.

      for (unsigned jLoop; insideRacheValud(b, j); ++j) {
sums[jLoop] += sums[jLoop - 1];
}


This should be your comment above the code.

// This is the original implementation.
// The code has been updated to prevent cache missing and optimize the
// the performance of the operation.
//
// This code is deliberately left here as a comment to provide a reference
// for generating unit tests and validation and to explain what the code
// is actually trying to achieve (as the optimized version is non-trivial to understand).
#if 0
double pairwisesum(double *a,unsigned n)
// a is clobbered.
{
unsigned i,j;
if (n)
{
for (i=1;i<n;i*=2)
for (j=0;j+i<n;j+=2*i)
a[j]+=a[j+i];
return a[0];
}
else
return 0;
}
#endif
`
• I'm not sure what you mean by "acoud", but the parentheses are not there for noise. They're there to make the summation pairwise. Without them, the summation would have more error. – Pierre Abbat Feb 22 at 19:29
• @PierreAbbat: Of course we are adding doubles so the order does matter. That would be a good comment to add to the code. – Martin York Feb 22 at 19:32
• So far I haven't seen a reason to go past C++17 (which introduced readers-writer locks, which I use). The original version is no longer in the code; I had to dig in git. I'll put it back in. – Pierre Abbat Feb 22 at 19:35