calculate prefix sum

Question

I have following code to accomplish prefix sum task:

  #include <iostream>
  #include<math.h>
  using namespace std;

  int Log(int n){
      int count=1;
      while (n!=0){
          n>>=1;
          count++;

      }
      return count;
  }
  int main(){
    int x[16]={39,21,20,50,13,18,2,33,49,39,47,15,30,47,24,1};
    int n=sizeof(x)/sizeof(int );
    for (int i=0;i<=(Log(n)-1);i++){
          for (int j=0;j<=n-1;j++){
              if (j>=(std::powf(2,i))){
                  int t=powf(2,i);
                  x[j]=x[j]+x[j-t];

              }
          }
     }
     for (int i=0;i<n;i++)
          cout<<x[i]<< "  ";

     return 0;
  } 

From this wikipedia page but i have got wrong result what is wrong? please help

Answer 1

I’m not sure what your code is supposed to do but implementing a prefix sum is actually pretty easy. For example, this calculates the (exclusive) prefix sum of an iterator range using an arbitrary operation:

template <typename It, typename F, typename T>
inline void prescan(It front, It back, F op, T const& id) {
    if (front == back) return;
    typename iterator_traits<It>::value_type accu = *front;
    *front++ = id;
    for (; front != back; ++front) {
        swap(*front, accu);
        accu = op(accu, *front);
    }
}

This follows the interface style of the C++ standard library algorithms.

To use it from your code, you could write the following:

prescan(x, x + n, std::plus<int>());

Are you trying to implement a parallel prefix sum? This only makes sense when you actually parallelize your code. As it stands, your code is executed sequentially and doesn’t gain anything from the more complex divide and conquer logic that you seem to employ.

Furthermore, there are indeed errors in your code. The most important one:

for(int i=0;i<=(Log(n)-1);i++)

Here, you’re iterating until floor(log(n)) - 1. But the pseudo-code states that you in fact need to iterate until ceil(log(n)) - 1.

Furthermore, consider this:

 for (int j=0;j<=n-1;j++)

This isn’t very usual. Normally, you’d write such code as follows:

for (int j = 0; j < n; ++j)

Notice that I used < instead of <= and adjusted the bounds from j - 1 to j. The same would hold for the outer loop, so you’d get:

for (int i = 0; i < std::log(n); ++i)

Finally, instead of std::powf, you can use the fact that pow(2, x) is the same as 1 << x (i.e. taking advantage of the fact that operations base 2 are efficiently implemented as bit operations). This means that you can simply write:

if (j >= 1 << i)
    x[j] += x[j - (1 << i)];

Answer 2

I think that std::partial_sum does what you want http://www.sgi.com/tech/stl/partial_sum.html

Answer 3

The quickest way to get your algorithm working: Drop the outer for(i...) loop, instead setting i to 0, and use only the inner for (j...) loop.

int main(){
    ...
    int i=0;
    for (int j=0;j<=n-1;j++){
        if (j>=(powf(2,i))){
            int t=powf(2,i);
            x[j]=x[j]+x[j-t];
        }
    }
    ...
}

Or equivalently:

for (int j=0; j<=n-1; j++) {
    if (j>=1)
        x[j] = x[j] + x[j-1];
}

...which is the intuitive way to do a prefix sum, and also probably the fastest non-parallel algorithm.

Wikipedia's algorithm is designed to be run in parallel, such that all of the additions are completely independent of each other. It reads all the values in, adds to them, then writes them back into the array, all in parallel. In your version, when you execute x[j]=x[j]+x[j-t], you're using the x[j-t] that you just added to, t iterations ago.

If you really want to reproduce Wikipedia's algorithm, here's one way, but be warned it will be much slower than the intuitive way above, unless you are using a parallelizing compiler and a computer with a whole bunch of processors.

int main() {
    int x[16]={39,21,20,50,13,18,2,33,49,39,47,15,30,47,24,1};
    int y[16];
    int n=sizeof(x)/sizeof(int);

    for (int i=0;i<=(Log(n)-1);i++){
        for (int j=0;j<=n-1;j++){
            y[j] = x[j];
            if (j>=(powf(2,i))){
                int t=powf(2,i);
                y[j] += x[j-t];
            }
        }
        for (int j=0;j<=n-1;j++){
            x[j] = y[j];
        }
    }

    for (int i=0;i<n;i++)
        cout<<x[i]<< "  ";
    cout<<endl;
}

Side notes: You can use 1<<i instead of powf(2,i), for speed. And as ergosys mentioned, your Log() function needs work; the values it returns are too high, which won't affect the partial sum's result in this case, but will make it take longer.