Updating partial sums with Fenwick tree

From Polymath Wiki
Jump to: navigation, search

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 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 \lt x \lt 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 (-i)&i. (Actually, -i&i 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

#include <vector>

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

    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;