Range Tree, in depth: multi-dimensional orthogonal range queries, canonical nodes, fractional cascading and rank-space ties

The idea: a tree of trees

A 2-D range tree is a tree of trees. The primary structure is a balanced binary tree over the points sorted by x. Each node represents a contiguous run of that order: the root all n points, its children the left and right halves, and so on down to single points. Every node also stores an associated structure: the same points, sorted by y. For two dimensions the associated structure can be a sorted array; for d dimensions it is a (d-1)-dimensional range tree, which is where the recursion in the name comes from.

A query interval on x can be covered exactly by O(log n) nodes, the canonical nodes, the same decomposition a segment tree uses. Inside each canonical node every point already satisfies the x condition, so only y remains, and that is a one-dimensional query: two binary searches on the node's y-sorted array. The answer is the union over canonical nodes, and the nodes are disjoint, so counts add up without double counting.

Building and querying a 2-D range tree

The implementation below builds the tree bottom-up exactly like merge sort: each node's y list is the merge of its children's lists, so construction takes O(n log n) time. Nodes are numbered heap-style, root 1 and children 2i and 2i+1, and each covers a half-open range of x-ranks. The query first converts the x-interval to a rank interval with two binary searches, then collects canonical nodes, then binary-searches y in each.

from bisect import bisect_left, bisect_right

class RangeTree2D:
    """Static 2-D range tree over x-ranks. Ties in x are safe."""
    def __init__(self, points):
        self.pts = sorted(points)                  # by x, then y
        self.xs = [p[0] for p in self.pts]
        self.n = len(self.pts)
        self.ys = {}                               # node -> sorted list of (y, x)
        if self.n:
            self._build(1, 0, self.n)

    def _build(self, node, lo, hi):
        if hi - lo == 1:
            self.ys[node] = [(self.pts[lo][1], self.pts[lo][0])]
            return
        mid = (lo + hi) // 2
        self._build(2 * node, lo, mid)
        self._build(2 * node + 1, mid, hi)
        a, b = self.ys[2 * node], self.ys[2 * node + 1]
        out, i, j = [], 0, 0                       # merge, as in merge sort
        while i < len(a) and j < len(b):
            if a[i] <= b[j]:
                out.append(a[i]); i += 1
            else:
                out.append(b[j]); j += 1
        self.ys[node] = out + a[i:] + b[j:]

    def _canonical(self, node, lo, hi, qlo, qhi, acc):
        if qhi <= lo or hi <= qlo:
            return
        if qlo <= lo and hi <= qhi:
            acc.append(node)
            return
        mid = (lo + hi) // 2
        self._canonical(2 * node, lo, mid, qlo, qhi, acc)
        self._canonical(2 * node + 1, mid, hi, qlo, qhi, acc)

    def query(self, x1, x2, y1, y2):
        qlo = bisect_left(self.xs, x1)             # x-interval -> rank interval
        qhi = bisect_right(self.xs, x2)
        nodes, found = [], []
        if self.n and qlo < qhi:
            self._canonical(1, 0, self.n, qlo, qhi, nodes)
        for node in nodes:
            ys = self.ys[node]
            a = bisect_left(ys, (y1, float("-inf")))
            b = bisect_right(ys, (y2, float("inf")))
            found.extend((x, y) for y, x in ys[a:b])
        return found

For counting rather than reporting, replace the extend with count += b - a; the query then costs O(log^2 n) regardless of how many points match. This code was checked against brute force on thousands of random queries, including point sets with repeated x values and repeated points.

Worked example on eight points

Take eight points sorted by x: (2,7), (3,2), (5,5), (6,9), (8,1), (9,6), (11,4), (12,8), with x-ranks 0 to 7. The query is x in [3, 11] and y in [3, 8].

1. Rank conversion: bisect_left(xs, 3) is 1 and bisect_right(xs, 11) is 7, so the query covers ranks 1 to 6.
2. Decomposition: the root (ranks 0-7) is partial, so recurse. Node 2 (0-3) is partial; its child node 4 (0-1) is partial and yields leaf node 9 (rank 1); node 5 (2-3) is fully inside. Node 3 (4-7) is partial; node 6 (4-5) is fully inside; node 7 (6-7) is partial and yields leaf node 14 (rank 6). Canonical nodes: 9, 5, 6, 14.
3. Node 9 has y list [2]; nothing in [3, 8].
4. Node 5 has [5, 9]; the binary searches select 5, the point (5,5).
5. Node 6 has [1, 6]; they select 6, the point (9,6).
6. Node 14 has [4]; it selects (11,4).

The answer is (5,5), (9,6) and (11,4), found by visiting 11 nodes, of which 4 are canonical and each get 2 binary searches on short lists. A scan of the x-range would have examined six points; on a million points with a wide x-range and a narrow y-range the difference is the whole point of the structure.

Why the bounds hold

Space: every point appears in the y list of each ancestor of its leaf, one per level, and there are about log2 n + 1 levels, so the 2-D tree stores O(n log n) entries. The eight-point example stores 32 entries, eight points times four levels. Query: the decomposition visits O(log n) nodes and selects at most two per level as canonical, and each canonical node costs one O(log n) binary search, so reporting is O(log^2 n + k).

In d dimensions the associated structure of each node is a (d-1)-dimensional range tree. Each extra dimension multiplies space by a log factor and query time by another, giving O(n log^(d-1) n) space, O(n log^(d-1) n) preprocessing with care, and O(log^d n + k) query time. Fractional cascading, below, removes one log factor from the query, giving O(log^(d-1) n + k), so O(log n + k) in two dimensions.

Fractional cascading

Every canonical node repeats a binary search for the same y1 and y2 on a list that is a subset of its parent's list. Fractional cascading exploits that. If you know the position of y1 in the parent's list, you can find its position in each child's list in O(1), provided the parent stores, for each position, how many of the items before it came from the left child. Search the root once, then carry the two positions down.

class CascadedRangeTree2D(RangeTree2D):
    """Counting with fractional cascading: two binary searches total."""
    def __init__(self, points):
        super().__init__(points)
        self.left_cnt = {}
        for node, ys in self.ys.items():
            if 2 * node in self.ys:                # internal node
                left, cnt, i = self.ys[2 * node], [0], 0
                for item in ys:                    # did this merged item come from left?
                    if i < len(left) and left[i] == item:
                        i += 1
                    cnt.append(i)
                self.left_cnt[node] = cnt          # cnt[a] = left items among ys[:a]

    def count(self, x1, x2, y1, y2):
        qlo, qhi = bisect_left(self.xs, x1), bisect_right(self.xs, x2)
        if not self.n or qlo >= qhi:
            return 0
        root = self.ys[1]
        a = bisect_left(root, (y1, float("-inf")))
        b = bisect_right(root, (y2, float("inf")))
        return self._walk(1, 0, self.n, qlo, qhi, a, b)

    def _walk(self, node, lo, hi, qlo, qhi, a, b):
        if qhi <= lo or hi <= qlo or a == b:
            return 0
        if qlo <= lo and hi <= qhi:
            return b - a
        cnt, mid = self.left_cnt[node], (lo + hi) // 2
        la, lb = cnt[a], cnt[b]                    # follow pointers, no search
        return (self._walk(2 * node, lo, mid, qlo, qhi, la, lb) +
                self._walk(2 * node + 1, mid, hi, qlo, qhi, a - la, b - lb))

In the example the root's list is [1, 2, 4, 5, 6, 7, 8, 9] and its prefix counts are [0, 0, 1, 1, 2, 2, 3, 3, 4]: the 2, 5, 7 and 9 came from the left child. For y in [3, 8] the root positions are a = 2 and b = 7; the left child gets 1 and 3, the right child 1 and 4, and so on down, with no further searching. This version counts in O(log n). The full layered range tree applies the same idea to reporting.

Ties and rank space

Textbook range trees assume no two points share an x coordinate, then split nodes on coordinate values. Real data has ties everywhere: timestamps at second resolution, prices in cents, grid coordinates. With value-based splitting, equal x values can land on both sides of a split and the query either misses or double counts them. The implementation above sidesteps this by working in rank space: points are sorted by (x, y), nodes cover ranks, and the query turns the x-interval into ranks with bisect_left and bisect_right, which include every tied point exactly once. The same trick, applied to y by storing (y, x) pairs and searching with infinite sentinels, handles ties in the second dimension. Composite keys of coordinate and id are the general fix when points are identical.

Updates, memory and engineering

The structure is static. Inserting a point touches one y list per level, and a sorted array insert is O(n), so dynamic use needs a different associated structure, such as a balanced tree or a 2-D Fenwick tree over compressed coordinates for counting. A common compromise is the logarithmic method: keep a few static trees of doubling sizes, rebuild by merging, and query all of them. Deletions can be tombstones with periodic rebuilds.

Memory dominates at scale. A million points give about 21 levels, so about 21 million entries; at 4 bytes for y and 4 for a point index that is roughly 168 MB, and the cascade counts add about 84 MB more. Production layouts drop the per-node dictionaries used here for one flat array per level, which is the merge-sort tree, and store indices rather than tuples. Wavelet trees answer the same 2-D counting queries in O(n log n) bits instead of words. If your queries are offline, a sweep over x with a Fenwick tree on y answers all of them in O((n + q) log n) time with linear memory.

Failure modes

• Value-based splits with duplicate coordinates silently lose or double count points; use rank space.
• Open versus closed intervals mixed between dimensions produce off-by-one errors at boundaries. Decide once and test boundary points explicitly.
• Memory blow-up in three or more dimensions. A 3-D tree on a million points holds roughly n times 21 squared, over 400 million entries; check before building.
• Using it for dynamic data with array-based associated structures turns inserts into full rebuilds.
• Recursion depth is fine for balanced trees but iterative traversal is faster in Python and avoids overhead on large builds.

Trade-offs and alternatives

The range tree wins on worst-case query time; the k-d tree and R-tree win on memory and updates. In practice many systems use a k-d tree or R-tree and reach for a range tree, or its flat merge-sort-tree form, only when profiling shows worst-case queries dominate.

What to do next

1. Write down your dimensions, query mix, update rate and memory budget before choosing.
2. Implement the rank-space version above and test it against brute force with duplicates.
3. Measure memory as levels times n times bytes per entry before building at scale.
4. Switch to flat per-level arrays when the per-node version works.
5. Add fractional cascading if profiling shows the per-node binary searches dominate.
6. For dynamic data, use the logarithmic method or a Fenwick tree over compressed coordinates.
7. Compare against a k-d tree on your real queries, not on the asymptotics.