2D Binary Indexed Tree, in depth: sparse coordinates, the offline-planned BIT of sorted lists, sweeps, and dominance chains

Recap: the walks a 2D BIT makes

A 1D Fenwick tree stores, at index i, the sum of the last lowbit(i) elements ending at i, where lowbit(i) = i & -i. A prefix query walks down by clearing the lowest set bit; an update walks up by adding it. Each walk visits at most about log2(n) nodes. If the derivation is new, read the Fenwick tree in depth first.

The 2D version nests two walks: the outer walk over x chooses which rows of the table to touch, and inside each one an inner walk over y. Every operation costs about log n times log m steps. The planned variant keeps exactly this shape. It changes only what an inner index means: instead of a raw y coordinate, the inner index is the rank of y within the sorted list of y values that the outer node will ever see.

The planning pass

The structure needs every update point before the first update. That is the offline part. Queries can be arbitrary and online; only the update positions must be known. Building takes three steps.

1. Compress x. Sort the distinct x values; an x coordinate becomes its 1-based rank.
2. Simulate every update's outer walk. For each point (x, y), walk i from rank(x) upward with i += i & -i and append y to node i's list. This is exactly the set of nodes a real add will touch, so the list is complete by construction.
3. Sort and deduplicate each node's list, and allocate an inner array one longer than it.

An add at (x, y) then walks the same outer nodes and, in each one, finds y by binary search. A query for prefix(x, y), the sum over points with px at most x and py at most y, walks the outer nodes downward from the number of distinct x values at most x, and in each node starts the inner walk at the number of listed y values at most y. Query coordinates need not be in any list; bisect handles them.

import bisect

class OfflineBIT2D:
    """Point add / dominance sum over arbitrary coordinates, after a planning pass."""

    def __init__(self, points):
        # points: every (x, y) that will ever be updated, known before the first add
        self.xs = sorted({x for x, _ in points})
        n = len(self.xs)
        self.ys = [[] for _ in range(n + 1)]
        for x, y in points:
            i = bisect.bisect_left(self.xs, x) + 1
            while i <= n:                        # same walk an add will take
                self.ys[i].append(y)
                i += i & -i
        for i in range(1, n + 1):
            self.ys[i] = sorted(set(self.ys[i]))
        self.t = [[0] * (len(col) + 1) for col in self.ys]

    def add(self, x, y, v):
        i = bisect.bisect_left(self.xs, x) + 1
        while i < len(self.t):
            col, row = self.ys[i], self.t[i]
            j = bisect.bisect_left(col, y) + 1  # y is guaranteed present in col
            while j < len(row):
                row[j] += v
                j += j & -j
            i += i & -i

    def prefix(self, x, y):
        """Sum of values at points with px <= x and py <= y."""
        i = bisect.bisect_right(self.xs, x)
        s = 0
        while i > 0:
            col, row = self.ys[i], self.t[i]
            j = bisect.bisect_right(col, y)
            while j > 0:
                s += row[j]
                j -= j & -j
            i -= i & -i
        return s

    def rect(self, x1, y1, x2, y2):
        # integer coordinates, so x1 - 1 is the predecessor of x1
        return (self.prefix(x2, y2) - self.prefix(x1 - 1, y2)
                - self.prefix(x2, y1 - 1) + self.prefix(x1 - 1, y1 - 1))

Note the two bisect flavours. Ranks in an add use bisect_left plus one, because the point is present and we want its own position. Counts in a query use bisect_right, because we want how many listed values are at most the query value. Mixing them up is the most common bug in this structure, and it fails only on queries that land exactly on a stored coordinate.

Worked example: six points, one of them at a billion

Take six points with one outlier coordinate: (3, 7), (10, 2), (10, 9), (25, 4), (40, 7) and (10^9, 1). The distinct x values are 3, 10, 25, 40 and 10^9, ranks 1 to 5. The planning pass sends each point up its outer walk. The point at x = 3 (rank 1) visits nodes 1, 2 and 4; the two points at x = 10 (rank 2) visit 2 and 4; x = 25 (rank 3) visits 3 and 4; x = 40 (rank 4) visits only 4; the outlier (rank 5) visits only 5. After deduplication the node lists are:

That is 10 inner cells. A dense table over the compressed grid would need 5 by 5 = 25, and the raw grid is out of the question. After adding weight 1 at every point, prefix(10, 7) walks outer node 2 (covering x ranks 1 and 2), counts the listed y values at most 7 in [2, 7, 9], which is two, and sums inner cells to get 2: the points (3, 7) and (10, 2). Then it clears the low bit and stops. rect(5, 1, 40, 8) returns 3, for (10, 2), (25, 4) and (40, 7), and the rectangle over everything returns 6. Each value matched a brute-force count.

Memory and time

Each update point is listed in at most about log2(n) outer nodes, so total inner cells are O(N log N). Measured on random coordinates up to 10^9, the planning pass recorded 7.17 list entries per point at N = 10^4 and 8.78 at N = 10^5, before deduplication; the average is about half of log2(N) because a walk from a random rank is shorter than the worst case. At 10^5 points that is under a million 64-bit cells, a few megabytes, against 80 GB for a dense compressed table.

Time per operation is O(log n log N): the outer walk, and in each node a binary search plus an inner walk. The binary searches add a constant factor of two to three over a dense 2D BIT in practice, and the pointer-chasing between per-node arrays hurts cache behaviour. In a compiled language, pack all inner arrays into one flat buffer with an offset per node, and store the y lists in a second flat buffer with the same offsets.

When points are static: the sweep

If the points are static and all queries are known in advance, you do not need a 2D structure at all. Sweep x from left to right, insert points into a 1D Fenwick tree over compressed y as the sweep passes them, and split each rectangle query into two events: add the y-range count at x2, subtract it at x1 - 1. Memory is O(N) and time is O((N + Q) log N).

def offline_rect_counts(points, queries):
    """Static points, rectangle count queries: sweep x, 1D BIT over y."""
    ys = sorted({y for _, y in points})
    m = len(ys)
    bit = [0] * (m + 1)

    def add(j):
        while j <= m:
            bit[j] += 1
            j += j & -j

    def pre(y):
        j, s = bisect.bisect_right(ys, y), 0
        while j > 0:
            s += bit[j]
            j -= j & -j
        return s

    events = []                                   # (x, sign, y_lo, y_hi, qid)
    for qid, (x1, y1, x2, y2) in enumerate(queries):
        events.append((x2, +1, y1, y2, qid))
        events.append((x1 - 1, -1, y1, y2, qid))
    events.sort()
    pts = sorted(points)
    ans, p = [0] * len(queries), 0
    for x, sign, ylo, yhi, qid in events:
        while p < len(pts) and pts[p][0] <= x:
            add(bisect.bisect_left(ys, pts[p][1]) + 1)
            p += 1
        ans[qid] += sign * (pre(yhi) - pre(ylo - 1))
    return ans

On the six example points with queries (5, 1, 40, 8), the full range and (11, 3, 30, 8), it returns 3, 6 and 1. Use the sweep whenever it applies; it is simpler, faster and smaller. Reach for the planned 2D BIT when updates and queries interleave and a query must see exactly the updates before it.

Prefix maxima and dominance chains

Fenwick trees also answer prefix maxima, provided values at a position only ever increase. Rectangle maxima are impossible, because inclusion-exclusion needs subtraction, but dominance maxima are fine, and they are what dynamic programs need. A classic case is the longest chain of items strictly increasing in three keys (a, b, c): sort by a, and for each item ask for the best chain ending at an item with smaller b and smaller c.

import itertools

class OfflineBIT2DMax(OfflineBIT2D):
    """Same planning pass; cells hold maxima, and values may only increase."""

    def raise_to(self, x, y, v):
        i = bisect.bisect_left(self.xs, x) + 1
        while i < len(self.t):
            col, row = self.ys[i], self.t[i]
            j = bisect.bisect_left(col, y) + 1
            while j < len(row):
                row[j] = max(row[j], v)
                j += j & -j
            i += i & -i

    def best_below(self, x, y):
        """Max over points with px < x and py < y (strict on both keys)."""
        i, s = bisect.bisect_left(self.xs, x), 0
        while i > 0:
            col, row = self.ys[i], self.t[i]
            j = bisect.bisect_left(col, y)
            while j > 0:
                s = max(s, row[j])
                j -= j & -j
            i -= i & -i
        return s

def longest_chain(items):
    """Longest sequence strictly increasing in all of (a, b, c)."""
    tree = OfflineBIT2DMax([(b, c) for _, b, c in items])   # same planning, max instead of +
    best = 0
    for _, group in itertools.groupby(sorted(items), key=lambda t: t[0]):
        group = list(group)
        # query the whole group before updating, so equal a never chains to itself
        lens = [tree.best_below(b, c) + 1 for _, b, c in group]
        for (_, b, c), length in zip(group, lens):
            tree.raise_to(b, c, length)
        best = max([best] + lens)
    return best

OfflineBIT2DMax reuses the planning pass and changes two things: inner cells keep the maximum instead of the sum, and the query uses bisect_left at both levels so that it covers strictly smaller b and c. For the items (1, 5, 2), (2, 6, 3), (2, 7, 9), (3, 1, 1), (4, 8, 4) and (5, 9, 10) it returns 4, the chain through (1, 5, 2), (2, 6, 3), (4, 8, 4) and (5, 9, 10). The grouping step matters: without it, two items with equal a could chain.

Choosing a structure

Failure modes and testing

• An update at an unplanned point. The binary search lands on a neighbour and corrupts a different cell, silently. In production code, check that the found value equals y and raise otherwise.
• bisect_left versus bisect_right. Wrong flavour in the query undercounts points sitting exactly on the boundary. Test boundaries explicitly.
• Predecessor on non-integer keys. x1 - 1 assumes integers. For floats or timestamps, use half-open queries and count values strictly below x1 with bisect_left.
• Overflow. Top-level cells sum everything; use 64-bit integers in compiled languages.
• Max variant with decreasing values. A Fenwick max cannot lower a value. If a position's value can drop, switch to a segment tree.
• Testing. Keep a brute-force O(N) oracle and run randomized instances with duplicate coordinates, negative weights and ties; the code in this article passed hundreds of them.

What to do next

1. Implement OfflineBIT2D, then reproduce the worked example and its 10-cell layout.
2. Write the brute-force oracle and a randomized test with duplicates and boundary queries.
3. Solve one problem both ways, sweep and planned BIT, and time them at N = 10^5.
4. Derive OfflineBIT2DMax and test longest_chain against an O(N^2) dynamic program.
5. Flatten the per-node arrays into two buffers with offsets and measure the speed-up.
6. Read the persistent segment tree for the online alternative when update positions cannot be planned.