Square Root Decomposition, in depth: blocks, lazy tags, sorted blocks and rebuilding on a schedule

The idea and the cost model

Take an array of n elements and a block size B. Element i lives in block i // B, so there are about n / B blocks. For each block store a summary, for range sums simply the block's total. A query over [l, r] has three parts: the tail of l's block, the whole blocks strictly between, and the head of r's block. The partial parts touch at most 2B elements; the middle touches at most n / B summaries.

So a query costs O(B + n / B). The sum of those two terms is smallest when they are equal, at B = sqrt(n), giving O(sqrt n). For n = 1,000,000 that is about 1,000 blocks of 1,000 elements, and a query reads at most about 3,000 values instead of a million. A point update changes one element and one block summary: O(1). Building the structure is a single O(n) pass.

Worked example by hand

Use a = [5, 2, 7, 1, 3, 8, 6, 4, 9] with B = 3. The block sums are 14, 12 and 19.

1. Query sum of indices 1 to 7. Index 1 is in block 0 and index 7 in block 2. Left partial: indices 1 and 2, giving 2 + 7 = 9. Middle: block 1, sum 12. Right partial: indices 6 and 7, giving 6 + 4 = 10. Total 31, matching a direct sum.
2. Add 10 to every element in indices 2 to 6. Index 2 is a partial piece of block 0: update the element to 17 and the block sum to 24. Block 1 is covered entirely, so do not touch its elements: record a lazy tag of +10 and raise its sum by 10 times 3, to 42. Index 6 is a partial piece of block 2: update it to 16 and the sum to 29.
3. Query indices 1 to 7 again. Left partial: 2 + 17 = 19, plus block 0's lazy tag (zero) times two elements. Middle: block 1's sum, 42, which already includes the lazy add. Right partial: 16 + 4 = 20. Total 81. Check: the original 31 plus 10 added to five elements inside the range is 81.

Range add and range sum in code

The implementation keeps three arrays: the elements, the true sum of each block, and a lazy add per block that applies to every element in it but has not been written into them. Element reads in partial blocks therefore add the block's lazy value.

import math


class SqrtDecomp:
    """Range add and range sum over a list, blocks of size about sqrt(n)."""

    def __init__(self, a, block=None):
        self.n = len(a)
        self.B = block or max(1, math.isqrt(self.n))
        self.a = list(a)
        nb = (self.n + self.B - 1) // self.B
        self.sum = [0] * nb                    # true sum of each block
        self.lazy = [0] * nb                   # pending add for every element
        for i, v in enumerate(a):
            self.sum[i // self.B] += v

    def add(self, l, r, v):                    # inclusive l..r
        B = self.B
        bl, br = l // B, r // B
        if bl == br:
            for i in range(l, r + 1):
                self.a[i] += v
            self.sum[bl] += v * (r - l + 1)
            return
        for i in range(l, (bl + 1) * B):       # left partial block
            self.a[i] += v
        self.sum[bl] += v * ((bl + 1) * B - l)
        for b in range(bl + 1, br):            # whole blocks: O(1) each
            self.lazy[b] += v
            self.sum[b] += v * B
        for i in range(br * B, r + 1):         # right partial block
            self.a[i] += v
        self.sum[br] += v * (r - br * B + 1)

    def query(self, l, r):                     # inclusive l..r
        B = self.B
        bl, br = l // B, r // B
        if bl == br:
            return sum(self.a[l:r + 1]) + self.lazy[bl] * (r - l + 1)
        s = sum(self.a[l:(bl + 1) * B]) + self.lazy[bl] * ((bl + 1) * B - l)
        s += sum(self.sum[bl + 1:br])
        s += sum(self.a[br * B:r + 1]) + self.lazy[br] * (r - br * B + 1)
        return s

Partial loops stop at r, so a shorter final block needs no special case. This code was checked against a brute-force array on thousands of random operations with block sizes from 1 to 9; do the same with any variant you write, because off-by-one errors in the partial ranges are the classic bug.

Choosing the block size

B = sqrt(n) balances a query that costs B + n / B. When operations cost differently, balance the real costs instead. If partial-block work costs c₁ per element and whole-block work c₂ per block, the optimum is B = sqrt(n · c₂ / c₁). The sorted-block variant below is an example: whole blocks cost a binary search, so B grows to about sqrt(n log n).

Constants matter more than the formula suggests. Scanning a contiguous partial block is cache-friendly and vectorises well, while each block summary may be a separate cache line, so the measured optimum is often larger than sqrt(n). Benchmark a few sizes around the theoretical value on realistic data and pick the fastest; a factor of two either side is common.

Richer summaries: sorted blocks

The summary can be anything you can rebuild from a block in about O(B) and query quickly. A powerful choice is a sorted copy of the block, which answers questions such as how many values in a[l..r] are at most x, with point assignments allowed. A whole block answers with one binary search; partial blocks are scanned; an update removes the old value from its block's sorted list and inserts the new one.

import bisect, math


class CountLE:
    """Count values <= x in a[l..r], with point assignment. Blocks kept sorted."""

    def __init__(self, a):
        self.a = list(a)
        n = len(a)
        self.B = max(1, int(math.sqrt(n * max(1, math.log2(n or 1)))))
        self.blocks = [sorted(a[i:i + self.B]) for i in range(0, n, self.B)]

    def assign(self, i, v):                     # O(B)
        blk = self.blocks[i // self.B]
        blk.pop(bisect.bisect_left(blk, self.a[i]))
        bisect.insort(blk, v)
        self.a[i] = v

    def count_le(self, l, r, x):                # O(B + (n/B) log B)
        B, total = self.B, 0
        bl, br = l // B, r // B
        if bl == br:
            return sum(1 for v in self.a[l:r + 1] if v <= x)
        total += sum(1 for v in self.a[l:(bl + 1) * B] if v <= x)
        for b in range(bl + 1, br):
            total += bisect.bisect_right(self.blocks[b], x)
        total += sum(1 for v in self.a[br * B:r + 1] if v <= x)
        return total

A segment tree can answer the same query with a sorted list per node (a merge sort tree), but updates then touch O(log n) lists of growing size. The block version is simpler and handles updates gracefully. For the static case, a persistent segment tree answers it in O(log n).

Decomposing time and values

The same square-root balance applies to sequences of operations. Keep a structure that is fast to query but expensive to rebuild, such as a prefix-sum array, and buffer updates in a small pending list. Each query reads the prebuilt answer and corrects it by scanning the pending list; when the list reaches K entries, rebuild in O(n). With K = sqrt(n), queries cost O(sqrt n) and updates O(sqrt n) amortised.

class LazyPrefix:
    """Prefix sums with point adds: rebuild after K pending updates."""

    def __init__(self, a, K=None):
        self.a = list(a)
        self.K = K or max(1, math.isqrt(len(a)))
        self.pending = []                       # (index, delta) not yet in pre
        self._rebuild()

    def _rebuild(self):
        self.pre = [0]
        for v in self.a:
            self.pre.append(self.pre[-1] + v)
        self.pending.clear()

    def add(self, i, delta):                    # amortised O(n / K)
        self.a[i] += delta
        self.pending.append((i, delta))
        if len(self.pending) >= self.K:
            self._rebuild()

    def range_sum(self, l, r):                  # O(K)
        s = self.pre[r + 1] - self.pre[l]
        return s + sum(d for i, d in self.pending if l <= i <= r)

This pattern matters when the expensive structure is something without a convenient dynamic version, such as a precomputed distance table or a compressed index. The same idea in a third form splits values instead of positions: in a frequency problem, at most sqrt(n) distinct values can occur more than sqrt(n) times, so heavy values get precomputed tables and light values are handled by brute force. Reordering offline queries by block is the basis of the Mo algorithm.

Choosing between range structures

At n = 1,000,000, log n is about 20 and sqrt(n) is 1,000, so a tree is far ahead on paper. On the summaries a tree handles well, use the tree. Sqrt decomposition earns its place when the summary is awkward to merge, when the operation mix is strongly skewed, or when correctness under time pressure matters more than the last factor of ten.

Failure modes

• Off-by-one at block edges. Wrong partial ranges miss or double-count one element. Test against brute force with tiny block sizes, including B = 1 and B larger than n.
• Forgetting the lazy tag. Reading elements of a block with a pending add without adding it returns stale values. Either add the tag on every read or push it into the elements before touching a partial block.
• Same-block queries. When l and r fall in one block, the three-part split double-counts. Handle that case first, as the code does.
• Wrong block size for the mix. A sqrt(n) block with expensive whole-block work is slow. Balance real costs and benchmark.
• Overflow. Lazy tag times block length can overflow fixed-width integers in C++ or Java. Use 64-bit values.

What to do next

• Implement SqrtDecomp in your main language and test it against a brute-force list on random operations.
• Add a range-assign operation, which needs a second lazy tag and a rule for combining it with adds.
• Build CountLE and time it against a merge sort tree for n = 200,000.
• Benchmark block sizes from sqrt(n) / 4 to 4 sqrt(n) and record where the real optimum lands.
• Solve a distinct-values-in-range problem with the Mo algorithm to see block ordering applied to queries.