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 O(1) time and retrieval takes O(N), in the case (b) update takes O(N) and retrieval O(1) time. [math]\displaystyle{ N }[/math] is the length of the sequence. The Fenwick tree data structure (although it's not really a tree) has both update and retrieval taking [math]\displaystyle{ O(\log N) }[/math] time.

The idea is to have two actions on positive integers, say [math]\displaystyle{ A }[/math] and [math]\displaystyle{ B }[/math], so that orbit of [math]\displaystyle{ x }[/math] under action of [math]\displaystyle{ A }[/math], and orbit of [math]\displaystyle{ y }[/math] under action of [math]\displaystyle{ B }[/math] intersect in exactly one point if [math]\displaystyle{ x \leq y }[/math], and none otherwise. Then maintain an array of size [math]\displaystyle{ N }[/math] initialized to 0. On each update of a sequence element [math]\displaystyle{ i }[/math], update all the elements of the array with indices in the [math]\displaystyle{ A }[/math] orbit of [math]\displaystyle{ i }[/math] by that same amount. To retrieve the partial sum of elements up to [math]\displaystyle{ j }[/math], sum up all the elements from the array with indices in the [math]\displaystyle{ B }[/math] orbit of [math]\displaystyle{ j }[/math]. In the case (a) above [math]\displaystyle{ A }[/math] is identity and [math]\displaystyle{ B(j) = j-1 }[/math], for (b) it's [math]\displaystyle{ A(i) = i+1 }[/math] and [math]\displaystyle{ B }[/math] is the identity.

Here's the punchline: for a positive integer [math]\displaystyle{ n }[/math], if we call [math]\displaystyle{ 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]\displaystyle{ n }[/math] (so e.g. [math]\displaystyle{ r(6) = 2 }[/math] and [math]\displaystyle{ r(8) = 8 }[/math]), then [math]\displaystyle{ A(i) = i + r(i) }[/math] and [math]\displaystyle{ B(i) = i - r(i) }[/math] are two such actions. What's more, in C and its ilk, [math]\displaystyle{ r(i) }[/math] is coded as -i&i.

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

Here's a simple implementation in C++, add error checking as needed

#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; i<int(v_.size()); 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;
    }
};