View source for Updating partial sums with Fenwick tree

The two straightforward ways to code updating and retrieval of partial sums of sequences are: (a) maintain the array of sequence elements - on each update of a sequence element update just the sequence element, and each time you need a partial sum you compute it from scratch, (b) maintain the array of partial sums - on each update of a sequence element update all the partial sums this element contributes to. In the case (a) update takes <math>O(1)</math> time and retrieval takes <math>O(N)</math>, in the case (b) update takes <math>O(N)</math> and retrieval <math>O(1)</math> time, where <math>N</math> is the length of the sequence. The [http://www.cs.ubc.ca/local/reading/proceedings/spe91-95/spe/vol24/issue3/spe884.pdf Fenwick tree] data structure has both update and retrieval taking <math>O(\log N)</math> time in the worst case, while still using the same <math>O(N)</math> amount of space.

The idea is to have two actions on positive integers, say <math>A</math> and <math>B</math>, which satisfy the following:

'''Property 1''' For any two positive integers <math>x</math> and <math>y</math>, orbit of <math>x</math> under action of <math>A</math>, and orbit of <math>y</math> under action of <math>B</math> intersect in exactly one point if <math>x \leq y</math>, and nowhere otherwise. 

The general algorithm when we have <math>A</math> and <math>B</math> that satisfy the property above is as follows. We will maintain an array of size <math>N</math> with all elements initialized to 0. By updating something we mean adding a quantity to it, whether positive or negative. On each update of the <math>i^{th}</math> sequence element, update all the elements of the array with indices in the <math>A</math> orbit of <math>i</math> by that same amount. To retrieve the partial sum of elements <math>1</math> thru <math>j</math>, sum up all the elements from the array with indices in the <math>B</math> orbit of <math>j</math>. 

The two straightforward algorithms we discussed in the first paragraph already fit into this scheme. In the case (a) above, <math>A</math> is identity and <math>B(j) = j-1</math>, for (b) it's <math>A(i) = i+1</math> and <math>B</math> is the identity. It's obvious that in both cases the pair <math>A</math> and <math>B</math> satisfy Property 1.

Here's the punchline: for a positive integer <math>n</math>, if we call <math>r(n)</math> the integer we get when we isolate the rightmost, i.e. the least significant, set bit in the binary representation of <math>n</math> (so e.g. <math>r(6) = 2</math> and <math>r(8) = 8</math>), then <math>A(i) = i + r(i)</math> and <math>B(i) = i - r(i)</math> satisfy Property 1.

Proof: <math>A</math> increases its argument while <math>B</math> decreases it, so if <math>x=y</math> then <math>A^0(x) = B^0(y)</math>, and that's the only intersection. If <math>0 < x < y</math>, then binary representations of <math>x</math> and <math>y</math> are <math>y = a1b_y</math> and <math>x = a0b_x</math> where <math>a</math>, <math>b_x</math>, <math>b_y</math> are (possibly empty) binary sequences, and <math>b_x</math> and <math>b_y</math> are of same length <math>|b|</math>. We may have had to pad the representation of <math>x</math> by zeros to the left as needed for this. The action <math>B</math> unsets the rightmost set bit, so eventually <math>y</math> will under action of <math>B</math> arrive at <math>a10^{|b|}</math>, if it's not of that form to begin with. If <math>b_x</math> has no set bits, then one more action of <math>B</math> on <math>y</math> will make it equal to <math>x</math>. If <math>b_x</math> has set bits, then <math>A</math>, which turns a number of the form <math>z01^{m}0^{n}</math> into <math>z10^{m+n}</math> will eventually bring <math>x</math> to <math>a10^{|b|}</math>, which is where we left <math>y</math>. By the argument from the beginning, this is the only place the two orbits intersect. QED

What's more, in C and its ilk, <math>r(i)</math> is coded as <code>(-i)&i</code>. (Actually, <code>-i&i</code> will parse just fine since unary minus has higher precedence in C than bitwise and.)

Here's a simple implementation in C++, add error checking as needed
<pre>
#include <vector>

class fenwick_tree {
    std::vector<int> v_;

  public:
    fenwick_tree(int N) : v_(N+1) {}

    // neither add() nor get() will work 
    // for k outside of [1, N] range

    // add val to the k-th element
    void add(int k, int val) {
        for (int i=k, e=v_.size(); i<e; i += -i&i)
            v_[i] += val;
    }

    // get sum of elements 1 thru k
    int get(int k) const {
        int r=0;
        for (int i=k; i>0; i -= -i&i)
            r += v_[i];
        return r;
    }
};
</pre>