2D Fenwick Tree, in depth: rectangle sums with point and range updates in O(log n · log m)

The problem: rectangle sums that keep changing

Start with the static case. A 2D prefix-sum table P, where P[x][y] is the sum of a[1..x][1..y], answers any rectangle in O(1) by inclusion-exclusion. The trouble starts with updates: changing one cell a[x][y] changes every P entry to its lower right, up to n · m entries. If updates and queries interleave, you need a structure that trades a little query time for much cheaper writes.

The Fenwick version uses one integer per cell and no recursion. A segment tree wins when you need operations a Fenwick tree cannot express, such as range minimum with updates; see the segment tree article.

From one dimension to two

In one dimension, cell t[i] stores the sum of a over the half-open range (i - lowbit(i), i], where lowbit(i) = i & -i is the value of the lowest set bit. Updating index x walks upward with x += lowbit(x); a prefix query walks downward with x -= lowbit(x). Each walk visits at most about log2(n) + 1 cells.

The 2D tree applies the same decomposition independently to each axis. Cell t[i][j] stores the sum of a over the rectangle of rows (i - lowbit(i), i] and columns (j - lowbit(j), j]: a tree of trees, flattened into one table.

So every operation becomes two nested loops, an outer walk over i and, for each i, an inner walk over j, costing O(log n · log m).

What each cell stores

Read a few cells to make the rule concrete. Index 2 has lowbit 2, so t[2][2] covers rows 1-2 and columns 1-2: 3 + 0 + 5 + 6 = 14. Index 3 has lowbit 1, so t[3][4] covers row 3 alone across columns 1-4: 1 + 2 + 0 + 1 = 4. Index 4 has lowbit 4, so t[4][4] covers everything and equals the matrix total, 38.

Implementation

The implementation is short. It is 1-based because lowbit(0) is 0 and a walk starting at 0 would never move; index 0 is an unused sentinel row and column.

class Fenwick2D:
    def __init__(self, n, m):
        self.n, self.m = n, m
        self.t = [[0] * (m + 1) for _ in range(n + 1)]   # row 0 and col 0 unused

    def add(self, x, y, delta):          # a[x][y] += delta
        i = x
        while i <= self.n:
            j = y
            while j <= self.m:
                self.t[i][j] += delta
                j += j & -j
            i += i & -i

    def prefix(self, x, y):              # sum of a[1..x][1..y]
        s, i = 0, x
        while i > 0:
            j = y
            while j > 0:
                s += self.t[i][j]
                j -= j & -j
            i -= i & -i
        return s

    def rect(self, x1, y1, x2, y2):      # sum of a[x1..x2][y1..y2], inclusive
        return (self.prefix(x2, y2) - self.prefix(x1 - 1, y2)
                - self.prefix(x2, y1 - 1) + self.prefix(x1 - 1, y1 - 1))

Two details matter. The inner index j must be reset to y for every outer step; hoisting it out of the loop is the most common bug. And the structure stores sums, not values: to set a cell, read its current value with rect(x, y, x, y), then add the difference. The rectangle query is inclusion-exclusion: the full prefix, minus the bands above and to the left, plus the corner block subtracted twice.

Worked example

Query the sum of rows 2-3 and columns 2-4 of the example matrix: the cells 6, 3, 2 and 2, 0, 1, which total 14. The tree computes it as four prefixes.

The answer is 28 - 8 - 9 + 3 = 14. Check the first row of the table: t[3][4] is 4 (row 3) and t[2][4] is 24 (rows 1-2), so prefix(3, 4) = 28. The outer walk went 3 to 2 to 0; the inner walk from 4 stopped after one step each time.

Now apply a[3][3] += 5. The outer walk visits rows 3 and 4 (3 + 1 = 4, then 4 + 4 = 8, past n). For each, the inner walk visits columns 3 and 4. Four cells change: t[3][3], t[3][4], t[4][3] and t[4][4]. Re-running the query gives 19, exactly 5 more, because the updated cell lies inside the rectangle.

Building in linear time

Building by calling add for every cell costs O(n · m · log n · log m). A linear build exists. In 1D, copy a into t and then, for each i in increasing order, push t[i] into its parent t[i + lowbit(i)]. In 2D, run that pass along every row, then again along every column: the first pass makes each row a 1D tree over columns, the second combines rows.

def build(a):                       # a is n x m, 0-based input
    n, m = len(a), len(a[0])
    f = Fenwick2D(n, m)
    t = f.t
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            t[i][j] = a[i - 1][j - 1]
    for i in range(1, n + 1):       # pass 1: 1D build along each row
        for j in range(1, m + 1):
            p = j + (j & -j)
            if p <= m:
                t[i][p] += t[i][j]
    for i in range(1, n + 1):       # pass 2: 1D build along each column
        p = i + (i & -i)
        if p <= n:
            for j in range(1, m + 1):
                t[p][j] += t[i][j]
    return f

A single pass pushing each cell to its row, column and diagonal parents double-counts unless the order and subtraction are exactly right. Two separable passes are simpler; the code above was checked against cell-by-cell insertion on random matrices.

Range updates

Some workloads add a constant to a whole rectangle: paint a region of a grid, apply a discount to a price band in a date range. Two variants cover them.

Range add, point query. Store the 2D difference array instead of the values. Adding v to rows x1..x2 and columns y1..y2 becomes four point updates: +v at (x1, y1), -v at (x1, y2 + 1), -v at (x2 + 1, y1) and +v at (x2 + 1, y2 + 1), skipping any that fall outside the table. The value of a cell is then prefix(x, y) over the difference tree. One tree, same costs.

Range add, range sum. Summing the difference array twice gives the prefix sum of the values, and expanding that double sum yields four terms. Keep four Fenwick trees over the difference array d, storing d, d · i, d · j and d · i · j. Then:

# S(x, y) = sum of values a[1..x][1..y] after any number of range adds
S(x, y) = (x+1)(y+1) * A(x,y) - (y+1) * B(x,y) - (x+1) * C(x,y) + D(x,y)
#   A = prefix over d,  B = prefix over d*i,  C = prefix over d*j,  D = prefix over d*i*j

def point(x, y, v):          # one corner of a range add, written to all four trees
    A.add(x, y, v); B.add(x, y, v * x); C.add(x, y, v * y); D.add(x, y, v * x * y)

A range add writes the four corners into all four trees, sixteen point updates, and a rectangle sum is inclusion-exclusion over S, sixteen prefix walks. Constant factors are large but the bounds stay O(log n · log m). The formula was verified against a brute-force grid; a common mistake is x instead of x + 1 in the weights.

Memory and sparse coordinates

Memory is the real constraint. A dense table holds (n + 1)(m + 1) integers: a 4096 by 4096 grid of 64-bit sums is 128 MiB, and a 100,000 by 100,000 coordinate space is impossible. Three techniques handle large, sparse coordinate spaces.

• Coordinate compression. If all points are known up front, sort the distinct x values and distinct y values and map them to ranks. A query rectangle maps to rank ranges by binary search. This shrinks n and m to the number of distinct coordinates, which helps when either axis is small.
• Offline per-node compression. When the point set is known but the grid of distinct values is still too large, give each outer node i its own sorted list of the y values that will ever reach it, and an inner 1D tree sized to that list. Total memory is O(N log N) for N points, and each inner access becomes a binary search plus a 1D walk.
• Hash-map cells. For online sparse updates, key t by (i, j) in a hash map; memory is O(U · log n · log m) for U updates, at a large constant-factor cost.

A fourth option is often better than all of them: if queries can be answered offline, sweep one axis in sorted order and keep a 1D Fenwick tree over the other. Counting points in rectangles then needs O(N) memory and O((N + Q) log N) time.

Performance and operational notes

In compiled languages, lay the table out as one contiguous array indexed i · (m + 1) + j, not as an array of row arrays. The inner loop walks along a row, so keep j as the fastest-moving index. Use 64-bit cells: the top-level cells hold sums of millions of values. Minimum and other non-invertible operations do not support rectangle queries, because inclusion-exclusion needs subtraction.

In services, one add touches up to log n · log m cells, so a concurrent reader can see a half-applied update. Serialise access with a lock, or batch updates and swap in a rebuilt tree.

Failure modes

• Off-by-one at zero. Calling add(0, y) loops forever because lowbit(0) is 0. Validate 1 ≤ x ≤ n at the API boundary.
• Inner index not reset. Sums are correct for small rectangles and wrong for larger ones, which makes the bug intermittent in tests.
• Wrong inclusion-exclusion corner. Using x1 instead of x1 - 1 drops the first row of the rectangle. Test single-cell rectangles against a brute-force matrix.
• Overflow in the four-tree variant. The d · i · j tree grows with the coordinates and overflows long before the plain tree.
• Memory explosion. Allocating a dense table for a sparse coordinate space fails at startup or, worse, under load. Size the table from the data, not from the coordinate range.

Trade-offs

Choose the 2D Fenwick tree when the operation is invertible (sum, count, XOR), the grid is dense and of moderate size, and updates interleave with queries. Choose a static prefix-sum table when the data never changes. Choose a 2D segment tree or a quadtree when you need minimum, maximum or lazy assignment. Choose an offline sweep with a 1D tree when all queries are known in advance. And for spatial search over points rather than aggregation over cells, use a k-d tree or, for overlapping intervals on one axis, an interval tree.

The structure generalises to d dimensions with d nested loops and O(log^d n) per operation, but memory grows as the product of all axes, so beyond three dimensions it is rarely the right tool.

What to do next

1. Implement Fenwick2D from memory, then write a test that compares rect against a brute-force matrix on random updates and random rectangles, including single cells and full-grid queries.
2. Reproduce the worked example by hand: trace prefix(3, 4) and the update at (3, 3) cell by cell.
3. Implement the two-pass build and assert it produces the same table as repeated add.
4. Add the range-add, point-query variant with the four-corner difference update.
5. Implement the four-tree range-add, range-sum variant and verify it against brute force; check the (x + 1) weights explicitly.
6. Estimate memory for your real coordinate space before choosing dense, compressed or offline.
7. If your queries are known in advance, try the offline sweep with a 1D tree and compare speed.