DEV Community: Bimal Kshetri

Building an LRU Cache in Python: The Data Structure Behind Every Cache

Bimal Kshetri — Sun, 16 Aug 2026 16:39:35 +0000

A cache is a small box holding copies of expensive things: rows you already fetched from a database, files you already read from disk, results you already computed. Filling it is easy. The whole difficulty is what happens when it is full — something has to go, and throwing out the wrong thing means paying for the expensive work again.

Least recently used, or LRU, is the policy almost everything settles on: discard the entry nobody has touched for the longest. It is a bet that the recent past predicts the near future, and on real workloads that bet pays.

The interesting part is the implementation, because both of the obvious data structures fail. A dict finds a value instantly but has no idea which key is stalest. A list kept in recency order knows exactly which key is stalest, but finding an entry inside it means walking it. What makes every operation O(1) is a hash map and a doubly linked list wired together — and once you have seen it you will recognise it inside CPython's functools.lru_cache, inside MySQL's buffer pool, and behind the eviction setting of every cache server you have configured.

The idea

An LRU cache has a fixed capacity — say three entries — and two operations.

get(key) returns the stored value, or a default if the key is not there.
put(key, value) stores a value, evicting the least recently used entry first if the key is new and the cache is already full.

Both operations count as a use. That single rule is what makes the structure work: after a successful get or any put, that key becomes the most recently used one, and every entry that was ahead of it slides one place closer to the exit. Entries that were already colder than it do not move.

So on every operation the cache must answer three questions, and it must answer all three in constant time or the whole thing is pointless:

Is this key present, and what is its value?
Given a key that is present, mark it as the most recently used.
Which entry is the least recently used, so it can be dropped?

Why a dict on its own fails

A Python dict answers question 1 perfectly and question 3 not at all, because a dict has no notion of "oldest use".

You can bolt one on: store a counter alongside each value and bump it on every access, which handles question 2 as well. But question 3 then means finding the smallest counter, and with no ordering maintained anywhere that is a scan of every entry. One eviction from a 100,000-entry cache reads 100,000 counters.

Why a list in recency order fails

The opposite approach is a plain list of keys, most recently used first. Question 3 becomes free — the answer is the last element. Questions 1 and 2 become expensive, because to move a key to the front you first have to find it, and finding a value in a list means comparing your way along it.

It gets worse than the scan. Deleting from position i in a Python list shifts every later element one slot left, and insert(0, key) shifts every single element one slot right. One read of one key costs three separate linear passes. The code below counts the scan alone: a thousand reads of a thousand-entry cache examine a million list positions, and that total grows with the square of the cache size.

The two structures together

Look again at what each question needs. Question 1 is a hash map's entire job — that is what dicts are for.

Question 2 is "remove this item from the middle of an ordered sequence and reattach it at the front". In an array-backed list that is O(n) because of the shifting. In a doubly linked list it is a handful of pointer assignments, because each node holds a reference to both of its neighbours and can splice itself out without anyone searching for it. The catch is the phrase if you already hold the node. Question 3 is "give me the last node", which the same list answers immediately as long as you keep a reference to the end.

The trick that joins them: the hash map does not store the values, it stores the nodes.

One dict lookup turns a key into the exact node, and from that node the list operations are pure pointer work with no searching at all. The dict never stores any ordering. The list is never searched. Each structure does only the one thing it is fast at, and together they cover all three questions in constant time.

The sentinels

One more piece separates clean code from a nest of special cases.

A linked list normally forces you to ask, on every insert and removal: is this the first node? The last? The only one? Is the list empty? Each answer needs different pointer updates, and each is a place to get it wrong.

The fix is to permanently allocate two nodes holding no data — a head sentinel and a tail sentinel — and keep every real node strictly between them. Now every real node has a neighbour on both sides, the first node is always head.next, the last is always tail.prev, and an empty cache is just the two sentinels pointing at each other. Every branch disappears. Two nodes of wasted memory buy you code with no edge cases in it.

Watching it work

Take a cache with capacity 3 and run seven operations through it. The rightmost column is the recency list, most recently used first — exactly what the linked list holds.

Operation	Returns	Cache, most recent first
`put("A", 1)`	—	`A`
`put("B", 2)`	—	`B A`
`put("C", 3)`	—	`C B A`
`get("A")`	`1`	`A C B`
`put("D", 4)`	— evicts `B`	`D A C`
`get("B")`	`None`	`D A C`
`get("C")`	`3`	`C D A`

Two rows are worth slowing down on.

get("A") moves A from the back to the front. Before the call the list is C B A, so A is the last real node and its neighbours are B and the tail sentinel. The unlink is two assignments: B's next now points at the tail sentinel, and the tail sentinel's prev now points at B. The relink is four more: A's prev becomes the head sentinel, A's next becomes C, the head sentinel's next becomes A, and C's prev becomes A. B and C never move in memory. Only pointers change.

put("D", 4) evicts B. The cache already holds three entries, so something must go, and the coldest entry is whatever sits immediately before the tail sentinel — B, untouched since it was inserted. Notice that get("A") on the previous line is what saved A: without that read, A would have been at the back and A would have gone instead. That is the entire policy, in one step.

The last two rows cover the remaining cases. get("B") misses, because B was evicted, and a miss changes nothing — you cannot promote an entry that is not there. get("C") hits and pulls C to the front, so D and A each slide one place towards eviction.

The code

Start with the version that is correct but wrong: a dict for the values, a list for the order.

class NaiveLRUCache:
    """An LRU cache built from a dict plus a list kept in recency order.

    It behaves correctly. It is still the wrong design, and the counter
    attached to it shows exactly why.
    """

    def __init__(self, capacity: int) -> None:
        self.capacity = capacity
        self.values: dict = {}
        self.recency: list = []  # most recently used first, coldest last
        self.scanned = 0  # list positions examined, so the cost is visible

    def get(self, key, default=None):
        if key not in self.values:
            return default
        position = self.recency.index(key)  # a linear scan — this is the problem
        self.scanned += position + 1
        del self.recency[position]  # shifts every later element one step left
        self.recency.insert(0, key)  # shifts every element one step right
        return self.values[key]

    def put(self, key, value) -> None:
        if key in self.values:
            self.get(key)  # already present, so this is just a recency refresh
        else:
            if len(self.values) == self.capacity:
                del self.values[self.recency.pop()]  # the tail is the coldest key
            self.recency.insert(0, key)
        self.values[key] = value


naive = NaiveLRUCache(capacity=1000)
for number in range(1000):
    naive.put(number, number * number)

# Read every key once, in order. The puts left key 0 at the far end, and each
# read drags the next key you are about to ask for to the far end in turn, so
# every read is a full-length scan.
for number in range(1000):
    naive.get(number)

print(f"1,000 reads of a 1,000-entry cache: {naive.scanned:,} list positions examined")

1,000 reads of a 1,000-entry cache: 1,000,000 list positions examined

A million positions examined for a thousand reads, and that is only the index scan — the del and the insert(0, ...) each shift up to a thousand more pointers on top. Multiply the cache size by a thousand and the per-read cost multiplies by a thousand too.

Here is the real thing. It is longer, and every extra line buys constant time.

from collections.abc import Hashable
from typing import Any


class _Node:
    """One cache entry, and its two neighbours in the recency list.

    The key is stored on the node as well as in the dict. Without it, finding
    the coldest node would not tell you which dict entry to delete.
    """

    __slots__ = ("key", "value", "prev", "next")

    def __init__(self, key: Any = None, value: Any = None) -> None:
        self.key = key
        self.value = value
        self.prev: "_Node | None" = None
        self.next: "_Node | None" = None


class LRUCache:
    """A fixed-capacity cache that discards the least recently used entry.

    Two structures kept in step:

    * a dict mapping each key to its node, for O(1) lookup;
    * a doubly linked list of those nodes in recency order, most recently
      used first, for O(1) reordering and O(1) eviction.

    Both get and put count as a use and move the entry to the front.
    """

    def __init__(self, capacity: int) -> None:
        if capacity < 1:
            raise ValueError("capacity must be at least 1")
        self.capacity = capacity
        self.entries: dict = {}
        # Two sentinel nodes holding no data. They exist so that every real
        # node always has a neighbour on both sides, which removes every
        # "is this the first or last node?" branch from the methods below.
        self.head = _Node()
        self.tail = _Node()
        self.head.next = self.tail
        self.tail.prev = self.head

    def _unlink(self, node: _Node) -> None:
        """Detach a node by making its two neighbours point past it."""
        node.prev.next = node.next
        node.next.prev = node.prev

    def _link_front(self, node: _Node) -> None:
        """Insert a node between the head sentinel and the current first node."""
        first = self.head.next
        node.prev = self.head
        node.next = first
        self.head.next = node
        first.prev = node

    def get(self, key: Hashable, default: Any = None) -> Any:
        node = self.entries.get(key)
        if node is None:
            return default
        self._unlink(node)
        self._link_front(node)
        return node.value

    def put(self, key: Hashable, value: Any) -> None:
        node = self.entries.get(key)
        if node is not None:
            # Overwriting an existing key is a use, not an insertion, so
            # nothing is evicted and the size does not change.
            node.value = value
            self._unlink(node)
            self._link_front(node)
            return
        if len(self.entries) == self.capacity:
            coldest = self.tail.prev  # the node immediately before the tail sentinel
            self._unlink(coldest)
            del self.entries[coldest.key]
        node = _Node(key, value)
        self.entries[key] = node
        self._link_front(node)

    def __len__(self) -> int:
        return len(self.entries)

    def __contains__(self, key: Hashable) -> bool:
        return key in self.entries  # deliberately not a use: no reordering

    def keys_mru_first(self) -> list:
        """Walk the list, for demonstrations. Real callers never need this."""
        keys, node = [], self.head.next
        while node is not self.tail:
            keys.append(node.key)
            node = node.next
        return keys


cache = LRUCache(capacity=3)


def show(label: str) -> None:
    print(f"{label:<20} {' '.join(cache.keys_mru_first())}")


for letter, number in [("A", 1), ("B", 2), ("C", 3)]:
    cache.put(letter, number)
    show(f"put {letter}={number}")

show(f"get A -> {cache.get('A')}")
doomed = cache.keys_mru_first()[-1]
cache.put("D", 4)
show(f"put D=4 (drop {doomed})")
show(f"get B -> {cache.get('B')}")
show(f"get C -> {cache.get('C')}")

put A=1              A
put B=2              B A
put C=3              C B A
get A -> 1           A C B
put D=4 (drop B)     D A C
get B -> None        D A C
get C -> 3           C D A

That is the table from the walkthrough, produced by the code rather than by hand.

How the code maps to the idea

self.entries holds nodes, not values. This is the one line to remember. self.entries.get(key) hands you the node, and from the node you get the value, the two neighbours, and the ability to move it — without touching any other entry.

The sentinels erase the edge cases. _unlink is two assignments with no if in sight, because node.prev and node.next are never None for a real node. _link_front is four assignments for the same reason. Write either method without sentinels and you need branches for "the list is empty", "the node is first", "the node is last" and "the node is the only one".

A promotion is the two methods back to back. _unlink then _link_front: six pointer writes, in that order, whether the cache holds 3 entries or 3 million. They are not always paired — eviction calls _unlink on the coldest node and never relinks it, and a brand-new node goes straight to _link_front with nothing to unlink.

The node stores its own key, which looks redundant until you evict. Eviction starts from the list — self.tail.prev — and ends at the dict, and the only way from a node to its dict entry is the key the node carries. Drop that field and the dict grows for ever while the list stays capped.

put checks for an existing key before it checks capacity. If a full cache evicted first and only then noticed the key was already present, it would have thrown away a live entry to make room for one that needed no room.

get takes a default instead of returning None. If None is a legal value to cache, a bare None return is indistinguishable from a miss and your code recomputes that entry every time. Passing a private sentinel object as the default is the standard way out.

__contains__ does not promote, because checking whether something is cached is not the same as using it. That is a design decision rather than a law, but make it on purpose and document it: a caller who assumes the opposite gets a different eviction pattern. And capacity < 1 is rejected up front, because a capacity of zero breaks the sentinel invariant: the first put sees a full cache, reads self.tail.prev for the coldest node, and gets the head sentinel — the one node that has no prev. _unlink then raises AttributeError on node.prev.next. Rejecting the capacity is cheaper than teaching _unlink to defend itself.

The version you would actually write at work

collections.OrderedDict already contains exactly this doubly linked list, implemented in C. move_to_end does the promotion and popitem does the eviction, so the whole cache collapses to a dozen lines.

from collections import OrderedDict


class OrderedDictLRU:
    """The same cache in a fraction of the lines, using OrderedDict.

    Note the flipped convention: here the most recently used entry sits at
    the *end*, because that is where assigning a new key puts it.
    """

    def __init__(self, capacity: int) -> None:
        self.capacity = capacity
        self.entries: OrderedDict = OrderedDict()

    def get(self, key: Hashable, default: Any = None) -> Any:
        if key not in self.entries:
            return default
        self.entries.move_to_end(key)  # O(1): unlinks and relinks one node
        return self.entries[key]

    def put(self, key: Hashable, value: Any) -> None:
        self.entries[key] = value
        self.entries.move_to_end(key)  # a new key is already last; an old one is not
        if len(self.entries) > self.capacity:
            self.entries.popitem(last=False)  # drop from the cold end

    def keys_mru_first(self) -> list:
        return list(reversed(self.entries))


compact = OrderedDictLRU(capacity=3)
for letter, number in [("A", 1), ("B", 2), ("C", 3)]:
    compact.put(letter, number)
compact.get("A")
compact.put("D", 4)
compact.get("B")
compact.get("C")
print("OrderedDict:", compact.keys_mru_first())
print("from scratch:", cache.keys_mru_first())

OrderedDict: ['C', 'D', 'A']
from scratch: ['C', 'D', 'A']

Same sequence, same final order. The trap is the move_to_end call inside put: assigning to a key that already exists updates the value but does not reorder it, so that line is doing real work even though it looks redundant for new keys. A plain dict can do this too, since deleting a key and reinserting it moves it to the end and dicts have kept insertion order since Python 3.7 — but move_to_end relinks one node instead of leaving a tombstone in the dict's table for a later resize to clean up.

The version you should usually reach for

If what you are caching is the return value of a function, do not build a cache at all. functools.lru_cache is a decorator that wraps any function whose arguments are hashable, and it is the same structure again — a dict of keys to links, plus a circular doubly linked list of those links.

from functools import lru_cache


@lru_cache(maxsize=128)
def fib(n: int) -> int:
    """Naive recursive Fibonacci: exponential uncached, linear with a cache."""
    return n if n < 2 else fib(n - 1) + fib(n - 2)


print(fib(60))
print(fib.cache_info())

1548008755920
CacheInfo(hits=58, misses=61, maxsize=128, currsize=61)

Naive fib(n) makes 2 × F(n+1) − 1 calls, so fib(60) undecorated is about five trillion of them. With the decorator, 61 misses do the real work and 58 hits come straight from the cache. cache_info() is the part people forget: it tells you whether the cache is earning its memory, and a hit count near zero means you are paying for bookkeeping and getting nothing. cache_clear() empties it.

functools.cache is the unbounded variant, identical to lru_cache(maxsize=None). With no ceiling there is nothing to evict, so it drops the linked list entirely and becomes a plain dict lookup — faster, and dangerous in exactly the way you would expect.

from functools import cache


@cache  # shorthand for lru_cache(maxsize=None): unbounded, nothing is ever evicted
def collatz_steps(n: int) -> int:
    """How many halve-or-triple steps it takes to reach 1."""
    if n == 1:
        return 0
    return 1 + collatz_steps(3 * n + 1 if n % 2 else n // 2)


print(max(collatz_steps(start) for start in range(1, 1000)))
print(collatz_steps.cache_info())

178
CacheInfo(hits=998, misses=2228, maxsize=None, currsize=2228)

999 starting values, and the cache holds 2,228 entries — the recursion cached every intermediate value it passed through as well. On a long-running process fed unbounded input, that is a memory leak with a decorator on it. Use functools.cache when the set of possible arguments is small and known, and lru_cache(maxsize=...) when it is not.

Complexity

get: O(1). A hit is one dict lookup, then _unlink (two pointer writes) and _link_front (four). Six writes, every time, with no loop anywhere — the count does not depend on how many entries the cache holds. A miss is one dict lookup and nothing else.

put: O(1). Existing key: one dict lookup, one value assignment, the same six pointer writes. New key into a full cache: one dict lookup, two writes to unlink the coldest node, one dict delete, one node allocation, one dict insert, four writes to link the new node at the front. A fixed amount of work, still no loop.

Eviction: O(1). This is what the tail sentinel buys. self.tail.prev is the least recently used node; nothing is searched, compared or sorted to find it.

The honest caveat is the dict. Python dict operations are O(1) on average, not in the worst case: if every key hashes to the same bucket, lookups degrade to a linear probe sequence, which needs pathological or attacker-chosen keys. Amortisation is involved too — a dict that has seen many inserts and deletes occasionally rebuilds its table, O(capacity) work on one unlucky insert, spread over all the ones that did not trigger it.

Compare the two implementations directly:

Cache capacity	Naive list: work per read	Hash map + linked list
1,000	up to 3,000 element visits	6 pointer writes
100,000	up to 300,000 element visits	6 pointer writes
1,000,000	up to 3,000,000 element visits	6 pointer writes

The naive column is index, then the shift from del, then the shift from insert(0, ...), each up to one full pass.

Space: O(capacity). Not O(number of distinct keys ever seen) — that is the whole selling point of a bounded cache. Each entry costs one _Node, one dict slot, and the key and value objects themselves. The __slots__ declaration matters more than it looks: on the CPython build used here, sys.getsizeof reports 64 bytes for a four-slot node against 48 bytes plus a 104-byte instance dictionary — 152 bytes — for the same class without it. At a capacity of 100,000 entries that one line saves about 8.8 MB. The two sentinels add two more nodes, for ever, regardless of capacity.

When to use it, and when not to

Use LRU when access has temporal locality — the same keys coming back soon after they were last used — and entries cost roughly the same to produce and take roughly the same space. Web sessions, database rows, rendered templates, compiled patterns and decoded images are all recency-friendly.

Do not hand-roll it in Python. Reach for functools.lru_cache if you are caching function results, and OrderedDict if you need a key-value store with your own hooks. Write the node-and-sentinel version only when you need something the standard library will not give you: an eviction callback, per-entry time-to-live, a budget in bytes rather than entries, or a policy that is not quite LRU. Otherwise the C implementations win on speed and on lines of code.

Do not use LRU for scan-shaped workloads. Its one spectacular failure mode is cycling repeatedly through slightly more distinct keys than the cache can hold. With capacity 1,000 and a loop over 1,001 keys, every request evicts the entry you are about to ask for next and the hit rate is exactly zero. Random eviction on that same loop hits about 99% of the time, because a randomly chosen victim is almost never the key you need next — LRU is beaten here by picking without looking. This is not hypothetical; it is what a full table scan does to a database cache. Look instead at MRU, at segmented LRU, or at admission policies that refuse to cache a key on its first sighting.

Do not put an in-process cache in front of shared state you do not own. Eight worker processes means eight independent caches, eight copies of the memory, and eight chances to serve a value another worker has already invalidated. Shared caching needs a shared cache: Redis or memcached.

This implementation is not thread-safe. Two threads interleaving inside _unlink corrupt the list. Every operation must hold one lock — CPython's own pure-Python lru_cache allocates an RLock for exactly this reason, with the comment "because linkedlist updates aren't threadsafe".

Other eviction policies, honestly

LRU is a default, not a law. Each alternative beats it on some real workload.

FIFO evicts whatever was inserted longest ago and never promotes on a read, so there is no bookkeeping on the hot path at all — just a queue. Measurably worse than LRU on most workloads, and measurably cheaper. CPython's re module takes this deal: compiled patterns live in a plain dict capped at 512 entries, and when it fills, the oldest inserted entry is deleted. Insertion-ordered dicts make that a one-liner, and pattern reuse is not recency-shaped enough to justify more.
LFU evicts the least frequently used entry. It beats LRU when popularity is stable and loses badly when it shifts, because an item that was hugely popular last week keeps its enormous count and squats in the cache for ever. Practical LFU implementations therefore decay their counters over time.
Random picks a victim at random. Zero metadata, zero per-access work, immune to the scan pathology, and on many workloads within a few percentage points of LRU's hit rate. If bookkeeping is your bottleneck, this is a serious answer.
MRU evicts the most recently used entry. Backwards, until you meet the cyclic scan: looping over 1,001 keys with room for 1,000, MRU keeps throwing away the same slot and hits on about 99% of requests where LRU hits on none. Random eviction lands within a point of it on that same loop, so MRU is the specialist here, not the only escape.
Belady's optimal algorithm evicts the entry whose next use is furthest in the future. It needs the future, so it cannot be implemented — but it can be computed afterwards on a recorded trace, which makes it the yardstick every other policy is measured against.
Hybrids are what modern systems ship: ARC, 2Q, segmented LRU and W-TinyLFU (used by the Caffeine cache library for Java) all mix recency with frequency and add scan resistance.

Where it shows up in the real world

functools.lru_cache is this post's data structure in the standard library, and memoising an expensive pure function is what most Python programmers use it for. Read functools.py and you will find a dict mapping keys to links, where each link is a four-element list holding prev, next, key and result, joined into a circular doubly linked list. The sentinel trick is there too, in a slicker form: a single root list initialised to point at itself, so the newest entry is always root[PREV] and the coldest is always root[NEXT]. On eviction it reuses the evicted link object as the new root rather than allocating a fresh one. Most CPython builds use a C implementation of the same design and fall back to that Python version when it is unavailable.

Redis exposes LRU as maxmemory-policy allkeys-lru — and does not implement exact LRU. A true recency list would cost two extra pointers on every key plus a list update on every access, across millions of keys. Instead each object header carries a 24-bit clock, and on eviction Redis samples a few random keys (maxmemory-samples, default 5), evicts the oldest of the sample, and keeps a pool of good candidates between rounds. The Redis documentation publishes hit-rate comparisons showing that sampling 10 keys lands very close to exact LRU. Redis 4.0 added allkeys-lfu beside it. The honest lesson: the exact structure is what you learn, and an approximation is what ships when memory per key is the constraint.

MySQL's InnoDB buffer pool splits its LRU list into a young and an old sublist, with the boundary 3/8 of the way down by default (innodb_old_blocks_pct is 37). Newly read pages enter at that midpoint rather than at the head, so a one-off full table scan fills the old sublist and drains out again without displacing the hot pages — a deliberate patch for the scan pathology above.

CPU caches approximate as well. Exact LRU across a 16-way set-associative cache needs enough state to encode an ordering of 16 lines, which is 45 bits per set; the tree-based pseudo-LRU used in real silicon needs 15. Multiply by every set in every cache in every core and the reason is obvious. Operating systems do the same thing one level up: Linux keeps active and inactive page lists with reference bits rather than a strict ordering, and PostgreSQL's shared buffers use a clock sweep with usage counters.

Memcached has used a segmented LRU since version 1.5, with hot, warm and cold segments and a background crawler moving items between them, again because plain LRU spent too much work promoting on every read. Browsers and CDNs evict cached responses under memory and disk pressure, where recency is the dominant signal, usually mixed with object size and remaining freshness.

Common mistakes

Not promoting on get. Leave the _unlink and _link_front pair out of get and you have quietly built a FIFO cache. It works, it passes casual tests, and it evicts entries you are using constantly. Reading is a use.

Updating a value without moving it. With OrderedDict, self.entries[key] = value on an existing key updates the value and leaves the order untouched. The unconditional move_to_end after the assignment is what makes it correct.

Evicting from the list but not the dict. The dict is the only thing keeping the evicted node alive, so forgetting del self.entries[coldest.key] gives you a dict that grows without limit and a len() that no longer matches the list. This is the bug the key field on the node exists to prevent.

Evicting before checking whether the key is present. A put on an existing key needs no space, so evicting first throws away a live entry for nothing and leaves the cache one entry short.

Returning None to mean "not cached". If None is a value you legitimately store, every read of it counts as a miss and the expensive work runs every time. Use an explicit sentinel object as the default.

Assuming functools.lru_cache is free on methods. The cache lives on the function object, shared by every instance, and each entry holds a reference to the self it was called with — so a decorated method keeps its instance alive for as long as the entry survives. For per-instance caching use functools.cached_property, or hold the cache on the instance.

Caching on unhashable arguments. lru_cache builds its key from the arguments, so a list or dict argument raises TypeError. Convert to a tuple or a frozenset at the boundary.

Practice

Add a peek(key) method that returns a value without promoting it, and show that keys_mru_first() is unchanged after calling it.
Add an on_evict callback, invoked with the key and value of every entry the cache drops, and use it to count evictions over a workload.
Add resize(new_capacity) that evicts from the tail until the cache fits, and returns how many entries it dropped.
Run one fixed sequence of key requests through your LRU cache and through a FIFO variant (identical, but get does not promote), and compare the hit counts. Then find a sequence where FIFO wins.
Implement an LFU cache that evicts the least frequently used entry, breaking ties by recency, and find one request pattern where it beats LRU and one where it is far worse.

Summary

An LRU cache is the clearest example in this series of two data structures covering each other's weakness. The hash map cannot order anything; the linked list cannot search. Point the map at the list's nodes and every operation — lookup, promotion, eviction — becomes a fixed number of pointer writes, whatever the capacity. The sentinels keep that code free of branches, and the key stored on each node keeps the two structures in step.


Difficulty	Hard
`get`	O(1) average — one dict lookup, then six pointer writes
`put`	O(1) average — the same, plus at most one eviction
Worst case	O(n) — only if the dict degenerates on pathological keys
Eviction	O(1) — always the node before the tail sentinel
Space	O(capacity) — one node and one dict slot per entry
Data structure	Hash map plus a doubly linked list with two sentinels
Ordering	By recency, most recently used first
Thread-safe	No as written — wrap every operation in one lock
Use it when	Access repeats, and memory has a hard ceiling
Avoid it when	The workload cycles through more keys than the cache holds
Real-world use	`functools.lru_cache`, Redis, the InnoDB buffer pool, CPU caches (approximated)
Python equivalent	`functools.lru_cache` / `functools.cache`, or `OrderedDict.move_to_end`

Write it once from scratch so the mechanism is yours, then use functools.lru_cache for the rest of your career. And keep the failure mode in mind: the moment your access pattern turns into a scan, the most-recently-used entry is the one you should be throwing away.

Keep reading

Hash Tables in Python — where the O(1) lookup half of this structure comes from, and what "average case" is hiding.
Linked Lists in Python — singly, doubly and circular lists built from scratch, including the sentinel trick used above.
Queues and Deques in Python — collections.deque is the same doubly linked list, written in C.
Dynamic Programming Explained — memoisation is the single biggest use of lru_cache, and this is why it works.
Big O Notation — the counting arguments behind every complexity claim here.

Prefix Sums in Python: Answering Range Queries in Constant Time

Bimal Kshetri — Sun, 16 Aug 2026 16:39:29 +0000

You have an array of a million daily sales figures and a dashboard that keeps asking "what did we take between day 12,000 and day 40,000?". Answering that by adding up 28,000 numbers is 28,000 additions per question. Ask forty different questions and you have done a million additions, most of them re-adding the same values you already added a moment ago.

Prefix sums remove that repetition with one move: compute every cumulative total once, up front, then answer any range sum with a single subtraction. Building the table costs O(n). Every query afterwards costs O(1) — two array reads and a minus sign — and that cost does not change whether the range covers 3 elements or 3 million.

The idea takes thirty seconds to write. The off-by-one takes most people an hour to stop getting wrong, so that is where this post spends its time first. After that, the same trick extends in three directions: rectangles in a 2D grid, range updates instead of range queries, and the hash-map version that counts subarrays summing to a target in linear time.

The idea

Take an array, and write down a second array that answers one question at every position: what is the total of everything before here?

That second array is the prefix sum array. For nums = [3, 1, 4, 1, 5, 9, 2, 6] it looks like this:

index      0    1    2    3    4    5    6    7    8
nums            3    1    4    1    5    9    2    6
prefix     0    3    4    8    9   14   23   25   31

Read it a column at a time. prefix[0] is 0, because the total of no numbers is zero. prefix[1] is 3, the total of the first one element. prefix[5] is 14, the total of the first five: 3 + 1 + 4 + 1 + 5. The rule that builds the whole thing is one line — prefix[i + 1] = prefix[i] + nums[i] — so each entry is its left neighbour plus one new value.

Two details look pedantic and are not.

The prefix array holds n + 1 entries, not n. Eight input numbers give nine cumulative totals. The extra one is the leading zero.

prefix[i] does not include nums[i]. It is the total of everything strictly to the left of index i. So a prefix entry describes the gap before an element rather than the element itself, which is exactly why the array needs one more slot than the input: there are nine gaps around eight numbers.

Now the payoff. The sum of nums[left..right], inclusive at both ends, is:

prefix[right + 1] - prefix[left]

prefix[right + 1] counts everything up to and including nums[right]; prefix[left] counts everything strictly before nums[left]. Both start from the front of the array, so subtracting cancels the shared front section and leaves exactly the elements from left to right.

Why the leading zero matters

Suppose you skip the sentinel and build a prefix array the same length as the input, where prefix[i] is the total of the first i + 1 elements. The query then reads prefix[right] - prefix[left - 1], and that expression is broken at left = 0.

In most languages prefix[-1] throws. In Python it does something worse: it returns the last element, which is the total of the entire array. The range sum comes back wrong, no exception is raised, and nothing points at the cause. You then patch it with if left == 0 and carry that special case everywhere.

The leading zero deletes the special case instead of handling it. prefix[0] is the sum of an empty range, which genuinely is zero, so left = 0 needs no different treatment from left = 5. That is the shape of a sentinel: one extra entry so the ordinary rule covers the boundary too.

Watching it work

Take the same array and ask for nums[2..5] — the values 4, 1, 5, 9, which add up to 19.

Look up two numbers:

prefix[6] = 23. That is 3 + 1 + 4 + 1 + 5 + 9, the first six elements, which is everything up to and including nums[5].
prefix[2] = 4. That is 3 + 1, everything strictly before nums[2].

Subtract: 23 − 4 = 19. The 3 + 1 at the front appears in both totals, so it cancels, and what is left is 4 + 1 + 5 + 9. One subtraction, no loop.

Now the awkward case, nums[0..3]: values 3, 1, 4, 1, total 9. The query is prefix[4] - prefix[0], which is 9 − 0 = 9. Nothing special happened. The zero at the front absorbed the boundary, exactly as designed.

One more, the single-element range nums[7..7]: prefix[8] - prefix[7] is 31 − 25 = 6, which is nums[7]. It works because adjacent prefix entries differ by exactly the element between them.

The code

The naive version first, so there is something to measure against.

def range_sum_naive(nums: list[int], left: int, right: int) -> int:
    """Add up nums[left..right] inclusive, one element at a time."""
    total = 0
    for index in range(left, right + 1):
        total += nums[index]
    return total


nums = [3, 1, 4, 1, 5, 9, 2, 6]
print(range_sum_naive(nums, 2, 5))
print(range_sum_naive(nums, 0, 3))

19
9

Correct, and the cost of every call is proportional to the width of the range. Here is the version that pays that cost once.

def build_prefix(nums: list[int]) -> list[int]:
    """Return prefix[i] = the sum of the first i values of nums.

    The result has len(nums) + 1 entries. prefix[0] is 0 — the sum of no
    values at all — and that sentinel is what makes every query a single
    subtraction with no special case for a range starting at index 0.
    """
    prefix = [0] * (len(nums) + 1)
    for index, value in enumerate(nums):
        prefix[index + 1] = prefix[index] + value
    return prefix


def range_sum(prefix: list[int], left: int, right: int) -> int:
    """Sum of nums[left..right] inclusive, from a prefix array."""
    return prefix[right + 1] - prefix[left]


prefix = build_prefix(nums)
print(prefix)
print(range_sum(prefix, 2, 5))
print(range_sum(prefix, 0, 3))
print(range_sum(prefix, 7, 7))
print(range_sum(prefix, 0, 7))

[0, 3, 4, 8, 9, 14, 23, 25, 31]
19
9
6
31

Those are the three cases from the walkthrough plus the whole array, and none needed a branch.

You rarely have to write build_prefix yourself. The standard library has had this since Python 3.2, and the initial argument that supplies the sentinel since 3.8:

from itertools import accumulate

print(list(accumulate(nums)))
print(list(accumulate(nums, initial=0)))
print(list(accumulate(nums, initial=0)) == build_prefix(nums))

[3, 4, 8, 9, 14, 23, 25, 31]
[0, 3, 4, 8, 9, 14, 23, 25, 31]
True

accumulate is a C-level loop, so it is faster than the Python for loop above and it is the version to reach for at work. Write build_prefix by hand only while you are learning it.

How the code maps to the idea

[0] * (len(nums) + 1) is the n + 1 rule made concrete. Allocating the exact size up front also means the loop only ever assigns to slots that already exist, so no append and no resizing.

prefix[index + 1] = prefix[index] + value is the one-line rule from the walkthrough. The + 1 on the left is what shifts the whole array right by one and leaves the sentinel untouched at index 0. Every entry is computed from the entry immediately before it, which is why one pass is enough — the work is never repeated.

prefix[right + 1] - prefix[left] is asymmetric on purpose, and that asymmetry is the whole off-by-one. The left end is exclusive in the prefix array and inclusive in the range, so it needs no adjustment; the right end is inclusive in the range, so it needs the + 1 to move past it. For half-open ranges — nums[left:right], matching Python slicing — the expression is the symmetric prefix[right] - prefix[left]. Pick one convention, put it in the docstring, and never mix the two.

Edge cases fall out of the arithmetic. An empty input gives [0], and there is no valid query to run against it. A range covering everything is prefix[n] - prefix[0], which the zero collapses to just prefix[n].

The values do not have to be positive. Nothing above assumes it — negative numbers and floats both work, though see the float warning further down.

Complexity

Build: O(n). The loop body runs exactly once per input element and does one addition and one store each time. Eight elements, eight additions. A million elements, a million additions. There is no nesting and no repeated work, so the count is exactly n.

Query: O(1). Two list index operations and one subtraction. Python lists are contiguous arrays of pointers, so indexing is a constant-time address calculation, not a walk. Three operations regardless of n and — the interesting part — regardless of how wide the range is. A range of 999,999 elements costs the same as a range of one.

Space: O(n). One extra list of n + 1 integers. That is a real cost, not a rounding error: for 10 million values you are holding a second 10-million-entry list.

There is no best or worst case to separate out here. The build loop runs exactly n iterations whatever the values are — no early exit, no branch that depends on the data — and a query does the same two lookups and one subtraction every time. Best, average and worst are the same bound, which is unusual and is exactly what makes the technique easy to reason about.

The trade only pays off when there are enough queries. Counting it exactly, for q queries whose ranges have total length L:

Approach	Operations
Naive	L additions, worst case q × n
Prefix sums	n additions to build, then q subtractions

queries = [(0, 7), (2, 5), (1, 6), (3, 3), (0, 3)]

for left, right in queries:
    assert range_sum(prefix, left, right) == range_sum_naive(nums, left, right)

naive_additions = sum(right - left + 1 for left, right in queries)
print(f"naive:  {naive_additions} additions")
print(f"prefix: {len(nums)} additions to build + {len(queries)} subtractions "
      f"= {len(nums) + len(queries)}")

naive:  23 additions
prefix: 8 additions to build + 5 subtractions = 13

Eight elements and five queries is already a win, and the gap widens fast. With n = 1,000,000 and 100,000 queries averaging 30,000 elements each, the naive approach performs about 3 billion additions and the prefix approach performs 1.1 million operations — roughly 2,700 times fewer.

The other direction is worth stating too. For a single query over the whole array, the naive loop does n additions while the prefix version does n additions plus an allocation plus a subtraction — strictly worse. This is an amortisation: cheap queries bought with an upfront payment, and the break-even point is when L, the total length of all the queries, exceeds n + q.

Two dimensions: sums over a rectangle

The same idea works on a grid, and the payoff is bigger because a naive rectangle sum is quadratic in the rectangle's side length.

Take this 4 × 4 grid, and ask for the total of the rectangle from row 1 to row 2 and column 1 to column 3:

Build a table P where P[r][c] is the sum of every cell strictly above row r and strictly left of column c — the whole top-left block. Same convention as before: pad with a zero row and a zero column, so P is 5 × 5 for a 4 × 4 grid.

Each entry is built from three neighbours:

P[r+1][c+1] = grid[r][c] + P[r][c+1] + P[r+1][c] - P[r][c]

The block above and the block to the left both contain the block that is above and to the left, so that corner block gets counted twice and has to come off once. This is inclusion-exclusion, and it appears again in the query.

To read a rectangle from (top, left) to (bottom, right) inclusive, take the big block ending at its far corner, cut off the band above, cut off the band to the left, then add back the top-left block you have now removed twice:

P[bottom+1][right+1] - P[top][right+1] - P[bottom+1][left] + P[top][left]

def build_prefix_2d(grid: list[list[int]]) -> list[list[int]]:
    """Return prefix[r][c] = the sum of the block above and left of (r, c).

    Both dimensions get a leading zero row/column, so no query needs a
    boundary check.
    """
    rows, cols = len(grid), len(grid[0])
    prefix = [[0] * (cols + 1) for _ in range(rows + 1)]
    for row in range(rows):
        for col in range(cols):
            # The block above and the block to the left overlap in the corner
            # block, so it is counted twice and has to come off once.
            prefix[row + 1][col + 1] = (
                grid[row][col]
                + prefix[row][col + 1]
                + prefix[row + 1][col]
                - prefix[row][col]
            )
    return prefix


def rectangle_sum(prefix: list[list[int]], top: int, left: int,
                  bottom: int, right: int) -> int:
    """Sum of the inclusive rectangle from (top, left) to (bottom, right)."""
    return (
        prefix[bottom + 1][right + 1]
        - prefix[top][right + 1]
        - prefix[bottom + 1][left]
        + prefix[top][left]
    )


grid = [
    [3, 0, 1, 4],
    [5, 6, 3, 2],
    [1, 2, 0, 1],
    [4, 1, 0, 1],
]
prefix_2d = build_prefix_2d(grid)
for row in prefix_2d:
    print(row)
print(rectangle_sum(prefix_2d, 1, 1, 2, 3))
print(rectangle_sum(prefix_2d, 0, 0, 3, 3))
print(rectangle_sum(prefix_2d, 2, 2, 2, 2))

[0, 0, 0, 0, 0]
[0, 3, 3, 4, 8]
[0, 8, 14, 18, 24]
[0, 9, 17, 21, 28]
[0, 13, 22, 26, 34]
14
34
0

Check the first answer by hand: the rectangle holds 6, 3, 2 on row 1 and 2, 0, 1 on row 2, which is 14. The formula gives 28 − 8 − 9 + 3 = 14. The second answer, 34, is the total of all sixteen cells; the third is the single cell at (2, 2), which is 0.

Cost. Building touches each of the rows × cols cells once and does three additions and subtractions per cell, so the build is O(rows × cols) — linear in the size of the grid, which is the least possible since every cell has to be read at least once. Each query is four lookups and three arithmetic operations: O(1), independent of the rectangle's area. Space is O(rows × cols) for the second table, and that is the real limit — a 10,000 × 10,000 grid needs a 100-million-entry prefix table.

The inverse problem: difference arrays

Prefix sums make many range reads cheap. Turn the problem around — many range writes, one read at the end — and the same relationship solves it backwards.

Say you have 8 slots, all starting at zero, and three updates: add 2 to indices 1 through 4, add 3 to indices 0 through 2, add 1 to indices 5 through 7. Done naively each update loops over its range, so k updates of average width w cost k × w writes.

Instead, record only the changes. Keep an array diff where diff[i] is how much the value jumps between index i - 1 and index i. Adding amount across [left, right] needs exactly two jumps: a step up at left, and a step back down just past right.

diff[left]      += amount     # from here on, everything is `amount` higher
diff[right + 1] -= amount     # past the right edge, take it back

Two writes per update, whatever the width. Then a single running sum over diff reconstructs the real values — and a running sum is a prefix sum, which is the whole point: the difference array is the inverse operation, so summing it undoes it.

Note that diff needs n + 1 slots again, this time for the opposite reason: an update whose right edge is the last index writes to diff[n], one past the end of the real data. That slot exists only to be ignored.

def build_difference(length: int, updates: list[tuple[int, int, int]]) -> list[int]:
    """Record each range update as two boundary changes, in O(1) per update."""
    diff = [0] * (length + 1)
    for left, right, amount in updates:
        diff[left] += amount        # from here on, every value is `amount` higher
        diff[right + 1] -= amount   # past the right edge, give it back
    return diff


def resolve(diff: list[int], length: int) -> list[int]:
    """Turn a difference array into the real values with one running sum."""
    return list(accumulate(diff[:length]))


def apply_updates_naive(length: int, updates: list[tuple[int, int, int]]) -> list[int]:
    """The obvious version: touch every index inside every range."""
    values = [0] * length
    for left, right, amount in updates:
        for index in range(left, right + 1):
            values[index] += amount
    return values


updates = [(1, 4, 2), (0, 2, 3), (5, 7, 1)]
diff = build_difference(8, updates)
print(diff)
print(resolve(diff, 8))
print(resolve(diff, 8) == apply_updates_naive(8, updates))

[3, 2, 0, -3, 0, -1, 0, 0, -1]
[3, 5, 5, 2, 2, 1, 1, 1]
True

Follow the running sum across [3, 2, 0, -3, 0, -1, 0, 0]: 3, 5, 5, 2, 2, 1, 1, 1. Adding the three updates by hand gives [3, 3, 3, 0, 0, 0, 0, 0] plus [0, 2, 2, 2, 2, 0, 0, 0] plus [0, 0, 0, 0, 0, 1, 1, 1], and those totals match — which is what the True on the last line asserts.

Cost. k updates at two writes each is O(k), then one pass to resolve is O(n): O(n + k) total, against O(n × k) in the worst case for the naive loop. The catch is that you cannot read a value until you have resolved, so this only works when every update arrives before any read. That restriction is exactly what makes it cheap.

The interview classic: subarray sum equals k

"How many contiguous subarrays add up to k?" is the problem that makes prefix sums click, because the naive answer is O(n²) and the fix is one hash map.

Every subarray nums[i..j] has sum prefix[j + 1] - prefix[i]. So asking for subarrays that sum to k is asking for pairs of prefix entries that differ by exactly k, with the earlier one first. Rewrite the condition:

prefix[j + 1] - prefix[i] = k, therefore prefix[i] = prefix[j + 1] - k

So walk the array once, keeping a running prefix sum. At each step, the number of subarrays ending here is the number of earlier prefix values equal to running - k. Keep a count of every prefix value seen so far in a dictionary, and that lookup is O(1) on average.

Try it on [3, 4, 7, 2, -3, 1, 4] with k = 7. The running prefix sums, sentinel included, are [0, 3, 7, 14, 16, 13, 14, 18], and exactly three ordered pairs differ by 7:

Those three pairs are the three subarrays: [3, 4], [7], and [7, 2, -3, 1].

from collections import defaultdict


def count_subarrays_with_sum(nums: list[int], target: int) -> int:
    """Count the contiguous subarrays of nums whose values add up to target."""
    counts: defaultdict[int, int] = defaultdict(int)
    counts[0] = 1        # the empty prefix, so a whole prefix can match on its own
    running = 0
    found = 0
    for value in nums:
        running += value
        # Every earlier prefix equal to running - target closes a subarray here.
        found += counts[running - target]
        counts[running] += 1
    return found


sample = [3, 4, 7, 2, -3, 1, 4]
print(count_subarrays_with_sum(sample, 7))
print(count_subarrays_with_sum([1, 1, 1], 2))
print(count_subarrays_with_sum([0, 0, 0], 0))

3
2
6

[1, 1, 1] with k = 2 has two answers, the first pair and the second. [0, 0, 0] with k = 0 has six, because all six subarrays sum to zero — three of length 1, two of length 2, one of length 3. That last case is what catches people who try to count with pointers instead.

Watch the counter run:

seen: defaultdict[int, int] = defaultdict(int)
seen[0] = 1
running = 0
for index, value in enumerate(sample):
    running += value
    print(f"i={index}  value={value:>3}  prefix={running:>3}  "
          f"want {running - 7:>3}  seen {seen[running - 7]} time(s)")
    seen[running] += 1

i=0  value=  3  prefix=  3  want  -4  seen 0 time(s)
i=1  value=  4  prefix=  7  want   0  seen 1 time(s)
i=2  value=  7  prefix= 14  want   7  seen 1 time(s)
i=3  value=  2  prefix= 16  want   9  seen 0 time(s)
i=4  value= -3  prefix= 13  want   6  seen 0 time(s)
i=5  value=  1  prefix= 14  want   7  seen 1 time(s)
i=6  value=  4  prefix= 18  want  11  seen 0 time(s)

The hit at i=1 is the one that needs counts[0] = 1. There, running is 7 and the algorithm looks for an earlier prefix of 0 — the sentinel, the empty prefix before the array starts. Without seeding the map you lose every subarray beginning at index 0, and the bug is invisible on test data whose answers all start later.

Cost. One pass, with one dictionary lookup and one dictionary write per element. Python dictionaries give O(1) average lookup, so the whole thing is O(n) average time and O(n) space. The map holds at most 2n + 1 entries: each step stores the prefix it just computed, and reading counts[running - target] on a defaultdict also inserts a zero for the value it failed to find.

Average is not a guarantee, and it is worth knowing why. In CPython a non-negative integer hashes to itself reduced modulo 2**61 - 1, and that is not randomised, so prefix sums chosen to share a hash really do drive each dictionary operation towards O(n) and the whole scan towards O(n²). Ordinary data never does that, but quote the O(n) as an average, not as a worst case.

This matters because the sliding window technique solves the same problem in O(1) space when every value is positive, since growing the window then always grows the sum. Add one negative number and that monotonicity is gone and the window breaks. Prefix sums plus a hash map do not care about signs at all.

When to use it, and when not to

Use prefix sums when the data is fixed and you will query it many times. Log analysis, dashboard ranges over historical data, sensor readings, precomputed lookup tables. The pattern to recognise: "given an array and q queries, each asking for a range".

Use a difference array when the updates all come before the reads — bulk range increments, interval booking counts, timeline aggregation.

Do not use it when the values change between queries. Updating nums[i] invalidates every prefix entry from i + 1 onwards, so a single element change costs O(n) to repair. If reads and writes interleave, reach for a Fenwick tree (binary indexed tree) or a segment tree: both give O(log n) point updates and O(log n) range sums. A prefix array is the degenerate case of those structures where you have decided updates never happen and bought O(1) queries with that promise.

Do not use it for a handful of queries. One or two range sums over a large array are cheaper computed directly than by building an n-entry table first.

Be careful with floats. Prefix sums of floating-point numbers accumulate rounding error along the whole array, and then the query subtracts two large, nearly equal totals — catastrophic cancellation, where the leading digits cancel and the error hiding beneath them becomes the answer's leading digits. For an accurate total of floats use math.fsum, which adds them exactly and rounds only once at the end. For range sums of floats, know that the error is proportional to the magnitude of the whole prefix, not just the range you asked about.

Integer overflow is not a Python problem — Python integers grow as needed — but it is the classic bug when you port this to C++, Java or Rust. A prefix sum over 100,000 values of up to a billion reaches 10^14, which overshoots the signed 32-bit maximum by a factor of about 46,000.

Where it shows up in the real world

Integral images in computer vision. A summed-area table is a 2D prefix sum, introduced by Frank Crow in 1984 for texture mapping and made famous by the Viola-Jones face detector in 2001. Every feature Viola-Jones scores is a difference of rectangle sums, and the integral image makes each of those rectangle sums four table lookups whatever the rectangle's size. Its detector holds a few thousand features arranged as a cascade of stages, and because most windows are rejected by the first stages only a handful of features are ever evaluated on a given window. A rectangle sum that costs the same at every scale is what made real-time face detection possible on 2001 hardware. OpenCV exposes it as cv2.integral to this day.

Running totals in databases and dataframes. SUM(amount) OVER (ORDER BY day ROWS UNBOUNDED PRECEDING) in PostgreSQL, SQL Server or BigQuery is a prefix sum computed by the engine in one pass. So are pandas.Series.cumsum and numpy.cumsum.

Weighted random choice in the Python standard library. random.choices(population, weights=...) calls itertools.accumulate on the weights to build a prefix array, then uses bisect to binary-search it for each pick. That is a prefix sum turning a weighted draw into a binary search over cumulative totals.

Parallel computing. Prefix sum is one of the fundamental parallel primitives, known there as scan. CUDA's Thrust library ships inclusive_scan and exclusive_scan, and they are the standard way to compute output offsets when many threads write variable-length results into one buffer: each thread's write position is the prefix sum of the sizes before it.

Common mistakes

Dropping the + 1 in the query. prefix[right] - prefix[left] gives you nums[left..right-1] — one element short at the right end. Fix: decide whether your range is inclusive or half-open, write it in the docstring, and test a single-element range, which fails loudly under both mistakes.

Building the prefix array without the sentinel. Then left = 0 needs prefix[-1], and in Python that silently reads the last element instead of raising. Always allocate n + 1 and start at zero.

Making the difference array length n. An update ending at the last index writes diff[n], so a length-n array raises IndexError on exactly the ranges that reach the end — often the last case anyone tests.

Forgetting the corner in 2D. Both the build step and the rectangle query need a corner term, and the signs are opposite: the build subtracts the block it counted twice, the query adds it back. Leave it out of the query and the answer is too small by the top-left block, which is correct-looking whenever top or left is 0 and wrong everywhere else. Test a rectangle that touches neither edge.

Missing counts[0] = 1 in the subarray counter. Every subarray that starts at index 0 is silently skipped. Test with [k] — a one-element array whose only element is the target, whose answer must be 1.

Mutating the input after building the prefix array. The prefix array is a snapshot. It goes stale the moment nums changes, and nothing warns you. If the data is mutable, either rebuild or use a Fenwick tree.

Practice

Given an array and a list of (left, right) queries, return the average of each range in O(1) per query.
Find an index where the sum of everything to its left equals the sum of everything to its right, in one pass over the array.
Given an array of only 0s and 1s, find the longest subarray with equal counts of each — map every 0 to −1 and look for two prefix sums that are equal.
Given a list of flight bookings, each (first_flight, last_flight, seats), return the total seats booked on every flight using a difference array.
Count the submatrices of a grid whose values sum to a target, by fixing a pair of rows and running the 1D subarray counter along the columns between them.

Summary

Prefix sums are the cheapest possible trade: one linear pass up front buys constant-time range sums forever after, as long as the data stops changing. Store n + 1 cumulative totals with a zero at the front, and every query is prefix[right + 1] - prefix[left] with no boundary cases to remember. The same relationship read backwards gives difference arrays, and the same relationship read as "pairs of prefix values k apart" gives the linear-time subarray counter.


Difficulty	Easy
Build time	O(n) — exactly one addition per element
Query time	O(1) — two lookups and one subtraction, any range width
Update time	O(n) — changing one value invalidates every later prefix
Space	O(n) — a second array of n + 1 totals
2D build / query	O(rows × cols) / O(1) — four lookups per rectangle
Difference array	O(1) per range update, O(n) once to resolve
Subarray counting	O(n) average time and O(n) space with a hash map
Data structure	List / array with O(1) indexing
Use it when	The array is fixed and total query length exceeds n
Avoid it when	Values change between queries — use a Fenwick or segment tree
Real-world use	Integral images in face detection, SQL running totals, `random.choices`
Python equivalent	`itertools.accumulate(nums, initial=0)`

Learn the n + 1 convention once and it will keep paying out, because the same off-by-one shows up in difference arrays, in 2D tables, and in every sliding-window variant that tracks a running total.

Keep reading

The Sliding Window Technique — the other way to make subarray problems linear, and where it stops working.
The Two Pointers Technique — the third member of this family of linear-scan tricks.
Hash Tables in Python — why the dictionary lookup in the subarray counter is O(1), from the inside.
Arrays and Dynamic Arrays — why indexing a Python list is constant time in the first place.
Big O Notation — the counting arguments used above, done properly.

The Sliding Window Technique in Python: Subarray Problems Made Linear

Bimal Kshetri — Sun, 16 Aug 2026 16:38:57 +0000

A list of 100,000 numbers contains 5,000,050,000 contiguous runs. If your first instinct for "find the best subarray" is to look at each of them, you have written an algorithm that will not finish this afternoon. A sliding window answers the same question in fewer than 200,000 steps, because it never touches any element more than twice: once when it joins the window, once when it leaves.

A window is just a contiguous run, held by two indices: left and right. The technique is a rule about how those indices are allowed to move. The right edge steps forward and one element joins the window. The left edge steps forward and one element leaves. Nothing else ever happens, neither edge ever goes backwards, and the quantity you care about — a sum, a character count, a set of distinct values — is patched up as elements join and leave instead of being recomputed from the elements still inside.

Two things about this pattern get skipped almost everywhere, and both are here. The first is why the variable-size version is still linear despite having a while loop nested inside a for loop, which looks quadratic and is not. The second is the precondition: sliding windows are only correct when shrinking from the left can actually repair an invalid window. Break that condition — an array containing negative numbers is the standard way to break it — and the algorithm does not crash, it just returns the wrong answer.

The idea

Fixed-size windows: add one, remove one

Start with the simplest version, where the width is handed to you. Given the readings [2, 1, 5, 1, 3, 2], what is the largest sum of 3 consecutive values?

The obvious approach adds up each window from scratch. There are n - k + 1 windows of width k in a list of n values, and each one costs k reads, so the bill is (n - k + 1) * k reads. For n of 100,000 and k of 1,000 that is 99,001,000 reads to answer one question about a list you could scan nine hundred and ninety times over in the same effort.

The waste is obvious once you line two neighbouring windows up. The window at indices 1 to 3 shares indices 1 and 2 with the window at indices 0 to 2. Only one value left and one value arrived. So carry the sum forward:

subtract the value that just fell off the left edge,
add the value that just arrived at the right edge.

Two operations per step, regardless of how wide the window is.

Variable-size windows: expand, then shrink

Most real problems do not tell you the width. Instead they give you a rule the window must obey — no repeated characters, at most 2 distinct values, a sum no greater than 10 — and ask for the longest (or shortest) run that obeys it.

The loop is always the same three steps:

Move right forward by one. The new element joins the window.
While the window breaks the rule, move left forward by one, evicting the leftmost element.
The window is valid again. Record it if it beats the best seen so far.

Step 2 is where the correctness lives, and it rests on one assumption that is worth stating out loud because almost nobody does:

If a window is invalid, every window that contains it is invalid too.

That is the monotonicity precondition, and it does two separate jobs. For a fixed right edge it makes validity a threshold: there is some index where the window becomes valid, everything at or after it is valid, and everything before it is not — so shrinking is guaranteed to reach that threshold rather than overshooting past it. For a fixed left edge it says an invalid window cannot be rescued by taking in more on the right, and that is what lets left stay put: once a left edge has been evicted, no window starting there will ever be valid again, at any later right edge. Both halves are needed, and both fall out of the one sentence above.

Check it on "no repeated characters": if a run contains two copies of the letter a, then any longer run containing it also contains both copies, so it is invalid too. It holds. Check it on "sum at most 10" with non-negative values: adding values on either side can only make the sum bigger, never smaller, so a run that is already over the limit stays over. It holds — as long as the values are non-negative. Hold that thought.

Watching it work

The fixed window, by hand

Take [2, 1, 5, 1, 3, 2] with k of 3. The first window is added up in full; every window after that is one subtraction and one addition.

Window	Values	Sum
indices 0 to 2	2, 1, 5	8, added in full
indices 1 to 3	1, 5, 1	8 − 2 + 1 = 7
indices 2 to 4	5, 1, 3	7 − 1 + 3 = 9
indices 3 to 5	1, 3, 2	9 − 5 + 2 = 6

The answer is 9. The whole list was read once, plus three extra reads for the values leaving — 9 reads against the brute force's 12. At six elements that is a shrug. At 100,000 elements with k of 1,000 it is 199,000 reads against 99,001,000.

The variable window, by hand

Now the harder shape. Find the length of the longest substring of "abcabcbb" with no repeated character.

Keep a dict last_seen mapping each character to the most recent index it appeared at, and a start index for the left edge. For each character you read, if it was last seen inside the current window, the left edge jumps to just past that earlier copy. That is a shrink of several positions in one move, which is fine — left only ever goes forward.

i=0, read a. Window "a", length 1.
i=1, read b. Window "ab", length 2.
i=2, read c. Window "abc", length 3. Best so far.
i=3, read a. Last seen at index 0, which is inside the window, so start jumps to 1. Window "bca", length 3.
i=4, read b. Last seen at index 1, inside. start jumps to 2. Window "cab".
i=5, read c. Last seen at index 2, inside. start jumps to 3. Window "abc".
i=6, read b. Last seen at index 4, inside. start jumps to 5. Window "cb", length 2.
i=7, read b. Last seen at index 6, inside. start jumps to 7. Window "b", length 1.

The answer is 3. Count the movement: the right edge took eight steps, and the left edge moved on five of those steps, travelling 0 to 1 to 2 to 3 to 5 to 7 — seven positions in total. Fifteen moves for an eight-character string, against a ceiling of 2n = 16. No substring was ever built or rescanned.

The code

The naive fixed-size version first, so there is something honest to compare against.

def max_sum_naive(values: list[int], k: int) -> int:
    """Largest sum of k consecutive values, recomputing every window in full."""
    if not 1 <= k <= len(values):
        raise ValueError("k must be between 1 and len(values)")

    best = sum(values[0:k])
    for start in range(1, len(values) - k + 1):
        best = max(best, sum(values[start:start + k]))
    return best


readings = [2, 1, 5, 1, 3, 2]
print(max_sum_naive(readings, 3))

Correct, and it re-reads k values every step for no reason. The window version keeps the running sum instead.

def max_sum_window(values: list[int], k: int) -> int:
    """Largest sum of k consecutive values, in a single pass.

    The first window is added up in full. After that each step adds the value
    entering on the right and subtracts the one leaving on the left, so keeping
    the sum correct costs two operations no matter how wide the window is.
    """
    if not 1 <= k <= len(values):
        raise ValueError("k must be between 1 and len(values)")

    # Summed by index rather than sum(values[:k]) so that no k-element slice is
    # ever built; the whole point is that this function holds nothing but a total.
    window_sum = sum(values[index] for index in range(k))
    best = window_sum

    for right in range(k, len(values)):
        window_sum += values[right] - values[right - k]
        best = max(best, window_sum)

    return best


print(max_sum_window(readings, 3))
print(max_sum_window([7], 1), max_sum_window([-4, -1, -9], 2))

9
7 -5

The gap between them is arithmetic, not opinion. The naive version reads (n - k + 1) * k values; the window version reads k for the first window and 2 for every step after.

print(f"{'n':>9} {'k':>7} {'naive reads':>15} {'window reads':>13}")
for size, width in [(1_000, 100), (100_000, 1_000), (1_000_000, 10_000)]:
    naive_reads = (size - width + 1) * width
    window_reads = width + 2 * (size - width)
    print(f"{size:>9,} {width:>7,} {naive_reads:>15,} {window_reads:>13,}")

        n       k     naive reads  window reads
    1,000     100          90,100         1,900
  100,000   1,000      99,001,000       199,000
1,000,000  10,000   9,900,010,000     1,990,000

Now the variable-size version, plus a traced copy that prints the walkthrough above so you can check the hand trace against the machine.

def longest_unique(text: str) -> int:
    """Length of the longest substring of text with no repeated character.

    last_seen maps a character to the most recent index it appeared at. When
    the character entering on the right is already inside the window, the
    window's left edge jumps to just past that earlier copy.
    """
    last_seen: dict[str, int] = {}
    start = 0
    best = 0

    for index, character in enumerate(text):
        previous = last_seen.get(character, -1)
        if previous >= start:
            # The repeat is inside the window, so every window that keeps both
            # copies is invalid. Skip left past the older one in one move.
            start = previous + 1
        last_seen[character] = index
        best = max(best, index - start + 1)

    return best


def longest_unique_traced(text: str) -> int:
    """The same algorithm, printing the window after every character read."""
    last_seen: dict[str, int] = {}
    start = 0
    best = 0
    for index, character in enumerate(text):
        previous = last_seen.get(character, -1)
        moved = ""
        if previous >= start:
            start = previous + 1
            moved = f"repeat at {previous}, start jumps to {start}"
        last_seen[character] = index
        best = max(best, index - start + 1)
        line = (f"i={index} read {character}  window [{start}:{index + 1}]"
                f" = {text[start:index + 1]:<4} best={best}")
        print(line + f"  {moved}" if moved else line)
    return best


print(longest_unique_traced("abcabcbb"))

i=0 read a  window [0:1] = a    best=1
i=1 read b  window [0:2] = ab   best=2
i=2 read c  window [0:3] = abc  best=3
i=3 read a  window [1:4] = bca  best=3  repeat at 0, start jumps to 1
i=4 read b  window [2:5] = cab  best=3  repeat at 1, start jumps to 2
i=5 read c  window [3:6] = abc  best=3  repeat at 2, start jumps to 3
i=6 read b  window [5:7] = cb   best=3  repeat at 4, start jumps to 5
i=7 read b  window [7:8] = b    best=3  repeat at 6, start jumps to 7
3

How the code maps to the idea

values[right] - values[right - k] is the add-one-remove-one update, written as a single expression. When the right edge sits at index right, the window covers the k indices ending there, so the value that just fell out is exactly k positions back. Getting that offset wrong by one is the classic bug in this line, and it does not raise — it silently sums the wrong window.

The first window is summed in full, outside the loop. That is the only place sum() appears. If you see sum() inside the loop, the running total is not running and the algorithm is back to quadratic.

if previous >= start is the guard that makes the left edge one-directional. Without it, start is set to previous + 1 whenever the character has ever been seen, including when that sighting was before the current window — which drags the left edge backwards and invents substrings that were never valid. The shortest input that exposes it is "abba":

def longest_unique_unguarded(text: str) -> int:
    """The same code with the `previous >= start` guard removed."""
    last_seen: dict[str, int] = {}
    start = 0
    best = 0
    for index, character in enumerate(text):
        if character in last_seen:
            start = last_seen[character] + 1
        last_seen[character] = index
        best = max(best, index - start + 1)
    return best


def longest_unique_brute(text: str) -> int:
    """Check every substring. Used only to confirm the fast version."""
    best = 0
    for start in range(len(text)):
        for end in range(start + 1, len(text) + 1):
            piece = text[start:end]
            if len(set(piece)) == len(piece):
                best = max(best, end - start)
    return best


for sample in ["abba", "abcabcbb", "pwwkew", "bbbbb", ""]:
    print(f"{sample!r:<10} window={longest_unique(sample)} "
          f"unguarded={longest_unique_unguarded(sample)} "
          f"brute={longest_unique_brute(sample)}")

'abba'     window=2 unguarded=3 brute=2
'abcabcbb' window=3 unguarded=3 brute=3
'pwwkew'   window=3 unguarded=3 brute=3
'bbbbb'    window=1 unguarded=1 brute=1
''         window=0 unguarded=0 brute=0

Trace "abba" to see it. At i=2 the second b arrives, so start moves to 2 and the window is "b". At i=3 the a arrives; it was last seen at index 0, which is behind the left edge, so it is not in the window at all and nothing should move. The unguarded version sets start back to 1, claims the window "bba" has length 3, and reports 3 for a string whose answer is 2. Every other test agrees, which is exactly why this bug survives code review.

Edge cases fall out of the arithmetic. An empty string never enters the loop and returns 0. A string of identical characters shrinks on every step and returns 1. In the fixed-size version, k equal to len(values) gives an empty range(k, len(values)), so the initial sum is the answer.

Two more windows worth memorising

Minimum window substring

Given a text and a target, find the shortest substring of the text containing every character of the target, counting repeats. This is the problem that teaches the need/have counter: a dict of required counts, a dict of current counts, and a single integer matched recording how many distinct required characters are currently present at full strength. The window is valid exactly when matched equals the number of distinct characters in the target — one integer comparison, not a dict comparison.

There is one structural difference from every window so far. For a longest answer you shrink until the window is valid, then record. For a shortest answer you record while the window is still valid, then shrink to try for better. Same loop, opposite placement of the recording line.

def min_window(text: str, target: str) -> str:
    """Shortest substring of text containing every character of target.

    need holds the required count of each character. matched counts how many
    distinct required characters currently sit in the window at full strength,
    so the window is valid exactly when matched equals len(need).
    """
    if not target or len(target) > len(text):
        return ""

    need: dict[str, int] = {}
    for character in target:
        need[character] = need.get(character, 0) + 1

    have: dict[str, int] = {}
    matched = 0
    left = 0
    best_start, best_length = 0, len(text) + 1

    for right, character in enumerate(text):
        if character in need:
            have[character] = have.get(character, 0) + 1
            if have[character] == need[character]:
                matched += 1

        # Shrink while the window is still valid: a valid window can only get
        # shorter by losing characters from the left, so record before evicting.
        while matched == len(need):
            if right - left + 1 < best_length:
                best_start, best_length = left, right - left + 1
            leaving = text[left]
            if leaving in need:
                have[leaving] -= 1
                if have[leaving] < need[leaving]:
                    matched -= 1
            left += 1

    return "" if best_length > len(text) else text[best_start:best_start + best_length]


print(min_window("ADOBECODEBANC", "ABC"))
print(repr(min_window("a", "aa")), repr(min_window("a", "a")), repr(min_window("ab", "b")))

BANC
'' 'a' 'b'

On "ADOBECODEBANC" the window becomes valid three separate times. The first valid window is "ADOBEC" at indices 0 to 5, length 6. It is recorded, then the leading A is evicted and the window is invalid again until the A at index 10 arrives. That second stretch shrinks all the way down to "CODEBA", length 6 — a tie, so nothing is recorded — before its C is evicted. The final C at index 12 makes it valid a last time, and shrinking from "ODEBANC" gives "EBANC" at length 5 and then "BANC" at length 4.

Note the have[leaving] < need[leaving] test on eviction. Dropping a B when the window holds two of them does not invalidate anything, so matched must only fall when the count actually drops below what is required.

Longest run with at most k distinct values

The same skeleton with a plain counter. The rule is len(counts) being at most k, and the one detail that catches people is that a key whose count reaches zero must be deleted, because len() counts keys, not non-zero keys.

def longest_at_most_k_distinct(values: list[int], k: int) -> int:
    """Length of the longest run containing at most k distinct values."""
    if k <= 0:
        return 0

    counts: dict[int, int] = {}
    left = 0
    best = 0

    for right, value in enumerate(values):
        counts[value] = counts.get(value, 0) + 1
        while len(counts) > k:
            leaving = values[left]
            counts[leaving] -= 1
            if counts[leaving] == 0:
                # A key left at zero would still be counted by len(counts).
                del counts[leaving]
            left += 1
        best = max(best, right - left + 1)

    return best


def window_travel(values: list[int], k: int) -> tuple[int, int, int]:
    """Return (best length, steps right took, steps left took)."""
    counts: dict[int, int] = {}
    left = 0
    best = 0
    right_steps = left_steps = 0

    for right, value in enumerate(values):
        right_steps += 1
        counts[value] = counts.get(value, 0) + 1
        while len(counts) > k:
            leaving = values[left]
            counts[leaving] -= 1
            if counts[leaving] == 0:
                del counts[leaving]
            left += 1
            left_steps += 1
        best = max(best, right - left + 1)

    return best, right_steps, left_steps


for sample, width in [([1, 2, 1, 3, 4], 2),
                      ([(index * index) % 13 for index in range(100_000)], 5)]:
    best, right_steps, left_steps = window_travel(sample, width)
    print(f"n={len(sample):>7,}  k={width}  best={best:>2}  right moved "
          f"{right_steps:>7,}  left moved {left_steps:>6,}  total "
          f"{right_steps + left_steps:>7,}  (2n = {2 * len(sample):,})")

print(longest_at_most_k_distinct([1, 2, 1, 3, 4], 2),
      longest_at_most_k_distinct([1, 1, 1], 1),
      longest_at_most_k_distinct([], 3))

n=      5  k=2  best= 3  right moved       5  left moved      3  total       8  (2n = 10)
n=100,000  k=5  best=10  right moved 100,000  left moved 99,992  total 199,992  (2n = 200,000)
3 3 0

collections.Counter will do the counting dict for you, but it does not delete keys that hit zero, so len() on a Counter after a decrement is not the distinct count. Either del the key yourself, as above, or track the distinct count in a separate integer.

This function is also the building block for the harder question "how many runs contain exactly k distinct values". Counting runs with exactly k distinct is counting runs with at most k, minus counting runs with at most k − 1 — the same window run twice, with a running total of window lengths instead of a maximum.

The precondition, and where sliding window breaks

Here is the failure everyone should see once. The problem: find the longest run whose sum is at most a given limit. The window version is four lines and it is correct on non-negative input.

def longest_at_most_sum(values: list[int], limit: int) -> int:
    """Longest run whose sum is at most limit. Correct only for non-negative values."""
    left = 0
    total = 0
    best = 0

    for right, value in enumerate(values):
        total += value
        while left <= right and total > limit:
            total -= values[left]
            left += 1
        best = max(best, right - left + 1)

    return best


def longest_at_most_sum_brute(values: list[int], limit: int) -> int:
    """Check every run. Slow, but right on any input."""
    best = 0
    for start in range(len(values)):
        running = 0
        for end in range(start, len(values)):
            running += values[end]
            if running <= limit:
                best = max(best, end - start + 1)
    return best


non_negative = [3, 1, 2, 4, 1]
with_a_negative = [2, 5, -5, 1]
print(longest_at_most_sum(non_negative, 6), longest_at_most_sum_brute(non_negative, 6))
print(longest_at_most_sum(with_a_negative, 3), longest_at_most_sum_brute(with_a_negative, 3))

agree = all(
    longest_at_most_sum(sample, 10) == longest_at_most_sum_brute(sample, 10)
    for sample in (
        [(index * 7 + offset * 5) % 9 for index in range(30)] for offset in range(200)
    )
)
print("non-negative samples agree:", agree)

3 3
2 4
non-negative samples agree: True

Two hundred non-negative inputs agree. [2, 5, -5, 1] with a limit of 3 does not: the window says 2, the truth is 4.

Trace it. At right = 1 the window is [2, 5] with a sum of 7, over the limit, so the shrink loop evicts the 2 (sum 5, still over) and then the 5 (sum 0), leaving left at 2 with an empty window. From there the best it can ever report is a run starting at index 2, and the answer — the whole array, which sums to exactly 3 — starts at index 0. That start was thrown away and left cannot go back.

The monotonicity precondition is what failed. With a negative in the array, removing a value from the left can increase the sum, so "invalid" no longer implies "invalid for every earlier left edge". There is no threshold to shrink towards, and the entire justification for a one-directional left pointer collapses. No exception, no warning, just a smaller number than the right one.

The repair is not a patched window, it is a different algorithm. Write prefix[i] for the sum of the first i values; then the sum of a run from i to j is prefix[j + 1] - prefix[i], and the question becomes a search over prefix values rather than a walk over window edges. Prefix sums covers that shift. The related question "is there a run summing to exactly S", which sliding windows also cannot answer once negatives appear, is solved in O(n) by storing every prefix sum in a hash table and looking up prefix - S at each step.

The same caution applies whenever the rule is not monotone under shrinking. "Longest run whose product is at most P" breaks the moment a zero or a value below 1 appears, for exactly the same reason.

Complexity

Fixed window: O(n) time, O(1) space. The first window costs k reads. Each of the remaining n - k steps costs exactly two reads and one comparison. Total: k + 2(n - k), which is at most 2n, so O(n). The naive version is (n - k + 1) * k, which is O(n · k) — and when k is proportional to n, that is O(n²).

Variable window: O(n) time, amortised. This is the argument to be able to give, because the code contains a while loop inside a for loop and that shape is quadratic in most other contexts.

Count the two edges separately instead of counting loop iterations:

right advances exactly once per outer iteration, so it takes exactly n steps over the whole run.
left starts at 0, only ever increases, and never passes n. So left takes at most n steps in total, across the entire run — not per outer iteration, in total.

Each step of either edge does a constant amount of work: one dict update, one comparison, one arithmetic operation. So the total work is bounded by n + n = 2n constant-cost steps, which is O(n).

One caveat worth naming rather than glossing over. A dict update is O(1) on average, not guaranteed — a set of keys that all hash to the same bucket drags a single lookup out to O(n). So the counting windows are O(n) expected, not O(n) come what may. The sum-based windows hold no dict and carry no such caveat.

The reason the nested loop does not multiply is that the two loops share one budget. A nested loop is quadratic when the inner count is independent of the outer loop's progress — an inner loop that runs n times for each of n outer iterations. Here every inner iteration permanently consumes one of the n moves left will ever make. A single outer iteration can trigger a long shrink, but that shrink is stolen from the budget of every later iteration. The measurement above says it plainly: on 100,000 elements, right moved 100,000 times and left moved 99,992 times, for 199,992 total moves against a ceiling of 200,000.

Space. O(1) for the sum-based windows: two indices and a running total. The counting windows hold a dict, bounded by the number of distinct values that can be inside a valid window — at most k + 1 keys for the at-most-k-distinct version, at most one key per distinct character for longest_unique (128 for plain ASCII text, so effectively constant), and at most one key per distinct character of the target for min_window.

Against the brute force, on a list of n items:

Items	Runs to check	Window pointer moves
100	5,050	200
1,000	500,500	2,000
10,000	50,005,000	20,000
100,000	5,000,050,000	200,000

The left column is n(n + 1) / 2, the number of contiguous runs in a list of n items. Ten times the data costs a hundred times the work there and ten times the work on the right.

When to use it, and when not to

Use it when three things are true: the answer is a contiguous run, the property you are testing can be updated in constant time when one element joins and one leaves, and shrinking from the left can restore validity. Sums, counts, character frequencies and "at most k of something" all satisfy the second condition. Maximums do not, which is the one common exception worth knowing — when the largest value in the window is the one leaving, you cannot recover the new largest in O(1) without extra structure. That problem, sliding window maximum, is solved with a monotonic collections.deque holding indices in decreasing value order, and it is still O(n) overall because each index is pushed and popped once.

Do not use it when the answer is not contiguous. "Longest increasing subsequence" and "longest common subsequence" allow gaps, and no window can express a gap. Those are dynamic programming problems.

Do not use it when the answer is a pair of positions rather than a range. That is the two pointers technique, a close relative where the two indices usually move towards each other rather than both rightwards, and where what sits between them is irrelevant.

Do not use it when the rule is not monotone under shrinking, as shown above. Sums with negatives, products with zeros or fractions, and any rule where removing an element can make things worse are all disqualified.

Do not use it for arbitrary range queries either. If the question is "the sum of indices 300 to 700" asked a thousand times over a fixed array, a window has nothing to slide along; build a prefix-sum array once and answer each query in O(1).

Where it shows up in the real world

TCP flow control. The receiver advertises a window size, and the sender may have at most that many unacknowledged bytes in flight. As acknowledgements arrive the window's left edge advances and new bytes become sendable on the right. It is the same two-edge structure, running on every connection your machine has open.

DEFLATE, the algorithm behind gzip, zlib and PNG. Its LZ77 stage encodes repeated data as back-references into a sliding window of the previous 32 KB of output. Data older than 32 KB has fallen out of the window and can no longer be referenced, which is precisely the trade that keeps the compressor's memory constant regardless of file size.

Rolling hashes. Rabin-Karp slides a fixed-width window over the text and updates the hash by removing the leaving character's contribution and adding the entering one's — the add-one-remove-one update on a hash instead of a sum. Content-defined chunking in backup and deduplication systems uses the same rolling window to decide where to cut a file.

Stream processing. Apache Flink and Kafka Streams both expose windowing as a first-class operator: tumbling windows that never overlap, and sliding or hopping windows that advance by less than their width. The aggregation kept per window is maintained incrementally as records enter and expire.

API rate limiting. The sliding-window-counter algorithm approximates a rolling limit by blending the current fixed window's count with a weighted share of the previous window's, which gives per-client memory of two integers instead of one timestamp per request. Cloudflare has published this approach for its rate limiter.

Signal processing. The short-time Fourier transform behind every spectrogram slides a fixed-width window along a waveform, usually with overlap, and transforms each position independently.

Common mistakes

Recomputing the property inside the loop. Writing sum(values[left:right + 1]) or len(set(text[left:right + 1])) inside the loop is the single most common way to write a sliding window that is secretly O(n²). The whole technique is the incremental update; if you recompute, you have kept the code and thrown away the algorithm.

Using if instead of while to shrink. One eviction is often not enough. On [1, 2, 1, 3] with k of 2, when the 3 arrives the window holds three distinct values; evicting index 0 removes a 1 but another 1 remains at index 2, so the count is still 3 and a second eviction is required.

Letting the left edge move backwards. The "abba" bug above. Any expression that assigns left a value without taking a maximum against its current value should be treated as suspect until you have proved the new value cannot be smaller.

Leaving zero-count keys in the dict. counts[value] -= 1 without the del when it reaches zero means len(counts) keeps counting a value that is no longer in the window, and the window never shrinks enough.

Recording the answer in the wrong place. For a longest-window problem, record after the shrink loop, when the window is valid. For a shortest-window problem, record inside the shrink loop, before evicting. Swapping them gives answers that are plausible and wrong.

Mixing up the two length conventions. With both ends inclusive the length is right - left + 1. With right exclusive it is right - left. Pick one and use it in every line, including the initial best.

Assuming a window applies because the problem mentions subarrays. Check the monotonicity precondition first. If shrinking cannot repair an invalid window, no amount of careful coding will save it.

Practice

Given a list of numbers and a width k, return the largest average of any k consecutive values.
Given a string, find the length of the longest substring containing at most two distinct characters.
Given a list of non-negative integers and a target, find the length of the shortest run whose sum is at least the target, returning 0 if none exists.
Given two strings, report every index in the first where an anagram of the second begins, using a fixed-width window of character counts.
Given a list and a width k, return the maximum of every window of that width, keeping the whole thing O(n) with a collections.deque of indices held in decreasing value order.

Summary

A sliding window replaces recomputation with repair. Fixed-size windows patch a running total with one addition and one subtraction; variable-size windows grow on the right and shrink on the left until the rule holds again. The linear bound comes from counting edge movements rather than loop iterations — each of the two edges makes at most n forward moves over the whole run, so the nested while costs 2n steps in total and not n² — and the correctness comes from a precondition that is easy to state and easy to forget: shrinking from the left must be able to restore validity. Verify that condition before you write the loop, because when it fails the code still runs.


Difficulty	Medium
Time (fixed window)	O(n) — k reads for the first window, then 2 per step
Time (variable window)	O(n) amortised — each edge makes at most n forward moves, so at most 2n steps total
Time (brute force)	O(n²) — a list of n items has n(n + 1) / 2 contiguous runs
Space	O(1) for sum windows; O(d) for counting windows, d = distinct values a valid window may hold
Precondition	Shrinking from the left must restore validity; fails on sums with negative values
Data structure	List or string with O(1) indexing, plus a dict for the counting variants
Use it when	The answer is a contiguous run and the property updates in O(1) as elements join and leave
Avoid it when	Values can be negative for a sum rule, or the answer allows gaps
Real-world use	TCP flow control, DEFLATE's 32 KB back-reference window, rolling hashes, Flink and Kafka Streams windowing
Python equivalent	`collections.deque(maxlen=k)` for a fixed-width buffer, `collections.Counter` for the counts

Keep reading

The Two Pointers Technique — the sibling pattern, for when the answer is a pair of positions rather than a run.
Prefix Sums — what to reach for when negatives break the window, and how to answer range queries in O(1).
Hash Tables — why the counting dicts above are O(1) per update on average, and the prefix-sum lookup trick for negatives.
Rabin-Karp — a rolling hash is the add-one-remove-one update applied to a hash function.
Big O Notation — the amortised counting argument used here, from first principles.

The Two Pointers Technique in Python: Turning O(n ) Into O(n)

Bimal Kshetri — Sun, 16 Aug 2026 16:38:51 +0000

A nested loop over every pair in a list of 1,000 items inspects 499,500 pairs. Two pointers walking that same list inspect at most 999. That gap is not a tuning detail — it is the difference between code you can run on a million items and code you cannot, and it usually costs about four lines.

Two pointers is a pattern rather than a single algorithm. You keep two indices into one sequence and move them under a rule that guarantees each index only ever travels in one direction. Because neither index ever backtracks, the total number of moves is bounded by the length of the sequence, and a problem that looked like it needed every pair collapses into a single sweep.

The catch is that the rule has to be provably safe. Moving a pointer throws away a large set of candidate answers unseen, and that is only legal if none of them could have been the answer. Most of this post is those proofs, because they are the hard part. Once you have the proof, the code writes itself.

The idea

Take a list. Put one index at position 0 and another at position len(values) - 1. Look at the pair they point at. Based on what you see, move one of them inwards. Repeat until they meet.

That is the first shape, converging pointers. The second puts both indices near the front and moves both rightwards under different conditions: one reads every element, the other advances only when something interesting happens. That is same-direction pointers, and the fast-and-slow pair used to find cycles in a linked list is the same shape with a fixed speed ratio.

Both shapes share the property that makes the pattern fast: each pointer moves one way only. In the converging shape left starts at 0 and only increases, right starts at len(values) - 1 and only decreases, and the loop stops the instant they meet. Every turn moves at least one of them, so the body runs at most n − 1 times. If the body costs a constant amount, the algorithm is O(n).

That counting argument is easy. The safety argument is not: why is it correct to move on rather than back up and try the pairs you skipped? The answer differs per problem, and it always comes from some ordering property of the data — which is why most converging-pointer problems need sorted input.

Watching it work

Given a list sorted in ascending order, find two positions whose values add up to a target. Take [1, 3, 4, 6, 8, 11] and a target of 10, with left at index 0 and right at index 5.

1 + 11 = 12. Too big. Move right in to index 4.
1 + 8 = 9. Too small. Move left in to index 1.
3 + 8 = 11. Too big. Move right in to index 3.
3 + 6 = 9. Too small. Move left in to index 2.
4 + 6 = 10. Found it, at indices 2 and 3.

Five steps. There are 15 distinct pairs in a six-element list, and a nested loop would have stumbled onto this answer on its tenth. But the count is not the interesting bit — the interesting bit is why the skipped pairs never needed checking at all.

Why discarding a whole row is safe

Draw every pair as a table. Rows are the left value, columns are the right value, each cell holds their sum. Only cells where the left index is smaller than the right index are real pairs; the rest are blank.

This table is sorted along both axes, and that follows directly from the list being sorted. Move right along a row: the left value is fixed while the right value grows, so the sums grow. Move down a column: the right value is fixed while the left value grows, so the sums grow again.

Now look at where the algorithm starts — the top-right corner. That cell is simultaneously the largest sum in its row and the smallest sum in its column. Everything follows from that:

If the corner is too big, it is the smallest entry in its column, so every other sum in that column is bigger still. Nothing there can be the target. Delete the whole column, which is exactly what right -= 1 does.
If the corner is too small, it is the largest entry in its row, so every other sum in that row is smaller still. Delete the whole row, which is exactly what left += 1 does.

Each step retires a full row or column without looking inside it. There are only n rows and n columns to retire, so the walk ends within 2n steps even though the table holds n(n − 1) / 2 cells. The path is a staircase from the corner down to the answer — the same walk used to search a matrix whose rows and columns are both sorted.

Keep both halves of that in your head. "Two pointers is O(n)" is a fact about the loop; "each move retires a row or a column" is why the loop is allowed to skip that much.

The code

Start with the version that needs no cleverness, so there is something to compare against.

# Postponed annotation evaluation, so `X | None` works on Python 3.9 too.
from __future__ import annotations


def pair_sum_brute_force(values: list[int], target: int) -> tuple[int, int] | None:
    """Return the indices of the first pair that adds to target, or None.

    Works on any list, sorted or not, because it checks every pair.
    """
    size = len(values)
    for left in range(size):
        for right in range(left + 1, size):
            if values[left] + values[right] == target:
                return left, right
    return None


numbers = [1, 3, 4, 6, 8, 11]
print(pair_sum_brute_force(numbers, 10))
print(pair_sum_brute_force(numbers, 100))

(2, 3)
None

Now the two-pointer version, plus a traced copy that prints the walkthrough above so you can check the hand trace against the machine.

def pair_sum(values: list[int], target: int) -> tuple[int, int] | None:
    """Find two indices in an ascending-sorted list whose values add to target.

    Returns the pair of indices, or None if no such pair exists. Runs in O(n)
    time and O(1) space. The list must already be sorted; this does not check.
    """
    left, right = 0, len(values) - 1

    while left < right:
        total = values[left] + values[right]
        if total == target:
            return left, right
        if total < target:
            # Even the largest partner still in play was too small, so no
            # remaining pair that uses values[left] can reach the target.
            left += 1
        else:
            # Even the smallest partner still in play was too large.
            right -= 1

    return None


def pair_sum_traced(values: list[int], target: int) -> tuple[int, int] | None:
    """The same algorithm, printing every step it takes."""
    left, right = 0, len(values) - 1
    while left < right:
        total = values[left] + values[right]
        if total == target:
            verdict = "found"
        elif total < target:
            verdict = "too small, move left in"
        else:
            verdict = "too big, move right in"
        print(f"left={left} right={right}  {values[left]:>2} + {values[right]:>2}"
              f" = {total:>2}  {verdict}")
        if total == target:
            return left, right
        if total < target:
            left += 1
        else:
            right -= 1
    return None


print(pair_sum_traced(numbers, 10))
print(pair_sum(numbers, 100))

left=0 right=5   1 + 11 = 12  too big, move right in
left=0 right=4   1 +  8 =  9  too small, move left in
left=1 right=4   3 +  8 = 11  too big, move right in
left=1 right=3   3 +  6 =  9  too small, move left in
left=2 right=3   4 +  6 = 10  found
(2, 3)
None

The gap between the two grows fast. Here both versions hunt for a target that does not exist in a sorted list of 1,000 even numbers, so neither can exit early, and the counters record how many pairs each one actually inspects.

def count_brute_force(values: list[int], target: int) -> int:
    """How many pairs the nested loops inspect."""
    looks = 0
    for left in range(len(values)):
        for right in range(left + 1, len(values)):
            looks += 1
            if values[left] + values[right] == target:
                return looks
    return looks


def count_two_pointers(values: list[int], target: int) -> int:
    """How many pairs the converging pointers inspect."""
    looks, left, right = 0, 0, len(values) - 1
    while left < right:
        looks += 1
        total = values[left] + values[right]
        if total == target:
            return looks
        if total < target:
            left += 1
        else:
            right -= 1
    return looks


evens = [2 * index for index in range(1000)]
print(count_brute_force(evens, 1), count_two_pointers(evens, 1))

499500 999

How the code maps to the idea

while left < right is strict, not left <= right. With <= the loop eventually reaches a state where both indices point at the same element, and pair_sum would report that element paired with itself — a real bug when the target is exactly twice some value in the list. Strictness also gives the termination proof: the gap shrinks by at least one every turn and the loop stops at zero.

The three-way branch is the staircase walk. Equality is the answer, total < target deletes a row, total > target deletes a column. Exactly one pointer moves per branch, which keeps both claims true: the loop always makes progress, and it never discards a row and a column in the same step.

Nothing is stored. Two integers and a sum, whatever the size of the list, so extra space is O(1). That is often the real reason to prefer it over a lookup table.

Edge cases fall out of the arithmetic. An empty list gives left = 0 and right = -1, so the loop never runs. A one-element list gives left = 0 and right = 0, likewise. Neither needs a guard.

The sorted precondition is not checked, deliberately: checking costs O(n) and defeats the point of a function you call inside a loop. On unsorted input it does not raise, it silently returns the wrong answer — which is worse, so document the precondition and mean it.

The other variants worth knowing

The pair sum is the cleanest example, but the shape shows up in a dozen disguises. These are the ones worth being able to write from memory.

Palindromes

The same converging walk, without the arithmetic. Compare the ends, step inwards, stop the moment they disagree. There is nothing to prove beyond the definition: a string is a palindrome exactly when every character matches its mirror, and the loop checks each mirror pair once.

def is_palindrome(text: str) -> bool:
    """True if text reads the same both ways, ignoring case and punctuation."""
    cleaned = [character.lower() for character in text if character.isalnum()]
    left, right = 0, len(cleaned) - 1

    while left < right:
        if cleaned[left] != cleaned[right]:
            return False
        left += 1
        right -= 1

    return True


for candidate in ["A man, a plan, a canal: Panama", "race a car", "ab", "a", ""]:
    print(f"{candidate!r:34} {is_palindrome(candidate)}")

'A man, a plan, a canal: Panama'   True
'race a car'                       False
'ab'                               False
'a'                                True
''                                 True

Writing this as cleaned == cleaned[::-1] is shorter and, in CPython, faster — the reversal runs in C. The two-pointer version wins on space, because it never builds the second copy, and it bails out on the first mismatch.

Reversing in place

Swap the ends, then close in. The loop runs len(values) // 2 times, exactly the minimum number of swaps needed.

def reverse_in_place(values: list) -> None:
    """Reverse a list by swapping the two ends and closing inwards."""
    left, right = 0, len(values) - 1
    while left < right:
        values[left], values[right] = values[right], values[left]
        left += 1
        right -= 1


letters = list("abcdef")
reverse_in_place(letters)
print("".join(letters))

digits = [1, 2, 3, 4, 5]
reverse_in_place(digits)
print(digits)

fedcba
[5, 4, 3, 2, 1]

In real Python you would write values.reverse(), which is this loop implemented in C. Write it once by hand anyway; it is the smallest version of the pattern there is.

Removing duplicates from a sorted list, in place

This is the same-direction shape. One pointer, read, visits every index. Another, last_kept, marks the end of the answer being built at the front of the same list. Because the list is sorted, equal values are always adjacent, so one comparison against the last value kept decides whether the current one is new.

def remove_duplicates(values: list[int]) -> int:
    """Compact a sorted list in place so that values[:kept] are its distinct
    values, and return kept.

    Everything from index kept onwards is left as whatever happened to be
    there and must be ignored. O(n) time, O(1) extra space.
    """
    if not values:
        return 0

    last_kept = 0
    for read in range(1, len(values)):
        # values is sorted, so a repeat can only ever sit next to the value it
        # repeats. One comparison against the last value kept is enough.
        if values[read] != values[last_kept]:
            last_kept += 1
            values[last_kept] = values[read]

    return last_kept + 1


sample = [1, 1, 1, 2, 3, 3]
kept = remove_duplicates(sample)
print(kept, sample[:kept], sample)

longer = [0, 0, 1, 1, 1, 2, 2, 3, 3, 4]
kept = remove_duplicates(longer)
print(kept, longer[:kept])

empty: list[int] = []
print(remove_duplicates(empty), empty)

3 [1, 2, 3] [1, 2, 3, 2, 3, 3]
5 [0, 1, 2, 3, 4]
0 []

Notice the printed tail: sample ends up as [1, 2, 3, 2, 3, 3]. The function does not shorten the list, it only guarantees the first kept entries are correct, and forgetting that is the most common bug in this variant. The write pointer can never overtake the read pointer — last_kept starts behind and advances at most once per turn — so nothing is overwritten before it has been read.

Fast and slow

Give the two same-direction pointers a fixed speed ratio instead of a condition and you get a different tool. Advance slow one node per turn and fast two, and when fast falls off the end, slow sits at the midpoint. Run the same pair inside a linked list that loops back on itself and they must collide, because fast closes the gap by exactly one node per turn, so inside a loop of length L they meet within L turns. That is Floyd's tortoise and hare.

class Node:
    """A singly linked list node, here only to demonstrate fast and slow."""

    def __init__(self, value: int) -> None:
        self.value = value
        self.next: Node | None = None


def build_list(values: list[int]) -> Node | None:
    head: Node | None = None
    for value in reversed(values):
        node = Node(value)
        node.next = head
        head = node
    return head


def middle_value(head: Node | None) -> int | None:
    """The value halfway along a linked list, found in a single pass.

    slow advances one node per turn and fast advances two, so fast reaches
    the end after slow has covered exactly half the nodes.
    """
    slow, fast = head, head
    while fast is not None and fast.next is not None:
        slow = slow.next
        fast = fast.next.next
    return slow.value if slow is not None else None


def has_cycle(head: Node | None) -> bool:
    """Floyd's tortoise and hare.

    fast closes the gap on slow by exactly one node per turn, so once both are
    inside a loop of length L they must collide within L turns.
    """
    slow, fast = head, head
    while fast is not None and fast.next is not None:
        slow = slow.next
        fast = fast.next.next
        if slow is fast:
            return True
    return False


odd_list = build_list([1, 2, 3, 4, 5])
even_list = build_list([1, 2, 3, 4])
print(middle_value(odd_list), middle_value(even_list), middle_value(None))
print(has_cycle(odd_list))

looped = build_list([1, 2, 3, 4, 5])
tail = looped
while tail.next is not None:
    tail = tail.next
tail.next = looped.next  # the last node now points back at the second
print(has_cycle(looped))

3 3 None
False
True

A list of even length has two middles and this version returns the second, which is why [1, 2, 3, 4] prints 3. Cycle detection is where the shape really earns its keep: the obvious alternative, a set holding every node already visited, costs O(n) memory. Two pointers cost two.

Container with most water

Treat each entry in a list of heights as a vertical line and pick the two lines that hold the most water between them. The area of a pair is the distance between them times the height of the shorter one, because water spills over the shorter side.

This is the variant worth studying, because the input is not sorted and the proof still works. Start at the two ends, the widest pair possible. Whichever line is shorter caps the area of every pair it belongs to, and every pair it has left is narrower than the one just measured. Its best remaining area is therefore strictly less than the area already recorded, so it can be discarded — safely, permanently, without looking.

def max_water(heights: list[int]) -> int:
    """The most water two of the vertical lines can hold between them.

    The area of a pair is the distance between the lines times the height of
    the shorter one, because water spills over the shorter side.
    """
    left, right = 0, len(heights) - 1
    best = 0

    while left < right:
        area = (right - left) * min(heights[left], heights[right])
        best = max(best, area)
        # The shorter line caps every rectangle it appears in, and every pair
        # it has left is narrower than this one, so it can never do better.
        if heights[left] < heights[right]:
            left += 1
        else:
            right -= 1

    return best


def max_water_brute_force(heights: list[int]) -> int:
    """Check all n(n - 1) / 2 pairs. Used here only to verify the fast one."""
    best = 0
    for left in range(len(heights)):
        for right in range(left + 1, len(heights)):
            best = max(best, (right - left) * min(heights[left], heights[right]))
    return best


print(max_water([1, 8, 6, 2, 5, 4, 8, 3, 7]))

agree = all(
    max_water(profile) == max_water_brute_force(profile)
    for profile in (
        [(index * index * 7 + offset * 31) % 29 + 1 for index in range(40)]
        for offset in range(300)
    )
)
print("matches the brute force on all 300 profiles:", agree)

49
matches the brute force on all 300 profiles: True

Checking against the quadratic version is not a proof, but it is the right habit. When a linear algorithm claims to match an exhaustive one, run both on a few hundred inputs before believing yourself.

3-sum: fix one, two-pointer the rest

Find every distinct triple that adds to zero. The exhaustive version is three nested loops, O(n³). The fix is to stop treating it as a three-dimensional problem: sort the list, fix the first value, and what remains is "find two values that add to minus the fixed one" — the pair sum you already have.

Sorting does double duty. It makes the inner sweep legal, and it puts equal values next to each other so duplicate triples are cheap to skip.

def three_sum(values: list[int]) -> list[tuple[int, int, int]]:
    """Every distinct triple of values that adds to zero, smallest first.

    Sorting is what makes the inner two-pointer sweep legal, and it is also
    what makes duplicates cheap to skip: equal values end up side by side.
    """
    numbers_sorted = sorted(values)
    size = len(numbers_sorted)
    triples: list[tuple[int, int, int]] = []

    for first in range(size - 2):
        if numbers_sorted[first] > 0:
            break  # three numbers all above zero cannot add to zero
        if first > 0 and numbers_sorted[first] == numbers_sorted[first - 1]:
            continue  # this value has already been the fixed element

        left, right = first + 1, size - 1
        while left < right:
            total = numbers_sorted[first] + numbers_sorted[left] + numbers_sorted[right]
            if total < 0:
                left += 1
            elif total > 0:
                right -= 1
            else:
                triples.append(
                    (numbers_sorted[first], numbers_sorted[left], numbers_sorted[right])
                )
                left += 1
                right -= 1
                while left < right and numbers_sorted[left] == numbers_sorted[left - 1]:
                    left += 1
                while left < right and numbers_sorted[right] == numbers_sorted[right + 1]:
                    right -= 1

    return triples


print(three_sum([-1, 0, 1, 2, -1, -4]))
print(three_sum([0, 0, 0, 0]))
print(three_sum([1, 2, 3]))
print(three_sum([-2, 0, 1, 1, 2]))

[(-1, -1, 2), (-1, 0, 1)]
[(0, 0, 0)]
[]
[(-2, 0, 2), (-2, 1, 1)]

The break on the first positive value is a small win worth understanding: once the smallest of the three is above zero, all three are, so no triple from here on can reach zero. On a list of mostly positive numbers that ends the outer loop almost immediately.

Complexity

The converging sweep is O(n). Each turn does a fixed amount of work — one addition, one comparison, one increment — and moves exactly one pointer. The two pointers share a single budget: the gap between them starts at n − 1 and shrinks by one per turn, so the loop runs at most n − 1 times. There is no hidden cost inside the body. That n − 1 is the worst case, and every miss pays it in full, because the loop only stops when the pointers meet. The best case is O(1) — the first pair tried, smallest plus largest, is already the target. A hit stops the moment it is found, so on average you pay about half the sweep, which is still O(n).

Space is O(1) for the sweep itself. Two indices and a running total, regardless of input size. Two of the variants above break that, and the break is in the code rather than in the idea. is_palindrome builds a cleaned list before it starts, so it holds O(n); skip the non-alphanumerics inline instead and it is O(1) again. three_sum calls sorted(), which returns a new list, so it holds O(n) as well — and then O(n) more for the triples it collects.

Sorting, when you need it, dominates. If the input arrives unsorted, the honest bound for the pair sum is O(n log n) for the sort plus O(n) for the sweep, which is O(n log n). The log comes from the sort: merge sort halves the problem repeatedly, giving log₂ n levels with all n items touched at each level. Timsort, the algorithm behind Python's sorted(), gets to the same bound a different way — it merges runs that are already in order rather than halving — but its worst case is also O(n log n). Sorting costs memory too, since sorted() builds a whole new list. Do not quote O(n) time or O(1) space for a function that sorts first.

3-sum is O(n²). The outer loop fixes each of n values in turn, and each one costs an O(n) sweep: n × n. The sorted() call adds O(n log n), which n² swallows. That is a big improvement on the O(n³) triple loop, but it is still quadratic — 10,000 inputs means roughly 100 million inner steps.

Fast and slow is O(n) time and O(1) space. fast advances two nodes per turn and falls off the end within n / 2 turns; with a cycle it collides within one extra lap.

Here is the comparison that matters, on the pair sum with the input already sorted:

Items	Pairs a nested loop checks	Steps two pointers take
100	4,950	99
1,000	499,500	999
10,000	49,995,000	9,999
100,000	4,999,950,000	99,999

Ten times the data costs a hundred times the work on the left and ten times on the right. At 100,000 items that is exactly 50,000 times fewer steps.

When to use it, and when not to

Use it when the input is sorted (or you can sort it cheaply), the answer is a pair of positions rather than a range, and extra space has to stay constant. Sorted arrays, strings compared from both ends, in-place compaction and merging two sorted sequences are all natural fits.

Do not use it on unsorted input just because the shape is familiar. For "do two values add to this target" on unsorted data, one pass with a set is better: O(n) expected time and no sort at all. It costs O(n) memory, and that is the trade. Two pointers only win when the data is already sorted or memory is genuinely tight. Hash tables explain why that lookup is O(1) on average — average is the whole guarantee, and keys chosen to collide can drag a single lookup back to O(n).

Do not use it when the answer is a contiguous run. That is the sliding window, which is a two-pointer method too — both ends move rightwards — but the question is different. A window cares about everything between its ends: the sum of the run, the number of distinct characters in it, whether it is still valid. Converging pointers care about the pair itself and mostly ignore what sits between them. If the answer is a subarray or substring, reach for a window; if it is two positions, reach for two pointers.

Do not use it when the property you rely on is not monotonic. The pattern rests entirely on "moving this pointer can only change things in one predictable direction". If moving left could send the sum either way, the safety proof collapses and the algorithm returns wrong answers without complaining.

Where it shows up in the real world

Merging sorted sequences. The merge step of merge sort is two pointers, one into each of two sorted halves, each advancing when its value is taken. The same walk is a sort-merge join in relational databases: PostgreSQL's Merge Join advances a cursor over each of two sorted inputs, taking whichever key is smaller. It is also how LSM-tree storage engines such as RocksDB and LevelDB compact several sorted files into one.

CPython's string stripping. str.strip() walks an index forward from the start while the character is whitespace, then walks a second index backward from the end under the same test, and slices between them. It is a converging pair with a very simple rule.

Mark-compact garbage collectors. Sliding compaction walks the heap with a read pointer and a write pointer, copying every live object down over the dead ones, exactly like the in-place deduplication above. The write pointer trails the read pointer, so no live object is overwritten before it has been moved.

Pollard's rho factorisation uses Floyd's tortoise and hare to detect the cycle in the pseudo-random sequence it generates. That is what lets it run in constant memory instead of storing every value it has seen.

Common mistakes

Feeding it unsorted input. pair_sum([8, 1, 3], 4) returns None. Trace it: 8 + 3 = 11 is too big so right moves to index 1, then 8 + 1 = 9 is still too big so right moves to index 0, and the loop ends — even though 1 + 3 is exactly 4. No exception, just a wrong answer. Sort first, or use a set.

Using left <= right. The loop then compares an element with itself. For pair_sum([2, 5, 7], 10) that would report indices 1 and 1 as a valid pair for the target 10.

Moving both pointers on a mismatch. It feels symmetric and it is wrong. On [2, 3, 4, 5] with target 8, the correct walk goes 2 + 5 = 7 (too small, move left) then 3 + 5 = 8 and finds it. Moving both jumps straight to 3 + 4 = 7, then the pointers cross and the answer is missed.

Forgetting a branch that moves nothing. Every path through the loop body must advance a pointer or return. A condition that leaves both indices untouched is an infinite loop, and it is the failure mode you will hit first when adapting the pattern to a new problem.

Treating the in-place result as a shortened list. remove_duplicates returns a count, and the entries past that count are stale. Always slice with the returned length.

Skipping duplicate handling in 3-sum. Without the continue for a repeated fixed value, [-1, -1, 0, 1, 1] yields the triple (-1, 0, 1) twice — once for each -1 in the input. The skips are not an optimisation, they are part of the specification.

Practice

Given a sorted list and a target, count how many index pairs add to that target, including when the same value appears more than once.
Move every zero in a list to the end while keeping the order of the non-zero values, in one pass with O(1) extra space.
Merge two already-sorted lists into one sorted list using one pointer into each, without calling sorted().
Given a list sorted in ascending order that may contain negatives, return the squares of its values in sorted order in O(n) — the trick is to fill the output from the back.
Given a list of bar heights, compute the total rainwater trapped between them using two converging pointers and a running maximum from each side.

Summary

Two pointers turns a quadratic search over pairs into a linear sweep by making each step retire an entire row or column of the pairs you never have to check. The loop bound is easy — neither index backtracks, so the body runs at most n times — and the real work is the safety proof that lets you move a pointer at all. Learn the staircase argument on the sorted pair sum, then notice the same shape in palindromes, in-place compaction, cycle detection and the container problem.


Difficulty	Medium
Time (sorted input)	O(n) — the gap between the pointers shrinks by one per turn
Time (unsorted input)	O(n log n) — the sort dominates the O(n) sweep
Time (3-sum)	O(n²) — n fixed values, each costing one linear sweep
Space	O(1) — two indices and a running value; O(n) as soon as you sort or copy first
In place	Yes, for the compaction and reversal variants
Precondition	Sorted input, or some other monotonic property to move on
Data structure	List / array with O(1) indexing; linked list for fast and slow
Use it when	The answer is a pair of positions and memory must stay constant
Avoid it when	Data is unsorted and a set is affordable, or the answer is a contiguous run
Real-world use	Merge joins, LSM compaction, `str.strip()`, mark-compact garbage collection
Python equivalent	`values.reverse()`, `heapq.merge()`, `bisect` for the sorted-search half

Keep reading

The Sliding Window Technique — the other two-index pattern, for when the answer is a contiguous run rather than a pair.
Binary Search — the other way to exploit sorted input, and what pairs with two pointers when you need a partner rather than a pair.
Merge Sort — its merge step is two pointers over two sorted lists, and it is where the O(n log n) sorting cost comes from.
Linked Lists — build the structure that fast and slow pointers traverse, including cycles.
Big O Notation — the counting arguments used throughout this post, from first principles.

Fast Exponentiation in Python: Computing Huge Powers in log n Steps

Bimal Kshetri — Sun, 16 Aug 2026 16:38:19 +0000

Ask Python for 2 ** 10**18 % 1000000007 and it will never answer. Not because the modulus is hard, but because ** is evaluated first, and 2 ** 10**18 is a number with 10^18 binary digits. Writing it down needs about 125 petabytes of memory. The machine will thrash and die long before the % runs.

Ask Python for pow(2, 10**18, 1000000007) and it answers instantly. Same value, 83 multiplications, and every value it carries along the way stays below the ten-digit modulus. The gap between those two lines is one of the largest practical speedups in all of elementary computing, and it comes from a single identity you already know: a^n for even n equals (a^(n/2))^2.

Halve the exponent, square the answer. Repeat. An exponent of 10^18 reaches zero in 60 halvings, not 10^18 decrements. That is fast exponentiation — also called binary exponentiation, exponentiation by squaring, or square-and-multiply — and it is the arithmetic that every RSA handshake and every Diffie-Hellman key exchange on the internet is built from.

The idea

Start with the obvious way to compute a^n: multiply by a, n times.

def naive_power(base: int, exponent: int) -> int:
    """Multiply base by itself, one multiplication per unit of the exponent."""
    result = 1
    for _ in range(exponent):
        result *= base
    return result


print(naive_power(3, 13))
print(naive_power(2, 10))

1594323
1024

That costs n multiplications. The cost is proportional to the value of the exponent, which is exponential in the number of digits you typed. A 19-digit exponent means a quintillion multiplications.

Now the identity. For any a and any even n:

a^n = (a^(n/2))^2

Read it right to left and it is a bargain. To get a^n, you do not need a^n. You need a^(n/2), and then one squaring. Half the problem, plus one multiplication.

Odd exponents cannot be halved cleanly, so peel one factor off first:

a^n = a × a^(n-1) when n is odd

Subtracting 1 from an odd number always leaves an even number, so an odd step is immediately followed by a halving. The two rules together, plus the base case a^0 = 1, define the whole algorithm:

n is 0: the answer is 1.
n is even: compute a^(n/2), square it.
n is odd: compute a^(n-1), multiply by a.

Every two steps at worst, the exponent halves. Halving 10^18 takes 60 steps. That is the entire performance argument.

Watching it work

Take 3^13. The naive loop would do 12 multiplications. Here is what the halving rules do instead, working downwards from the top:

3^13 — 13 is odd, so 3^13 = 3 × 3^12. Need 3^12.
3^12 — even, so 3^12 = (3^6)^2. Need 3^6.
3^6 — even, so 3^6 = (3^3)^2. Need 3^3.
3^3 — odd, so 3^3 = 3 × 3^2. Need 3^2.
3^2 — even, so 3^2 = (3^1)^2. Need 3^1.
3^1 — odd, so 3^1 = 3 × 3^0. Need 3^0.
3^0 — the base case, 1.

Now unwind, filling in real numbers on the way back up: 3^0 = 1, 3^1 = 3, 3^2 = 9, 3^3 = 27, 3^6 = 729, 3^12 = 531441, 3^13 = 1594323.

Six lines, and only five of them are genuine multiplications — 3 × 3^0 is a multiplication by 1, which costs nothing in principle. Twelve multiplications became five.

The same thing, seen in binary

There is a second way to look at this that explains the name square-and-multiply, and it is the one the iterative code implements.

Write the exponent in binary. 13 is 1101, which says 13 = 8 + 4 + 1. Exponents add when powers multiply, so:

3^13 = 3^8 × 3^4 × 3^1

You can build 3^1, 3^2, 3^4, 3^8 by starting at 3 and squaring repeatedly — each squaring doubles the exponent, so four squarings reach 3^16. Then keep only the ones whose bit is set and multiply them together.

6561 × 81 × 3 = 1594323. Two multiplications to combine them, three squarings to build them: five, matching the recursion exactly. That is not a coincidence — the recursion's even steps are the squarings and its odd steps are the kept bits, read in the opposite order.

The number of squarings is the number of bits in n minus one. The number of combining multiplications is the number of 1 bits minus one. Both are at most log2(n), which is where the bound comes from.

The code

The recursive version reads exactly like the three rules.

def fast_power(base: int, exponent: int) -> int:
    """base ** exponent, computed with the halving identity.

    An even exponent squares the half power. An odd exponent peels off one
    factor of base, which always leaves an even exponent behind.
    """
    if exponent == 0:
        return 1
    if exponent % 2 == 0:
        half = fast_power(base, exponent // 2)
        return half * half
    return base * fast_power(base, exponent - 1)


print(fast_power(3, 13))
print(fast_power(2, 10))
print(fast_power(5, 0))

1594323
1024
1

Add a print and you can watch the walkthrough happen for real:

def traced_power(base: int, exponent: int) -> int:
    """The same recursion, printing each power as the calls unwind."""
    if exponent == 0:
        print(f"{base}^0 = 1")
        return 1
    if exponent % 2 == 0:
        half = traced_power(base, exponent // 2)
        print(f"{base}^{exponent} = ({base}^{exponent // 2})^2 = {half * half}")
        return half * half
    smaller = traced_power(base, exponent - 1)
    print(f"{base}^{exponent} = {base} * {base}^{exponent - 1} = {base * smaller}")
    return base * smaller


traced_power(3, 13)

3^0 = 1
3^1 = 3 * 3^0 = 3
3^2 = (3^1)^2 = 9
3^3 = 3 * 3^2 = 27
3^6 = (3^3)^2 = 729
3^12 = (3^6)^2 = 531441
3^13 = 3 * 3^12 = 1594323

The iterative version walks the exponent's bits from the bottom up. It keeps two values: result, the answer accumulated so far, and current, which holds base raised to the next power of two.

def fast_power_iterative(base: int, exponent: int, trace: bool = False) -> int:
    """Square-and-multiply, driven by the bits of the exponent from the bottom up."""
    result = 1
    current = base       # after step k this holds base ** (2 ** k)
    remaining = exponent

    while remaining > 0:
        bit = remaining % 2
        if bit == 1:
            result *= current
        if trace:
            print(f"remaining={remaining:<3} bit={bit}  current={current:<5} result={result}")
        remaining //= 2
        if remaining > 0:
            # The last squaring would never be used, so skip it.
            current *= current

    return result


print(fast_power_iterative(3, 13))
fast_power_iterative(3, 13, trace=True)

1594323
remaining=13  bit=1  current=3     result=3
remaining=6   bit=0  current=9     result=3
remaining=3   bit=1  current=81    result=243
remaining=1   bit=1  current=6561  result=1594323

How the code maps to the idea

remaining % 2 is the current low bit of the exponent, and remaining //= 2 is the halving — a right shift by one place. So the loop runs once per bit of n, which for 10^18 is 60 times.

current is the repeated-squaring chain. It starts at base and squares every iteration, so it passes through base^1, base^2, base^4, base^8 and so on. Iteration k always has current equal to base ** (2 ** k), which is exactly the value the binary decomposition wants when bit k is set.

result *= current is the "keep this one" decision. It fires only when the bit is 1, which is why the loop's multiplication count depends on how many 1 bits n has.

The order inside the loop matters. You must fold current into result before squaring current, because the bit you just examined refers to the current value, not the next one. Swapping those two lines silently produces base ** (2 * n) instead.

The final squaring is skipped by the if remaining > 0 guard. Nothing reads current after the last iteration, and that skipped squaring is the most expensive one, on the biggest numbers.

The recursion's depth is logarithmic, not linear. Each odd step is followed by an even step, so fast_power recurses at most 2 × log2(n) deep — about 120 frames for a 10^18 exponent, nowhere near Python's default limit of 1000.

Counting the multiplications directly makes the growth impossible to argue with:

def multiplication_count(exponent: int) -> int:
    """Multiplications square-and-multiply performs, without performing them."""
    count = 0
    remaining = exponent
    while remaining > 0:
        if remaining % 2 == 1:
            count += 1          # fold base ** (2 ** k) into the result
        remaining //= 2
        if remaining > 0:
            count += 1          # square to reach the next power of two
    return count


print(f"{'n':>19}  {'bits in n':>9}  {'naive':>19}  {'fast':>4}")
for exponent in (13, 100, 1000, 10 ** 6, 10 ** 18):
    print(f"{exponent:>19}  {exponent.bit_length():>9}  {exponent - 1:>19}  "
          f"{multiplication_count(exponent):>4}")

                  n  bits in n                naive  fast
                 13          4                   12     6
                100          7                   99     9
               1000         10                  999    15
            1000000         20               999999    26
1000000000000000000         60   999999999999999999    83

A quintillion multiplications become 83, and fast never exceeds twice bits in n.

Modular exponentiation

Everything above still produces enormous numbers: 7^1000 looks harmless and has 846 digits. In cryptography you never want a^n itself — you want a^n mod m, and there the numbers can be kept small at every step.

The rule that makes this work is that modular arithmetic survives multiplication:

(x × y) mod m = ((x mod m) × (y mod m)) mod m

So you can reduce after every single multiplication instead of at the end. No intermediate value ever exceeds m^2, which for a 2048-bit modulus means numbers of at most 4096 bits rather than numbers of n × 2048 bits.

def fast_power_mod(base: int, exponent: int, modulus: int) -> int:
    """base ** exponent % modulus, reducing after every multiplication."""
    if modulus == 1:
        return 0                      # everything is congruent to 0 mod 1
    result = 1
    current = base % modulus
    remaining = exponent

    while remaining > 0:
        if remaining % 2 == 1:
            result = (result * current) % modulus
        remaining //= 2
        if remaining > 0:
            current = (current * current) % modulus

    return result


print(fast_power_mod(7, 1000, 13), pow(7, 1000, 13))
print(len(str(7 ** 1000)))
print(fast_power_mod(2, 10 ** 18, 1_000_000_007), pow(2, 10 ** 18, 1_000_000_007))

9 9
846
719476260 719476260

Three details earn their place. base % modulus up front handles a base larger than the modulus. The modulus == 1 guard covers the one input the loop gets wrong: every integer is congruent to 0 modulo 1, and an exponent of 0 never enters the loop, so the initial result = 1 would come straight back out. And the reduction happens on both the accumulator and the squaring chain — miss either and the numbers grow without limit.

Use the built-in

Python has this in the language. pow(base, exponent) is the same as base ** exponent, and the three-argument pow(base, exponent, modulus) is modular exponentiation implemented in C, using a windowed variant of the same square-and-multiply algorithm. Since Python 3.8, a negative exponent with a modulus gives the modular inverse: pow(a, -1, m) is the value x with a × x ≡ 1 (mod m), and it raises ValueError when no inverse exists, which happens exactly when a and m share a factor.

print(pow(3, 13))
print(pow(3, 13, 100))

private_exponent = pow(17, -1, 3120)
print(private_exponent, (17 * private_exponent) % 3120)

try:
    pow(6, -1, 9)
except ValueError as error:
    print(f"ValueError: {error}")

print(pow(7, 11, 13), pow(7, -1, 13))

1594323
23
2753 1
ValueError: base is not invertible for the given modulus
2 2

That 2753 is not a random number: it is the RSA private exponent from the textbook key with primes 61 and 53, where the modulus is 3233 and e is 17. The last line shows the other route to an inverse — by Fermat's little theorem, a^(p-2) ≡ a^(-1) (mod p) for prime p, so pow(7, 11, 13) and pow(7, -1, 13) agree.

Complexity

Multiplications: O(log n). The loop runs once per bit of n, and the number of bits is floor(log2(n)) + 1. Each iteration does at most two multiplications: one squaring, plus one fold into the result when the bit is set. So the total is at most 2 × log2(n) + 1, hit when every bit is set, and never below log2(n), the floor reached when n is a power of two and only one bit is set. The exact count is (bits - 1) + popcount(n), and the table above shows it: 10^18 has 60 bits and 24 of them set, giving 59 + 24 = 83.

The honest caveat: a multiplication is not O(1). Big O on multiplication counts is the right measure only when each multiplication has fixed cost. For pow(a, n, m) that holds — every operand stays below m, so each multiplication costs the same, and the total is O(log n) machine-word multiplications for fixed m.

Without a modulus it does not hold at all. a^n has about n × log2(a) bits, so the final squaring alone operates on numbers half that size, and CPython multiplies large integers with Karatsuba's algorithm at roughly O(d^1.585) for d-digit operands. The cost is then dominated by the last one or two multiplications, and it is governed by the size of the answer — which no algorithm can beat, because you have to write the answer down.

Space: O(1) for the iterative version — three integers, regardless of n. The recursive version adds O(log n) stack frames. With a modulus, every one of those integers stays below m, so the space is genuinely constant in n.

Matrix exponentiation, and Fibonacci in O(log n)

Nothing in the algorithm cares that the values are numbers. It needs multiplication, associativity and an identity element — nothing else. Matrices have all three.

Linear recurrences can be written as matrix powers. The Fibonacci step F(n+1) = F(n) + F(n-1) is exactly what the matrix [[1, 1], [1, 0]] does to the vector (F(n), F(n-1)). Raise that matrix to the nth power and its top-right entry is F(n).

Matrix = list[list[int]]


def matrix_multiply(left: Matrix, right: Matrix) -> Matrix:
    """Multiply two 2x2 integer matrices."""
    return [
        [left[0][0] * right[0][0] + left[0][1] * right[1][0],
         left[0][0] * right[0][1] + left[0][1] * right[1][1]],
        [left[1][0] * right[0][0] + left[1][1] * right[1][0],
         left[1][0] * right[0][1] + left[1][1] * right[1][1]],
    ]


def matrix_power(matrix: Matrix, exponent: int) -> Matrix:
    """The same square-and-multiply loop, with matrices instead of numbers."""
    result = [[1, 0], [0, 1]]        # the identity matrix plays the role of 1
    current = matrix
    remaining = exponent

    while remaining > 0:
        if remaining % 2 == 1:
            result = matrix_multiply(result, current)
        remaining //= 2
        if remaining > 0:
            current = matrix_multiply(current, current)

    return result


def fibonacci(index: int) -> int:
    """F(index), using O(log index) matrix multiplications."""
    return matrix_power([[1, 1], [1, 0]], index)[0][1]


print([fibonacci(n) for n in range(11)])
print(fibonacci(90))
print(len(str(fibonacci(1_000_000))))

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
2880067194370816120
208988

The loop body is identical to fast_power_iterative with * replaced by matrix_multiply and 1 replaced by the identity matrix. The last line computes the millionth Fibonacci number — 208,988 digits long — with 26 matrix multiplications, the count the table above gives for an exponent of 10^6. The usual iterative Fibonacci would need a million additions on numbers of that size.

The same construction handles any recurrence of the form T(n) = c1×T(n-1) + ... + ck×T(n-k): build the k-by-k companion matrix and raise it to the nth power in O(k^3 log n) operations.

When to use it, and when not to

Use it whenever the exponent is large and the modulus is present. That is the entire domain: cryptography, hashing, counting problems modulo a prime like 10^9 + 7, and linear recurrences via matrices. Under a modulus there is no reason ever to write the naive loop.

Do not write it yourself in Python. pow(base, exponent, modulus) does the same thing in C, and beats a Python-level loop by a wide margin. The value in the code above is knowing what pow costs and why, not shipping it.

Do not use it for cryptographic code you deploy. The textbook loop branches on the bits of the exponent, so its running time and power draw leak the secret exponent — this is the basis of Paul Kocher's 1996 timing attacks on Diffie-Hellman and RSA. Real implementations use constant-time variants such as the Montgomery ladder; OpenSSL's is BN_mod_exp_mont_consttime. Use a vetted library.

Do not bother for small fixed exponents. x * x * x beats any general routine for a cube.

Do not use it to build an exact huge power you then throw away. If you only need the last few digits or the size of the answer, use a modulus or math.log instead. 2 ** 10**9 is a valid Python expression and a bad idea.

Where it shows up in the real world

RSA. Both encryption and decryption are single modular exponentiations: the ciphertext is m^e mod n and the plaintext is c^d mod n. The common public exponent 65537 keeps the encrypting side cheap: it is 2^16 + 1 — binary 10000000000000001, two bits set — so encryption costs 16 squarings and one multiply. The private exponent is a full-length 2048-bit number, which is why signing is far slower than verifying.

Diffie-Hellman and its elliptic-curve variant. Finite-field Diffie-Hellman computes g^a mod p with a secret exponent. The elliptic-curve version, which is what a modern TLS handshake normally negotiates, runs the identical algorithm with point addition in place of multiplication, where it is called double-and-add.

Primality testing. The Miller-Rabin and Fermat tests are built on computing a^(n-1) mod n for candidate n. Every large prime used in cryptography was found by running modular exponentiation a few hundred times.

Rolling hashes. Rabin-Karp and Rabin fingerprints need base^(window - 1) mod m to drop the outgoing character. With a window of 10,000 that is 21 multiplications instead of nearly 10,000.

Modular inverses. RSA key generation derives the private exponent d as the inverse of e modulo (p-1)(q-1), which is one pow(e, -1, phi) call.

Linear recurrences. Matrix exponentiation is the standard way to answer "the nth term of this recurrence, modulo p, for n up to 10^18".

Common mistakes

Computing the power first, then the modulus. (a ** b) % m is the single most common bug here. It is mathematically correct and computationally hopeless — Python builds the entire a^b before reducing. Always pow(a, b, m).

Forgetting to reduce the squaring chain. Reducing result but not current keeps the answer correct and lets current grow to the full unreduced size, which throws away the whole point.

Squaring before multiplying into the result. Inside the loop, result *= current must come before current *= current. Reversed, every set bit contributes one power too many.

Returning 1 for a modulus of 1. pow(x, y, 1) is 0, because every integer is congruent to 0 modulo 1. A hand-written version that starts result = 1 and never multiplies returns 1 instead. Guard it.

Negative exponents. fast_power(2, -3) never reaches exponent == 0 with the halving rules as written, because -3 is odd and -4 // 2 is -2, then -1, then -2 again — it recurses until Python raises RecursionError. The iterative version fails more quietly: while remaining > 0 never runs, so it returns 1. Either raise for negative exponents or convert them to 1 / fast_power(base, -exponent).

Using math.pow. It returns a float and loses exactness above 2^53. math.pow(3, 40) is not an integer answer. Use ** or pow for integers.

Practice

Rewrite fast_power so it raises a ValueError on a negative exponent, and confirm the error fires.
Write the top-down version that walks bin(exponent) from the leading bit, squaring the result each step and multiplying by the base when the bit is 1, then check it agrees with fast_power_iterative.
Compute the last ten digits of 2^(10^18) using a modulus of 10^10, and explain why the answer is even.
Implement a modular inverse using Fermat's little theorem for a prime modulus, and compare it against pow(a, -1, m) for every a from 1 to 100 with modulus 101.
Use matrix_power to compute the nth term of T(n) = 2×T(n-1) + 3×T(n-2) with T(0) = 0 and T(1) = 1, modulo 10^9 + 7, for n equal to one million.

Summary

Fast exponentiation is one identity applied repeatedly: halve the exponent and square, peel a factor off when the exponent is odd. That turns n multiplications into (bits - 1) + popcount(n), which is never more than 2 × log2(n) + 1. Reduce modulo m after every multiplication and the numbers stay small too, which is what makes public-key cryptography possible at all. In Python you write pow(base, exponent, modulus) — but knowing what those three arguments cost is the difference between code that returns in a microsecond and code that never returns.


Difficulty	Medium
Best case	O(log n) — n a power of two: log2(n) squarings and one fold
Average case	O(log n) multiplications
Worst case	O(log n) multiplications — at most `2 × log2(n) + 1`, when every bit is set
Space	O(1) iterative — three integers; O(log n) stack frames if recursive
Exact cost	`(bits in n − 1)` squarings plus `popcount(n)` folds
Multiplication cost	Constant only under a modulus; without one, big-int multiplication dominates
Data structure	Integers, or anything with associative multiplication and an identity
Generalises to	Matrices (linear recurrences), polynomials, elliptic-curve points
Use it when	The exponent is large and you want it modulo something
Avoid it when	The exponent is small and fixed, or the exact unreduced power is astronomically large
Real-world use	RSA, Diffie-Hellman, Miller-Rabin primality, rolling hashes, modular inverses
Python equivalent	`pow(base, exponent, modulus)`, and `pow(a, -1, m)` for inverses since 3.8

Keep reading

The Euclidean Algorithm — the other halving-style number-theory algorithm, and the extended version that also computes modular inverses.
Binary Search — the same "halve the problem" move applied to a sorted list instead of an exponent.
Sieve of Eratosthenes — where the primes used as moduli come from.
Rabin-Karp String Matching — rolling hashes, which need a modular power of the base.
Big O Notation — why log n and n are separated by more than most people expect.

Sieve of Eratosthenes in Python: Every Prime Under a Million, Fast

Bimal Kshetri — Sun, 16 Aug 2026 16:38:13 +0000

There are 78,498 prime numbers below one million, and a laptop can produce every one of them in a few milliseconds. It can only do that if you stop asking the obvious question.

The obvious question is "is this number prime?", asked once per number. Answer it by trial division and listing the primes below a million costs about 68 million modulo operations. The sieve of Eratosthenes asks a different question — "which numbers can I rule out?" — and answers it with 2.1 million writes and not a single division. Same list of primes, thirty-two times less work, and the gap widens with every extra digit.

Eratosthenes ran the library at Alexandria in the third century BC and is better known for measuring the circumference of the Earth. His method for primes is 2,200 years old and is still what production code uses when it needs every prime up to a bound. What follows proves the two optimisations that most implementations write without justifying, derives the O(n log log n) running time from the sum of the reciprocals of the primes rather than asserting it, and covers the two variants you will actually reach for: a segmented sieve for ranges too large to hold in memory, and a smallest-prime-factor table that turns factorisation into a lookup.

The idea

Write out every whole number from 1 to n. Cross out 1 straight away — a prime has exactly two distinct divisors and 1 has one, so it does not qualify.

Now repeat two steps until you run out of work:

Find the smallest number still standing. It is prime.
Cross out every multiple of it.

Step 1 is the part that deserves an argument, and the argument is three lines. Suppose k is the smallest number above 1 that is still standing. Every prime below k was picked in an earlier round, and each of those rounds crossed out all of its own multiples. So k is not a multiple of any prime smaller than itself. A composite number always has a prime factor smaller than itself. Therefore k is not composite.

Notice what step 2 does not involve. Crossing out the multiples of 7 means visiting 14, 21, 28, 35 — repeated addition with a stride of 7. There is no division anywhere in the algorithm and no test of any individual number. Trial division asks n separate questions, each costing up to the square root of n in divisions; the sieve replaces all of them with a handful of straight walks through an array.

Watching it work

Take the numbers 1 to 30. Cross out 1, then start.

Round 1 — the smallest survivor is 2, so 2 is prime. Cross out 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30. That is 14 cells written.

Round 2 — the smallest survivor is now 3. Its first multiple worth visiting is 3 × 3 = 9, because 6 already went with the twos. Cross out 9, 12, 15, 18, 21, 24, 27, 30 — 8 cells, and half of them (12, 18, 24, 30) were already crossed out by 2. The sieve does not check first. Reading a cell to find out whether it is already marked costs at least as much as writing it, and writing "already crossed out" again is harmless.

Round 3 — the smallest survivor is 5. Cross out 25 and 30. Only 25 is new; 10, 15 and 20 went in earlier rounds and 30 went twice already.

Round 4 — the smallest survivor is 7, and this is where you stop. 7 × 7 = 49 is already past 30, and every multiple of 7 beyond 7 itself that is 30 or less is 7 times something smaller than 7, so it carries a factor below 7 and went long ago. The same holds for 11, 13 and every later survivor: nothing is left to mark, so everything still standing is prime.

Ten primes — 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 — for 24 writes and zero divisions.

The six-column layout is not decorative: every prime past 3 lands in the first or the fifth column, because the other four columns are exactly the numbers divisible by 2 or by 3. Past 5, only 8 of every 30 consecutive numbers can be prime at all — the observation that wheel factorisation uses to shrink a sieve's memory.

The code

Start with the baseline, so the improvement is measurable rather than asserted.

from math import isqrt, log


def is_prime_by_trial_division(number: int) -> bool:
    """True if number is prime, decided by dividing by every candidate up to its square root."""
    if number < 2:
        return False
    divisor = 2
    # If number = a * b with a <= b, then a * a <= number, so any divisor above
    # the square root can only appear paired with one below it.
    while divisor * divisor <= number:
        if number % divisor == 0:
            return False
        divisor += 1
    return True


def primes_by_trial_division(limit: int) -> list[int]:
    """Every prime below limit, decided one number at a time."""
    return [number for number in range(2, limit) if is_prime_by_trial_division(number)]


print(primes_by_trial_division(30))
print(len(primes_by_trial_division(1000)))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
168

That is correct and it is the right tool when you have exactly one number to test. As a way to produce a whole list it throws away everything it learns: the moment it discovers that 4 is even, it forgets, and asks the same question again at 6.

Here is the sieve. limit is exclusive, matching range, so sieve_of_eratosthenes(30) returns the primes below 30.

def sieve_of_eratosthenes(limit: int) -> list[int]:
    """Return every prime strictly below limit, in ascending order.

    Rather than testing numbers, this marks composites: each prime it finds
    strikes off that prime's own multiples in a single sweep of additions.
    """
    if limit <= 2:
        return []

    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False

    candidate = 2
    # Every composite below limit has a prime factor no larger than its own
    # square root, so once candidate squared reaches limit nothing is left to mark.
    while candidate * candidate < limit:
        if is_prime[candidate]:
            # Multiples of candidate below candidate squared already carry a
            # smaller prime factor and were struck off on an earlier sweep.
            for multiple in range(candidate * candidate, limit, candidate):
                is_prime[multiple] = False
        candidate += 1

    return [number for number, prime in enumerate(is_prime) if prime]


print(sieve_of_eratosthenes(30))
print(sieve_of_eratosthenes(2), sieve_of_eratosthenes(3), sieve_of_eratosthenes(0))
print(sieve_of_eratosthenes(10_000) == primes_by_trial_division(10_000))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
[] [2] []
True

Twenty-two lines, and the last one checks it against the baseline on every number below ten thousand.

In CPython the version above spends nearly all its time in the interpreter, one bytecode dispatch per write. Because the marking pattern is an arithmetic progression, you can hand the entire inner loop to C with a single slice assignment.

def sieve_flags(limit: int) -> bytearray:
    """The same sieve with the marking loop handed to C by slice assignment."""
    if limit <= 2:
        return bytearray(max(limit, 0))

    flags = bytearray(b"\x01") * limit
    flags[0:2] = b"\x00\x00"

    for candidate in range(2, isqrt(limit - 1) + 1):
        if flags[candidate]:
            start = candidate * candidate
            # bytearray(k) is k zero bytes, so one assignment clears the whole
            # arithmetic progression without a Python-level loop.
            flags[start::candidate] = bytearray(len(range(start, limit, candidate)))

    return flags


flags = sieve_flags(1_000_000)
print(sum(flags), flags[999_983], flags[999_984])

78498 1 0

Identical algorithm, identical number of writes, roughly forty times faster on a million-element sieve, and one byte per number instead of the eight a list of pointers costs. Converting the flags back into a list of primes then costs more than the sieve itself — if the question is "is 786,433 prime?", keep the flags and index into them.

How the code maps to the idea

The array is the sheet of paper, and the index is the number. is_prime[17] answers a question about 17: no hashing, no search, no comparison, just an offset into a contiguous block of memory. That is the decision the whole algorithm rests on, and it is also where the O(n) space cost comes from.

is_prime[0] = is_prime[1] = False handles the two numbers the loop never visits. Neither is prime and nothing crosses them out.

The outer loop stops at the square root. Any composite m can be written as a * b with a <= b, which forces a * a <= m, so m has a prime factor no larger than the square root of m — and therefore smaller than the square root of limit. Some earlier round has already crossed m out. Past that point every remaining round would sweep an empty range.

if is_prime[candidate] skips composites, and loses nothing by doing so. A composite candidate has a prime factor q smaller than itself, and every multiple of candidate is also a multiple of q, so q already dealt with all of them.

The inner loop starts at candidate * candidate, not 2 * candidate. Take any earlier multiple k * candidate with 2 <= k < candidate. Then k has a prime factor q with q <= k < candidate, so q ran in an earlier round, marking multiples of q from q * q upwards. And k * candidate is at least q * candidate, which exceeds q * q because candidate exceeds q. So it was already struck off: skipping it is correct, not merely convenient.

The inner loop step is candidate, which keeps the walk sequential in memory and free of arithmetic beyond one addition per cell.

Edge cases fall out of the single guard at the top. Limits of 0, 1 and 2 return an empty list before any array is allocated; a limit of 3 clears the first two flags, never enters the outer loop because 2 × 2 is not below 3, and answers [2].

Now count what the two approaches actually do, rather than trusting the description.

def count_trial_division(limit: int) -> int:
    """Modulo tests trial division performs to list every prime below limit."""
    tests = 0
    for number in range(2, limit):
        divisor = 2
        while divisor * divisor <= number:
            tests += 1
            if number % divisor == 0:
                break
            divisor += 1
    return tests


def count_sieve(limit: int, square_start: bool = True, stop_at_root: bool = True) -> int:
    """Cells the sieve writes to, with either optimisation switchable off."""
    writes = 0
    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False
    outer_bound = isqrt(limit - 1) + 1 if stop_at_root else limit
    for candidate in range(2, outer_bound):
        if is_prime[candidate]:
            start = candidate * candidate if square_start else 2 * candidate
            for multiple in range(start, limit, candidate):
                is_prime[multiple] = False
                writes += 1
    return writes


print(f"{'limit':>9}  {'modulo tests':>14}  {'sieve writes':>14}  {'ratio':>6}")
for limit in (1_000, 10_000, 100_000):
    tests = count_trial_division(limit)
    writes = count_sieve(limit)
    print(f"{limit:>9,}  {tests:>14,}  {writes:>14,}  {tests / writes:>6.1f}")

    limit    modulo tests    sieve writes   ratio
    1,000           5,287           1,409     3.8
   10,000         117,526          16,979     6.9
  100,000       2,745,693         193,076    14.2

The ratio roughly doubles for every tenfold increase in n, which is what you expect when one side grows like the square root of n and the other barely grows at all. Extend the same measurement to a million and trial division needs 67,740,403 modulo tests against the sieve's 2,122,046 writes — and a modulo is a much more expensive machine instruction than a byte store.

Complexity

Time: O(n log log n), best, average and worst alike. There is no input to be lucky or unlucky with: the only parameter is the bound, so a given n always writes to exactly the same cells. The three cases collapse into one count, and here is where that count comes from.

For each prime p with p * p < n, the inner loop writes to the multiples of p from p * p up to n. That is (n - p * p) / p + 1 writes, which is a little under n / p. Add it up over every prime the outer loop reaches:

total writes  ≈  n × (1/2 + 1/3 + 1/5 + 1/7 + ... + 1/p)   for primes p ≤ √n

So the entire running time hangs on one quantity: the sum of the reciprocals of the primes up to a bound. That sum diverges, but astonishingly slowly. Mertens proved in 1874 that

1/2 + 1/3 + 1/5 + ... + 1/p  =  ln ln x + M + o(1)     for primes p ≤ x

where M ≈ 0.2615 is the Meissel-Mertens constant. The double logarithm is not a coincidence, and you can see where it comes from without the full proof. The prime number theorem says there are about x / ln x primes below x, which makes the k-th prime roughly k ln k in size. Summing reciprocals of the primes is then summing 1 / (k ln k) over k, and that sum behaves like the integral of 1 / (k ln k), which is ln ln k. One logarithm comes from the thinning-out of the primes; the other comes from integrating a reciprocal.

Substituting x = √n gives ln ln √n = ln(ln n / 2) = ln ln n - ln 2, so the count of writes is about n × (ln ln n - ln 2 + M), and dropping the constants leaves O(n log log n).

That is a claim you can check.

MEISSEL_MERTENS = 0.2614972128


def predicted_writes_per_number(limit: int) -> tuple[float, float]:
    """(plain n log log n estimate, same estimate corrected for the square start)."""
    base_primes = sieve_of_eratosthenes(isqrt(limit - 1) + 1)
    # Mertens: the sum of 1/p over primes p <= x is ln ln x + M.
    plain = log(log(limit)) - log(2) + MEISSEL_MERTENS
    # The square start skips the first p - 1 multiples of every base prime p.
    corrected = plain - (sum(base_primes) - len(base_primes)) / limit
    return plain, corrected


print(f"{'n':>9}  {'measured':>8}  {'ln ln n - ln 2 + M':>19}  {'corrected':>9}")
for limit in (10_000, 100_000, 1_000_000):
    measured = count_sieve(limit) / limit
    plain, corrected = predicted_writes_per_number(limit)
    print(f"{limit:>9,}  {measured:>8.3f}  {plain:>19.3f}  {corrected:>9.3f}")

        n  measured   ln ln n - ln 2 + M  corrected
   10,000     1.698                1.789      1.685
  100,000     1.931                2.012      1.920
1,000,000     2.122                2.194      2.118

The last column matches the measurement to within a fifth of a percent at a million. The formula is not an analogy; it is the count.

It is also worth seeing how small log log n is. At n = 1,000,000 the sieve does 2.12 writes per number. At n = 1,000,000,000,000 it does about 2.9. For any n you will ever sieve, log log n sits between 1 and 4 — which is why people describe the sieve as "basically linear" and are not being sloppy when they do.

The baseline it beats: O(n √n), or better in practice. Trial division tests each of n numbers with up to √n divisions, giving the plain O(n √n) bound. The measured cost is lower, because composites exit at their first divisor and only the primes pay the full square root — about n / ln n of them, for a true cost near n^1.5 / ln n. Either way it is a power of n, and the sieve is a logarithm of a logarithm.

Neither of the two famous optimisations changes the complexity class. They are constant-factor wins, and it is worth knowing how large.

print("writes at n = 1,000,000, by which optimisations are switched on")
print(f"  both:                 {count_sieve(1_000_000):>9,}")
print(f"  no square start:      {count_sieve(1_000_000, square_start=False):>9,}")
print(f"  neither:              {count_sieve(1_000_000, False, False):>9,}")

writes at n = 1,000,000, by which optimisations are switched on
  both:                 2,122,046
  no square start:      2,197,837
  neither:              2,775,208

Starting at p * p instead of 2 * p saves 3.6% of the writes — real, but far less than its fame suggests, and the share shrinks as n grows. Running the outer loop all the way to n while starting at 2 * p costs 31% more, because you are then summing 1/p over every prime below n rather than below √n, and the difference between those two sums tends to ln 2. If you keep the square start, stopping at the square root saves no writes at all — the ranges past that point are empty — it saves the loop overhead of visiting the other n cells.

Space: O(n). One flag per number, whether or not you ever look at it. With a bytearray that is one byte per number: 1 MB for a million, 1 GB for a billion. A bit array cuts it eightfold and a wheel that skips multiples of 2, 3 and 5 cuts it by another factor of 30/8, but the shape of the cost does not change. This is the sieve's real limitation, and it is the reason the next section exists.

Two variants worth the extra code

Sieving a range that will not fit in memory

Suppose you want the primes just above one trillion. A flat sieve would need a trillion cells. But to sieve any block of numbers up to high, all you need are the primes up to √high — because that is the largest prime factor a composite in that range can be forced to have. For a trillion that is the 78,498 primes below a million, which fit comfortably in memory.

So: sieve the small primes once, then sweep the target range one block at a time, marking each block with those same primes.

def segmented_sieve(low: int, high: int) -> list[int]:
    """Every prime in the half-open range [low, high), one block at a time.

    Memory holds the primes up to the square root of high plus a single block,
    so high may be far beyond anything a flat sieve could store.
    """
    if high <= 2:
        return []
    low = max(low, 2)
    base_primes = sieve_of_eratosthenes(isqrt(high - 1) + 1)

    found: list[int] = []
    block_size = 1 << 16
    for block_start in range(low, high, block_size):
        block_stop = min(block_start + block_size, high)
        flags = bytearray(b"\x01") * (block_stop - block_start)

        for prime in base_primes:
            # First multiple of prime at or above block_start, but never below
            # prime squared, or a base prime would strike itself out.
            first = max(prime * prime, -(-block_start // prime) * prime)
            for multiple in range(first, block_stop, prime):
                flags[multiple - block_start] = 0

        found.extend(block_start + offset for offset, flag in enumerate(flags) if flag)
    return found


print(segmented_sieve(0, 30) == sieve_of_eratosthenes(30))
print(segmented_sieve(1_000_000, 1_000_100))
print(len(segmented_sieve(10**12, 10**12 + 10_000)))

True
[1000003, 1000033, 1000037, 1000039, 1000081, 1000099]
335

There are 335 primes in the ten thousand numbers starting at a trillion, found without ever allocating more than 65,536 flags plus the base primes. Working memory is O(√high + block size), plus whatever share of the output you keep. The max(prime * prime, ...) term is the same square-start rule as before, and it is what stops the block containing 2 from crossing 2 out.

Two costs to be aware of. Every block pays a pass over all 78,498 base primes even when most of them have no multiple inside it, so very small blocks are wasteful — 32 KB to 256 KB is the usual choice, sized to fit the CPU cache, which is also where the speed comes from. And you cannot ask for primes near 10^30 this way: √high grows too, and the base sieve becomes the bottleneck.

Sieving the smallest prime factor

Change one thing: instead of storing a flag, store the smallest prime that divides each index. The loop is the same shape, plus one condition so that the first prime to reach a cell keeps it.

def smallest_prime_factor_sieve(limit: int) -> list[int]:
    """spf[number] is the smallest prime dividing number, for every number below limit."""
    spf = list(range(limit))
    candidate = 2
    while candidate * candidate < limit:
        if spf[candidate] == candidate:  # nothing smaller divides it, so it is prime
            for multiple in range(candidate * candidate, limit, candidate):
                if spf[multiple] == multiple:  # not yet claimed by a smaller prime
                    spf[multiple] = candidate
        candidate += 1
    return spf


def factorise(number: int, spf: list[int]) -> list[int]:
    """Prime factors of number with multiplicity. Needs 1 <= number < len(spf)."""
    factors: list[int] = []
    while number > 1:
        prime = spf[number]
        factors.append(prime)
        number //= prime  # peel off one copy of that prime and continue
    return factors


spf = smallest_prime_factor_sieve(1_000_000)
print(factorise(999_999, spf))
print(factorise(999_983, spf))
print(factorise(360, spf), factorise(1, spf))

[3, 3, 3, 7, 11, 13, 37]
[999983]
[2, 2, 2, 3, 3, 5] []

A cell keeps the first prime that reaches it, and rounds run in increasing order of prime, so that first arrival is the smallest prime factor. For any number above 1, spf[number] == number is then exactly the test for primality; spf[0] and spf[1] are the two cells no round ever claims, so exclude them. That is why the table also tells you 999,983 is prime.

Factorising afterwards is O(log n): every step divides the number by a prime of at least 2, so it cannot take more than log2(n) steps — at most 20 for a number below a million. Trial-division factorisation of the same number takes up to √n divisions, so about a thousand. If you need to factorise many numbers below a fixed bound, this table is the answer, and the same sweep can build Euler's totient, the Möbius function or a divisor count at no extra asymptotic cost.

The price is memory: a Python list of a million distinct integers is roughly 40 MB, against 1 MB for the bytearray of flags. array("i", range(limit)) from the standard library brings it down to 4 MB.

When to use it, and when not to

Use it when you need every prime below a fixed bound, or a per-number table derived from factorisation, and that bound fits in memory. This is precomputation: pay O(n log log n) once at startup, then answer any number of "is it prime" or "factorise it" questions in constant or logarithmic time.

Do not use it to test a single number, especially a large one. Asking whether one 100-digit number is prime by sieving is not slow, it is impossible — the array would need more cells than there are atoms in the observable universe. The right tool is the Miller-Rabin test, which is a handful of modular exponentiations (see fast exponentiation for the log-time squaring trick it depends on) and is deterministic for every 64-bit input with a fixed set of seven bases. For one modest number, trial division by 2, 3 and then 6k ± 1 is fine and needs no memory at all.

Do not use a flat sieve for a narrow window far out. Primes between 10^12 and 10^12 + 10^4 want the segmented sieve above, which costs O(√high) memory instead of O(high).

Be careful with the "linear sieve". There is a well-known O(n) variant that marks each composite exactly once by pairing it with its smallest prime factor. The bound is genuinely better, but in practice it is often slower than the classic sieve, because it accesses memory out of order and loses the sequential cache behaviour that makes the plain version fast. Measure before you adopt it.

Python's standard library has no prime function at all. math.isqrt and math.gcd are there; primality and factorisation are not. If you would rather not maintain this yourself, SymPy provides sympy.sieve, sympy.primerange, sympy.isprime and sympy.factorint.

Where it shows up in the real world

RSA and Diffie-Hellman key generation. Producing a 2048-bit RSA key means finding two 1024-bit primes, and the practical way to find one is to test random odd candidates until one passes. A Miller-Rabin round on a 1024-bit number is a full modular exponentiation — thousands of times more expensive than a small-division check. So implementations sieve first. OpenSSL ships a static table of the first 2,048 primes (the largest is 17,863), and its candidate generator computes the candidate's remainder against every one of them, then walks forward through candidates by updating those remainders — a segmented sieve over the candidate window. By Mertens' third theorem, that table alone eliminates about 88% of odd candidates before a single expensive test runs.

Competitive programming and Project Euler. "Find the sum of all the primes below two million" is Project Euler problem 10, and it is a two-line answer with a sieve, against about a minute of waiting with the trial division above. Contest solutions routinely sieve to 10^6 or 10^7 during setup and then answer thousands of queries by lookup; the smallest-prime-factor table is the standard way to factorise a hundred thousand inputs inside a time limit.

Hash table sizing. The GNU C++ standard library picks the bucket count for unordered_map from a hard-coded list of primes, growing to the next entry when the load factor is exceeded. A prime bucket count keeps a poorly distributed hash function from collapsing onto a few buckets. The list is a compile-time constant rather than something computed at start-up — the kind of one-off table a sieve exists to produce.

Prime counting and number theory software. Kim Walisch's primesieve — a heavily optimised segmented, wheel-factorised sieve — is the standard tool for generating or counting primes anywhere in the 64-bit range, and is what research code reaches for instead of writing its own.

Common mistakes

Starting the inner loop at 2 * candidate and calling it a bug. It is not a bug, it is 3.6% more writes. The actual bug is starting at 2 * candidate and running the outer loop to n, which costs 31% more and is what most first attempts do.

Testing before writing. if is_prime[multiple]: is_prime[multiple] = False looks like it avoids work. It adds a read to every write and avoids nothing, because the write is idempotent.

Using a set of composites instead of a flat array. A set turns each mark into a hash computation and a probe into scattered memory, and costs an order of magnitude more space. The sieve is fast precisely because the index is the number and the writes are sequential.

Off-by-one at the limit. Decide whether your limit is inclusive or exclusive and put it in the docstring. range-style exclusive is the convention above; sieve_of_eratosthenes(30) therefore does not consider 30 itself.

Forgetting 0 and 1. [True] * limit claims both are prime. Clear them explicitly; nothing in the loop ever visits them.

Sieving inside a loop. Building a fresh sieve every time you need a primality check turns an O(n log log n) precomputation into a per-query cost. Sieve once, keep the array, index into it.

Practice

Sum every prime below two million, and check that the total is what you would expect from a table of prime sums.
Count the twin primes below one million — pairs p and p + 2 that are both prime.
Extend the sieve so that it also records, for every number below the limit, how many distinct prime factors it has.
Write nth_prime(k) by sieving up to the bound k * (ln k + ln ln k), which is known to exceed the k-th prime for k of at least 6, then indexing the result.
Use the segmented sieve to find the first gap of at least 100 between consecutive primes, and report where it starts.

Summary

The sieve of Eratosthenes wins by inverting the question. Testing numbers costs a square root per number; ruling out multiples costs one write per composite discovery, and the total number of those writes is n times the sum of the reciprocals of the primes up to √n — which Mertens tells us is ln ln n - ln 2 + M, a number between 1 and 4 for every input you will ever use. That is the whole argument, and the measured writes agree with it to a fraction of a percent.


Difficulty	Medium
Time	O(n log log n) in every case — n times the sum of 1/p over primes p up to `√n`
Space	O(n) — one flag per number; 1 MB per million with a `bytearray`
Baseline it replaces	O(n √n) trial division — 32 times more work at n = 1,000,000
Divisions performed	None — the inner loop only adds
Data structure	Flat array of flags, indexed by the number itself
Range variant	Segmented sieve — O(√high) memory, any window up to 64 bits
Factorisation variant	Smallest-prime-factor table — O(log n) factorisation afterwards
Use it when	You need every prime below a fixed n, or a per-number factor table
Avoid it when	Testing one large number — use Miller-Rabin instead
Real-world use	Small-prime rejection in RSA key generation, contest precomputation
Python equivalent	Nothing in the standard library; `sympy.sieve` / `sympy.primerange`

Write it once, keep the array, and stop asking numbers whether they are prime.

Keep reading

The Euclidean Algorithm — the other 2,000-year-old algorithm still in daily use, and the one behind RSA key setup.
Fast Exponentiation — modular exponentiation in log n steps, which is what Miller-Rabin runs after the sieve has done the cheap rejections.
Prefix Sums — the same trade in a different shape: precompute once, answer every query in constant time.
Big O Notation — how to build the counting arguments used above from scratch.
Arrays and Dynamic Arrays — why indexing a contiguous block is the fastest lookup there is, which is the sieve's whole advantage.

The Euclidean Algorithm in Python: The Oldest Algorithm Still in Daily Use

Bimal Kshetri — Sun, 16 Aug 2026 16:37:41 +0000

Euclid wrote this one down around 300 BC, in Book VII of the Elements. It is still the code that runs when Python reduces a fraction, and still the code that runs when a machine generates an RSA key. No other algorithm in this series is that old and that current at the same time.

The problem it solves is small: given two whole numbers, find the largest number that divides both of them exactly. The greatest common divisor of 48 and 18 is 6. The obvious approach is to try every candidate from the smaller number downwards until one divides both, which for two eighteen-digit numbers means up to a quintillion tests. Euclid's method answers the same question in at most 90 divisions, and usually far fewer.

That gap — from a quintillion to under a hundred — comes from a single observation about divisors that takes two lines to prove. Once you have it, least common multiples, modular inverses and RSA private keys all fall out of the same loop.

The idea

A common divisor of a and b is a number that divides both with no remainder. The common divisors of 48 and 18 are 1, 2, 3 and 6, so the greatest common divisor — the gcd — is 6.

Euclid's insight is that you can shrink the problem without changing the answer.

Take two positive numbers with a larger than b, and suppose some number d divides both. Then d divides their difference a - b too — a whole number of ds minus a whole number of ds is still a whole number of ds. The argument runs backwards just as easily: anything dividing b and a - b divides their sum, which is a.

So the pair (a, b) and the pair (a - b, b) have exactly the same set of common divisors. Not a similar set — the same one. Identical sets have identical largest members:

gcd(a, b) = gcd(a - b, b)

That is the whole algorithm. Keep replacing the larger number with the difference. The numbers shrink, the answer never moves, and when the two finally become equal you are looking at the gcd, because the greatest divisor a number shares with itself is the number.

From subtracting to dividing

Repeated subtraction is correct but wasteful: the gcd of 1,000,000 and 2 costs half a million identical subtractions.

Look at what those subtractions achieve. They strip out copies of b until what is left is smaller than b — which is division with remainder, spelled slowly. Doing it in one step gives Euclid's real rule:

gcd(a, b) = gcd(b, a mod b)

where a mod b is the remainder after dividing a by b. When the remainder hits zero, b divides a exactly, and the gcd of the pair is b.

Two properties make this a real algorithm rather than a nice identity.

It is correct at every step. Write a = q * b + r. The same divisor argument still applies, with q copies of b removed instead of one: r = a - q * b is a difference of multiples of any common divisor d, and q * b + r = a is a sum of them. Same set of common divisors, same greatest member.

It always terminates. The remainder is never negative and is always strictly smaller than b, so the second number of the pair strictly decreases every step and is bounded below by zero. A strictly decreasing sequence of non-negative integers has to reach 0. It cannot take longer than b steps; it will turn out to take dramatically fewer.

Watching it work

Start with the subtraction form on 48 and 18, since it is the version you can follow without doing any arithmetic in your head.

48 and 18 — 48 is larger, so replace it: 48 - 18 = 30. Now (30, 18).
30 and 18 — replace again: 30 - 18 = 12. Now (12, 18).
12 and 18 — now 18 is the larger one: 18 - 12 = 6. Now (12, 6).
12 and 6 — 12 - 6 = 6. Now (6, 6).
The two are equal. The answer is 6.

Four subtractions. Every intermediate pair — (30, 18), (12, 18), (12, 6), (6, 6) — has the same common divisors as the original: 1, 2, 3 and 6.

Now the division form on a harder pair, 1071 and 462. By subtraction this takes 11 steps. By division it takes three:

1071 = 2 * 462 + 147, so the remainder is 147. The pair becomes (462, 147).
462 = 3 * 147 + 21, remainder 21. The pair becomes (147, 21).
147 = 7 * 21 + 0, remainder 0. Done — the gcd is 21.

Notice the shift that happens between rows: the old b becomes the new a, and the remainder becomes the new b. That single move is the entire loop body. And notice the last non-zero remainder, 21 — that is always the answer. Check it: 1071 = 21 × 51 and 462 = 21 × 22, and 51 and 22 share nothing.

The code

Here is the subtraction version first, because it maps one-to-one onto the argument above.

def gcd_subtraction(a: int, b: int) -> int:
    """Greatest common divisor by repeated subtraction.

    Both inputs must be positive integers. Replacing the larger number with
    the difference leaves the set of common divisors untouched, so the answer
    never changes; the two numbers just get smaller until they meet.
    """
    while a != b:
        if a > b:
            a -= b
        else:
            b -= a
    return a


print(gcd_subtraction(48, 18))
print(gcd_subtraction(1071, 462))

6
21

Correct, but it is only ever a teaching version — feed it 1,000,000 and 2 and it performs 499,999 subtractions, and feed it a zero and it never returns at all.

The division form fixes both problems and is shorter:

def gcd(a: int, b: int) -> int:
    """Greatest common divisor by Euclid's division rule.

    Repeatedly replaces the pair (a, b) with (b, a mod b). The second value
    strictly shrinks every step, so it reaches 0, and the first value is the
    answer at that point. gcd(0, 0) is 0 by convention, matching math.gcd.
    """
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a


def gcd_recursive(a: int, b: int) -> int:
    """The same rule written as its own definition: gcd(a, b) = gcd(b, a mod b)."""
    if b == 0:
        return abs(a)
    return gcd_recursive(b, a % b)


print(gcd(1071, 462), gcd_recursive(1071, 462))
print(gcd(48, 18), gcd(18, 48))
print(gcd(0, 5), gcd(5, 0), gcd(0, 0))
print(gcd(-48, 18), gcd(270, 192))

Both versions are five lines of real work. The recursive one is the identity gcd(a, b) = gcd(b, a mod b) typed out verbatim, which makes it the better one to read; the iterative one is the better one to ship, since it uses no stack.

How the code maps to the idea

while b: is the termination test. In Python an integer is falsy only when it is 0, so this loop runs until the remainder reaches zero — exactly the stopping condition from the proof. When it exits, b is 0 and the answer sits in a.

a, b = b, a % b is the whole algorithm in one line. The right-hand side is fully evaluated before either assignment happens, so a % b still uses the old a. Split it into two statements and you destroy a before computing the remainder — that is the most common way this function gets broken.

The order of the arguments does not matter. If you call gcd(18, 48), the first iteration computes 18 mod 48, which is 18, so the pair becomes (48, 18) and the algorithm carries on normally. Passing them the wrong way round costs exactly one extra division, never a wrong answer.

abs on the way in handles negative inputs. Divisors do not care about sign — 6 divides both 48 and −48 — so the gcd is conventionally non-negative.

Zero needs no special case. gcd(5, 0) skips the loop and returns 5, which is right: every number divides 0, so the greatest divisor shared by 5 and 0 is 5. gcd(0, 5) takes one iteration to swap into the same state. gcd(0, 0) returns 0, the convention math.gcd also uses.

Now count the work instead of trusting the claim:

def count_subtraction_steps(a: int, b: int) -> int:
    """How many subtractions the naive form performs."""
    steps = 0
    while a != b:
        if a > b:
            a -= b
        else:
            b -= a
        steps += 1
    return steps


def count_division_steps(a: int, b: int) -> int:
    """How many divisions the modulo form performs."""
    steps = 0
    while b:
        a, b = b, a % b
        steps += 1
    return steps


print(f"{'pair':<20}{'subtractions':>14}{'divisions':>11}")
for first, second in [(48, 18), (1071, 462), (1_000_000, 2)]:
    pair = f"gcd({first}, {second})"
    print(f"{pair:<20}{count_subtraction_steps(first, second):>14}"
          f"{count_division_steps(first, second):>11}")

pair                  subtractions  divisions
gcd(48, 18)                      4          3
gcd(1071, 462)                  11          3
gcd(1000000, 2)             499999          1

That last row is the argument for division in one line: 499,999 steps against one.

Complexity

Time: O(log min(a, b)). Here is where the logarithm comes from.

Claim: when a is at least b, the remainder a mod b is smaller than a / 2. Two cases, and they cover everything. If b is at most a / 2, then the remainder is smaller than b, which is at most a / 2. If b is bigger than a / 2, then b fits into a exactly once, so the remainder is a - b, which is again smaller than a / 2. Either way the remainder loses at least half. (The loop guarantees the precondition: from the first division onwards, the new first number is the old second one, which is the larger of the pair.)

Now follow the pair through two steps. Starting from (a, b) you get (b, a mod b), then (a mod b, ...). So after two steps the first slot holds a value smaller than half of what it held before. A quantity that halves every two steps reaches 1 after about 2 * log2(a) steps, and the loop stops there. Within a division or two both numbers are at most the smaller input — straight away if you passed the larger one first, after one extra swap if you did not — so the bound is O(log min(a, b)) divisions.

Concretely: an eighteen-digit number is about 60 bits, so this argument caps the run at roughly 120 divisions — against 10¹⁸ trial divisions for the naive search. The true bound is sharper still, and it has a name.

The worst case is consecutive Fibonacci numbers, and this is not a curiosity — it is exactly the input that makes every quotient as small as possible. Each Fibonacci number is the sum of the two before it, so dividing one by its predecessor gives a quotient of 1 and a remainder equal to the one before that. Quotient 1 is the least progress a division step can make.

print(f"{'pair':<16}{'divisions':>10}")
previous, current = 1, 1
for _ in range(8):
    previous, current = current, previous + current
    pair = f"gcd({current}, {previous})"
    print(f"{pair:<16}{count_division_steps(current, previous):>10}")

worst_steps, worst_pair = 0, (0, 0)
for larger in range(2, 1000):
    for smaller in range(1, larger):
        steps = count_division_steps(larger, smaller)
        if steps > worst_steps:
            worst_steps, worst_pair = steps, (larger, smaller)
print(f"worst pair below 1000: gcd{worst_pair} takes {worst_steps} divisions")

pair             divisions
gcd(2, 1)                1
gcd(3, 2)                2
gcd(5, 3)                3
gcd(8, 5)                4
gcd(13, 8)               5
gcd(21, 13)              6
gcd(34, 21)              7
gcd(55, 34)              8
worst pair below 1000: gcd(987, 610) takes 14 divisions

Each rung of the Fibonacci ladder costs exactly one more division than the last, and the brute-force search over every pair below 1,000 lands on 987 and 610 — Fibonacci numbers again.

This is what Lamé's theorem (1844) formalises: if the algorithm needs n division steps on a pair whose smaller member is b, then b is at least the (n+1)-th Fibonacci number. Turned around, n is never more than five times the number of decimal digits of b — so an eighteen-digit pair takes at most 90 divisions, tightening the 120 the halving argument gave. The 14-step run above therefore needed a b of at least 610, which is exactly what the search found. Lamé's proof is generally counted as the first complexity analysis of an algorithm in history, written a century before there was a computer to run one.

The best case is O(1). When b already divides a, the first division leaves a remainder of 0 and the loop stops: gcd(48, 24) costs one division, and gcd(48, 0) costs none at all. Nothing can be cheaper than that, and it is a reachable case, not a theoretical one: any pair where one number is a multiple of the other lands on it.

The average case is logarithmic too, with a much smaller constant than the worst case. Count the divisions over random pairs and it works out at roughly two per decimal digit of the smaller number — about 35 divisions for an eighteen-digit pair, against the 90 Lamé permits. The halving argument still caps every single input, so no choice of inputs can push the average above O(log min(a, b)).

Space: O(1) for the iterative version — two integers and a temporary. The recursive version uses one stack frame per division, so O(log min(a, b)) frames; for any pair of 64-bit integers that is under 100 frames, comfortably inside Python's default recursion limit of 1,000.

One honest caveat: counting divisions assumes a division is one operation, which holds only while the numbers fit in a machine word. RSA keys are 2048 bits and up. Counting bit operations instead, the whole run costs about O(k²) for k-bit inputs with schoolbook division — practical, but the reason production libraries use the refinements below.

Two things you get for free

Least common multiple

The lowest common multiple of 4 and 6 is 12. There is no separate algorithm for it, because gcd and lcm are two halves of one fact:

a * b = gcd(a, b) * lcm(a, b)

The reason is visible in the prime factorisations. For each prime, the gcd takes the smaller exponent and the lcm takes the larger. The smaller plus the larger is the sum, and the sum is what the product a * b has. So the lcm is the product divided by the gcd.

def lcm(a: int, b: int) -> int:
    """Lowest common multiple, derived from the gcd.

    Divides before multiplying so the intermediate value never exceeds the
    answer, which matters in every language with fixed-width integers.
    """
    if a == 0 or b == 0:
        return 0
    return abs(a // gcd(a, b) * b)


print(lcm(4, 6), lcm(21, 6), lcm(1071, 462))
print(lcm(0, 7))

12 42 23562
0

Write it as a // gcd(a, b) * b, not a * b // gcd(a, b). The division is exact either way, but dividing first keeps the intermediate value no larger than the answer. In Python that is a habit; in C or Java it is the difference between a correct result and a wrapped one. The zero guard matters too — without it, lcm(0, 0) divides by a gcd of zero.

The extended Euclidean algorithm

Bézout's identity says something stronger than "the gcd exists": for any a and b there are integers x and y with

a * x + b * y = gcd(a, b)

Those coefficients are already hiding in the division trace — you just have to read it backwards. From the run on 1071 and 462:

21 = 462 - 3 * 147 (the second division, rearranged)
147 = 1071 - 2 * 462 (the first division, rearranged)

Substitute the second into the first: 21 = 462 - 3 * (1071 - 2 * 462) = 7 * 462 - 3 * 1071. So x = -3 and y = 7. Multiply it out if you like: −3,213 + 3,234 = 21.

The recursive implementation does that substitution automatically. Each call gets its child's coefficients and rewrites them in terms of its own two numbers.

def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return (g, x, y) such that a * x + b * y == g == gcd(a, b).

    Both inputs must be non-negative. There is no abs() on the way in here,
    unlike gcd above, because x and y have to match the numbers you passed;
    a negative input comes back with a negative g.

    The base case is free: gcd(a, 0) is a, and a * 1 + 0 * 0 == a. Every
    other level rewrites its child's answer in terms of its own two numbers.
    """
    if b == 0:
        return a, 1, 0
    g, x, y = extended_gcd(b, a % b)
    # The child solved b * x + (a mod b) * y == g. Substituting
    # a mod b == a - (a // b) * b and collecting terms gives the line below.
    return g, y, x - (a // b) * y


g, x, y = extended_gcd(1071, 462)
print(f"gcd = {g}, x = {x}, y = {y}")
print(f"1071 * ({x}) + 462 * {y} = {1071 * x + 462 * y}")

gcd = 21, x = -3, y = 7
1071 * (-3) + 462 * 7 = 21

The recursion descends exactly as far as the plain version does — same divisions, same bound, three extra arithmetic operations per level: the floor division, one multiplication and one subtraction. Extended Euclid is O(log min(a, b)) too, and because it is recursive it also carries the O(log min(a, b)) stack frames.

Modular inverses, and why RSA needs them

This is the payoff. In ordinary arithmetic the inverse of 17 is 1/17. Modular arithmetic has no fractions, so the inverse of a modulo m is the whole number x with a * x leaving a remainder of 1 when divided by m.

Such an x exists if and only if gcd(a, m) = 1, and the reason is Bézout. If the gcd is 1, then a * x + m * y = 1 for some x and y, and m * y vanishes modulo m — so a * x leaves remainder 1. If the gcd is bigger than 1, no multiple of a can ever land on 1, because everything in sight is divisible by that common factor.

So the extended algorithm is the modular inverse algorithm. Take x, reduce it into range, done.

def modular_inverse(value: int, modulus: int) -> int:
    """The number that multiplies `value` back to 1, working modulo `modulus`."""
    g, x, _ = extended_gcd(value % modulus, modulus)
    if g != 1:
        raise ValueError(f"{value} has no inverse modulo {modulus}")
    return x % modulus


public_exponent = 17
totient = 3120  # (61 - 1) * (53 - 1) for the textbook RSA key with n = 3233
private_exponent = modular_inverse(public_exponent, totient)
print(f"private exponent d = {private_exponent}")
print(f"17 * {private_exponent} mod 3120 = {(public_exponent * private_exponent) % totient}")

try:
    modular_inverse(6, 9)
except ValueError as error:
    print(error)

private exponent d = 2753
17 * 2753 mod 3120 = 1
6 has no inverse modulo 9

Those numbers are the textbook RSA key: primes 61 and 53, modulus 3,233, public exponent 17. The private exponent 2,753 is not chosen or searched for — it is the modular inverse of 17, produced by extended Euclid in five divisions. Generating an RSA key pair means picking a public exponent coprime to the totient (a gcd test) and then inverting it (extended Euclid). Both halves are this post.

When to use it, and when not to

Use it whenever you need a gcd, an lcm, or a modular inverse — in practice, fractions, ratios, cycle lengths and any modular arithmetic. It is one of the few algorithms with no real competition in its weight class.

But in Python, call the standard library. math.gcd is written in C, accepts any number of arguments since 3.9, and for large integers CPython switches to Lehmer's algorithm — a refinement that handles several quotients at a time from the leading digits alone, cutting the number of full-width divisions sharply. math.lcm arrived in 3.9, and since 3.8 pow(a, -1, m) gives you a modular inverse directly.

import math
from fractions import Fraction

print(math.gcd(1071, 462), math.lcm(4, 6, 10), math.gcd(12, 18, 24))
print(Fraction(1071, 462))
print(pow(17, -1, 3120))

for width, height in [(1920, 1080), (2560, 1440), (1366, 768)]:
    divisor = math.gcd(width, height)
    print(f"{width}x{height} -> {width // divisor}:{height // divisor}")

21 60 6
51/22
2753
1920x1080 -> 16:9
2560x1440 -> 16:9
1366x768 -> 683:384

Write it yourself when you need the Bézout coefficients — pow(a, -1, m) hands back the inverse but throws away x and y, which you need for solving equations of the form a * x + b * y = c. Or in a language with no gcd in its standard library. Or when you are learning, which is the best reason of all.

Do not use the subtraction form. Its only defence is hardware where division is unavailable or ruinously expensive, and even there the right answer is binary GCD (Stein's algorithm), which replaces division with subtraction, halving and parity tests.

Do not expect it to factorise anything. The gcd tells you what two numbers share, not what either one is made of. Finding the prime factors of a single number is a vastly harder problem — the one RSA's security rests on.

Where it shows up in the real world

Python's fractions.Fraction normalises every fraction on construction by dividing numerator and denominator by their gcd. That is why Fraction(1071, 462) printed 51/22 above. CPython's statistics module adds its partial sums as Fraction objects to stay exact, so gcds run on the way to a plain mean too.

RSA key generation. Choosing a public exponent means checking that it is coprime to (p - 1) * (q - 1), which is a gcd test, and the private exponent is that public exponent's modular inverse, which is extended Euclid. Both halves of generating a key pair are this algorithm.

pow(a, -1, m) in CPython is an extended-Euclid implementation. The same algorithm produces the inverse of q modulo p that a PKCS#1 RSA private key stores next to the two primes — the coefficient behind the Chinese remainder theorem shortcut that makes RSA decryption several times faster.

Aspect ratios. Dividing a screen's width and height by their gcd is how 1920×1080 becomes 16:9. It also exposes the odd ones: 1366×768 reduces no further than 683:384, because 683 is prime. A true 16:9 panel 768 pixels tall would be 1365.33 pixels wide, so that laptop resolution is famously almost 16:9 and never exactly.

Anything with repeating cycles. Two signals with periods 12 and 18 realign after lcm(12, 18) = 36 units. Gear ratios, polyrhythms in music, and traffic-light cycle planning are the same computation with different units.

Common mistakes

Passing zero to the subtraction version. With a = 0 and b = 5, the two are never equal, and b -= a subtracts nothing forever. The loop hangs. The division form has no such failure mode.

Returning the wrong variable. After while b: ends, b is 0 by definition. The answer is a. Returning b gives you a very confident zero.

Splitting the swap. Writing a = b followed by b = a % b computes the remainder of b divided by b, which is 0, so the function returns after one step with the wrong answer. Keep it as a single tuple assignment.

Porting the modulo naively. Python's % returns a non-negative result when the divisor is positive; C, Java, Go and JavaScript take the sign of the dividend instead, so -48 % 18 is −12 in those languages and 6 in Python. A direct port can return a negative gcd, and the subtraction form can loop forever. Take absolute values on entry.

Forgetting to reduce the Bézout coefficient. The x that comes back is often negative — for 17 modulo 3120 it is −367. A modular inverse must be reported in range, so return x % modulus.

Assuming an inverse always exists. It exists only when the gcd is 1. Check the returned g before using x, as modular_inverse does above; skipping the check gives you a number that quietly fails to invert anything.

Practice

Write gcd_of_list that returns the gcd of any number of integers, using functools.reduce over the two-argument version.
Reduce a fraction to lowest terms without using fractions, making sure a negative denominator moves its sign to the numerator.
Count the divisions for every pair below 10,000 and confirm that the worst pair is again two consecutive Fibonacci numbers.
Use extended_gcd to solve a * x + b * y = c for given a, b and c, reporting that no solution exists when the gcd does not divide c.
Implement binary GCD (Stein's algorithm) using only subtraction, halving and even/odd tests, and compare its step count against the division form on a few large pairs.

Summary

The Euclidean algorithm is the best return on five lines of code in all of computing. One observation — that a and b share their divisors with b and a mod b — turns a search over a quintillion candidates into fewer than a hundred divisions, and the extended version throws in Bézout coefficients and modular inverses at no extra asymptotic cost. Twenty-three centuries later, nothing has replaced it.


Difficulty	Easy
Best case	O(1) — `b` already divides `a`, so one division ends it
Average case	O(log min(a, b)) — the remainder at least halves every two steps
Worst case	O(log min(a, b)) — consecutive Fibonacci numbers, every quotient 1
Space	O(1) iterative; O(log min(a, b)) stack frames if recursive
Data structure	Two integers, nothing else
Handles negatives	Yes, with `abs` on entry; the gcd is taken as non-negative
Extended version	Same bound, returns `x` and `y` with `a * x + b * y = gcd(a, b)`
Use it when	You need a gcd, an lcm, a reduced fraction or a modular inverse
Avoid it when	Python has it already — and never use the subtraction form
Real-world use	RSA key generation, `Fraction` normalisation, aspect ratios, cycle lengths
Python equivalent	`math.gcd(a, b)`, `math.lcm(a, b)`, `pow(a, -1, m)`

Keep reading

Fast Exponentiation — the other half of modular arithmetic, and how RSA actually encrypts once it has these keys.
Sieve of Eratosthenes — the other ancient algorithm still in daily use, and the fastest way to list every prime below a limit.
Binary Search — the same logarithmic payoff from the same trick of discarding half the problem each step.
Big O Notation — the full treatment of the counting arguments used above.
Recursion and Backtracking — why the recursive form is safe here and dangerous elsewhere.

Rabin-Karp in Python: Finding Substrings With Rolling Hashes

Bimal Kshetri — Sun, 16 Aug 2026 16:37:34 +0000

Substring search done the obvious way re-reads the pattern at every position in the text. Rabin-Karp does something stranger: it turns every window of the text into a single number and compares numbers instead of strings. Comparing two integers costs the same whether the pattern is 3 characters or 300.

That trade only pays off because of one arithmetic trick. Once you know the number for the window starting at position i, you can get the number for the window starting at i + 1 in constant time — subtract the leading character's contribution, shift, add the new trailing character. That is the rolling hash, and it is the entire algorithm.

Be clear about where it is genuinely useful, though. For finding one string inside another in Python, text.find(pattern) beats anything you write here, and KMP has the better worst case. Rabin-Karp earns its place elsewhere: hunting ten thousand patterns at once, or wanting a cheap hash of every window for its own sake — deduplication, backup chunking, plagiarism detection.

The idea

Give every string of length m a number, by reading it as a number written in base B — exactly as 407 is 4 × 10² + 0 × 10 + 7 in base 10. Map the characters to digits (a is 0, b is 1, up to z as 25) and take the base to be the alphabet size, 26. Then "ban" is 1 × 26² + 0 × 26 + 13 = 689.

That mapping is exact: distinct strings of the same length get distinct numbers, because base-26 digits are unique. But you cannot afford the numbers — a 40-character window in base 256 is a 320-bit integer, and arithmetic on it is no longer constant time. So take it modulo a prime P that fits in a machine word. The hash is now small and fast, and equal hashes only mean the strings are probably equal. That word "probably" separates a correct implementation from a broken one.

Here is the baseline you are trying to beat:

def naive_search(text: str, pattern: str) -> list[int]:
    """Return every start index where pattern occurs in text, by direct comparison."""
    n, m = len(text), len(pattern)
    matches: list[int] = []
    for start in range(n - m + 1):
        # The slice copies m characters before the comparison even starts.
        if text[start:start + m] == pattern:
            matches.append(start)
    return matches


print(naive_search("banana", "ana"))
print(naive_search("aaaaa", "aa"))

[1, 3]
[0, 1, 2, 3]

There are n - m + 1 starting positions and each comparison can cost up to m characters, so the worst case is O(n·m). A character-by-character version behaves far better on English text, because most comparisons die on the first character — but "usually fine" is not a bound, and the pathological input ("aaaa…a" searched for "aaa…ab") really does cost the full product. The slice above never gets even the good case: text[start:start + m] copies m characters before the comparison begins, so it pays the full product on every input.

Rabin-Karp replaces that inner comparison with an integer comparison:

The rolling update

Computing each of those hashes from scratch costs O(m), and there are n - m + 1 of them, so you are back at O(n·m) and have gained nothing. The win comes from how much consecutive windows share.

Window i covers text[i] through text[i + m - 1]; window i + 1 covers text[i + 1] through text[i + m]. They share m - 1 characters, so going from one to the next is three digit operations:

Remove the leading digit. text[i] is the most significant digit of window i, contributing value(text[i]) · B^(m-1). Subtract exactly that.
Shift left one place. Multiply by B, moving every remaining character up one power. That is what sliding means.
Append the new digit. Add value(text[i + m]) in the now-empty units place.

As one line, everything taken modulo P:

hash(i + 1) = ((hash(i) - value(text[i]) * high_power) * B + value(text[i + m])) % P

where high_power is B^(m-1) % P, computed once before the loop. Two multiplications, a subtraction, an addition and one modulo — and not one of them depends on m. That is the O(1) step that makes the whole thing worth doing.

Watching it work

Take the text "banana" and the pattern "ana", with base 26 and modulus 101. The letter values are a = 0, b = 1, n = 13.

The pattern first: "ana" is 0 × 676 + 13 × 26 + 0 = 338, and 338 mod 101 = 35. That number is computed once and never again. The first window, "ban", is 1 × 676 + 0 × 26 + 13 = 689, and 689 mod 101 = 83. Not 35, so no match at position 0. Every window after that comes from the rolling update, with high_power = 26² mod 101 = 676 mod 101 = 70.

Roll from "ban" to "ana". The character leaving is b, worth 1, so subtract 1 × 70 = 70: 83 - 70 = 13. Shift: 13 × 26 = 338. The character entering is a, worth 0, so add nothing. 338 mod 101 = 35, which equals the pattern hash — so verify the characters, text[1:4] really is "ana", and report a match at index 1.

Roll from "ana" to "nan". Leaving is a, worth 0, so subtract nothing: still 35. Shift: 35 × 26 = 910. Entering is n, worth 13: 910 + 13 = 923, and 923 mod 101 = 14 since 101 × 9 = 909. Not 35, no match.

Roll from "nan" to "ana". Leaving is n, worth 13: 14 - 13 × 70 = 14 - 910 = -896. This is where hand-written implementations break, because the intermediate went negative. Python's % always returns a non-negative result for a positive modulus, so -896 % 101 is 13 and the arithmetic carries on. Shift and append: 13 × 26 = 338, plus 0, mod 101 = 35. Another hit, verified, match at index 3.

Here are those hashes computed the slow way, one full pass per window:

BASE = 26
MOD = 101


def letter_value(character: str) -> int:
    """Map 'a'..'z' onto the digits 0..25 of a base-26 number."""
    return ord(character) - ord("a")


def slow_window_hash(window: str) -> int:
    """Read a string as a base-26 number, kept small by a modulus. Costs O(len(window))."""
    value = 0
    for character in window:
        value = (value * BASE + letter_value(character)) % MOD
    return value


text, pattern = "banana", "ana"
m = len(pattern)

print("pattern", pattern, "->", slow_window_hash(pattern))
for start in range(len(text) - m + 1):
    print(start, text[start:start + m], "->", slow_window_hash(text[start:start + m]))

pattern ana -> 35
0 ban -> 83
1 ana -> 35
2 nan -> 14
3 ana -> 35

And here is the same sequence of numbers produced by rolling, at a fixed cost per window instead of a fresh pass over each one:

high_power = pow(BASE, m - 1, MOD)  # 26 ** 2 mod 101: the weight of the leading character
rolling = slow_window_hash(text[:m])

print("high_power", high_power)
print(text[:m], rolling)
for start in range(len(text) - m):
    leaving = letter_value(text[start])
    entering = letter_value(text[start + m])
    rolling = ((rolling - leaving * high_power) * BASE + entering) % MOD
    print(text[start + 1:start + 1 + m], rolling)

high_power 70
ban 83
ana 35
nan 14
ana 35

Identical hashes, one constant-time step each. Note pow(BASE, m - 1, MOD): Python's three-argument pow does modular exponentiation in O(log m) multiplications rather than m of them, the technique in fast exponentiation.

Spurious hits

Now break it deliberately. Keep base 26 but drop the modulus to 31, still prime but far too small. "ana" is 338 and 338 mod 31 = 28; "nan" is 13 × 676 + 0 + 13 = 8801 and 8801 mod 31 = 28 as well. Two windows sharing no arrangement of characters now hash identically.

That is a spurious hit: equal hashes, unequal strings. You cannot eliminate it, because you are mapping every possible m-character string into P buckets and there are vastly more strings than buckets. The only defence is the one the algorithm builds in: when the hashes match, compare the actual characters before reporting anything. Skip that check to save time and you have written a search that returns wrong answers on inputs you will never think to test.

A large prime does not remove collisions, it makes them rare — which is the whole point of the modulus choice:

Prime matters. A modulus sharing factors with the base throws information away. The extreme case is a power of the base, B = 256 with P = 2^24: that keeps only the window's last three bytes and discards every earlier character, so "...abc" and "...xyz abc" collide by construction. A prime shares no factors with a sensible base, so all P residues stay reachable.
Large matters. For a hash spread evenly over P values, two different windows collide with probability about 1/P. With P near 10^9, scanning a million-character text you expect roughly 10^6 / 10^9 = 0.001 spurious hits. With P = 31 you expect one every 31 windows.
The base should be at least the alphabet size, so distinct characters are distinct digits, and coprime to the modulus. 256 for bytes; 31 and 257 are common for text.

The code

The real implementation, with a byte-sized base and a prime near a billion:

def rabin_karp_search(text: str, pattern: str, base: int = 256,
                      modulus: int = 1_000_000_007) -> list[int]:
    """Return every start index where pattern occurs in text.

    Each window's hash is derived from the previous window's in constant time,
    and every hash match is verified against the real characters before it is
    reported, because two different windows can share a hash.
    """
    n, m = len(text), len(pattern)
    if m == 0:
        return list(range(n + 1))
    if m > n:
        return []

    high_power = pow(base, m - 1, modulus)
    pattern_hash = 0
    window = 0
    for index in range(m):
        pattern_hash = (pattern_hash * base + ord(pattern[index])) % modulus
        window = (window * base + ord(text[index])) % modulus

    matches: list[int] = []
    for start in range(n - m + 1):
        # The slice runs only when the hashes agree, so it is rare on real text.
        if window == pattern_hash and text[start:start + m] == pattern:
            matches.append(start)
        if start < n - m:
            leaving = ord(text[start])
            entering = ord(text[start + m])
            window = ((window - leaving * high_power) * base + entering) % modulus

    return matches


print(rabin_karp_search("banana", "ana"))
print(rabin_karp_search("abracadabra", "abra"))
print(rabin_karp_search("aaaa", "aa"))
print(rabin_karp_search("hello", "world"))
print(rabin_karp_search("hi", "longer pattern"))

[1, 3]
[0, 7]
[0, 1, 2]
[]
[]

How the code maps to the idea

The setup loop builds both hashes at once, walking the pattern and the first window together with Horner's method: multiply the running value by the base, add the next digit, reduce. That loop is the O(m) preprocessing cost, and the only place a hash is built from scratch.

high_power is B^(m-1) mod P — the weight of the leading character. It never changes, so it is computed once. Getting it wrong by one power (using B^m) is the most common bug here. The first window is still hashed from scratch, so a match at index 0 survives; every rolled hash after it is wrong, so the search quietly returns a short list rather than raising.

ord(character) is the digit mapping. Base 256 with ord handles ASCII directly. Non-ASCII survives too — ord goes well past 255, so digits can exceed the base, but the hash stays a deterministic function of the characters. For byte-level work, encode to UTF-8 first and hash the bytes.

The verification is the and in the condition. Python short-circuits, so the slice comparison only runs when the hashes already agree. On text with few matches it almost never executes, which is exactly why the average case is linear.

The guarded roll. The update sits behind if start < n - m because the last window has no text[start + m] to read. Without the guard, every search ends in an IndexError.

Edge cases. An empty pattern returns every position from 0 to n, matching the convention that an empty string occurs everywhere; it needs a special case because high_power is meaningless for m = 0. A pattern longer than the text returns before any arithmetic. Overlapping matches fall out free — the window advances one character, not m, which is why "aaaa" searched for "aa" reports 0, 1 and 2.

Spurious hits are measurable. This counts how often the hashes agree against how often the strings really match, on a 1,200-character text built by repeating "banana":

def count_hash_hits(text: str, pattern: str, modulus: int) -> tuple[int, int]:
    """Return (hash hits, verified matches) for the base-26 lowercase hash."""
    m = len(pattern)
    high_power = pow(BASE, m - 1, modulus)
    pattern_hash = window = 0
    for index in range(m):
        pattern_hash = (pattern_hash * BASE + letter_value(pattern[index])) % modulus
        window = (window * BASE + letter_value(text[index])) % modulus

    hits = verified = 0
    for start in range(len(text) - m + 1):
        if window == pattern_hash:
            hits += 1
            if text[start:start + m] == pattern:
                verified += 1
        if start < len(text) - m:
            window = ((window - letter_value(text[start]) * high_power) * BASE
                      + letter_value(text[start + m])) % modulus
    return hits, verified


haystack = "banana" * 200  # 1,200 characters holding 400 real occurrences of "ana"
for modulus in (31, 101, 1_000_000_007):
    hits, verified = count_hash_hits(haystack, "ana", modulus)
    print(f"modulus {modulus:>10}: {hits} hash hits, {verified} verified matches")

modulus         31: 600 hash hits, 400 verified matches
modulus        101: 400 hash hits, 400 verified matches
modulus 1000000007: 400 hash hits, 400 verified matches

Modulus 31 produces 200 wasted verifications — one for every "nan" in the text. Modulus 101 already separates those two windows. The billion-sized prime is what you would actually ship.

Complexity

Preprocessing: O(m). One pass over the pattern and one over the first window, plus pow(base, m - 1, modulus), which is O(log m) multiplications.

The scan: O(n) hash work. There are n - m + 1 windows, and each rolling update is a fixed two multiplications, one subtraction, one addition and one modulo. No part of it depends on m.

Verification: O(m) per hash hit. Two kinds of hit exist. Real matches: if the pattern occurs k times you compare k · m characters. That cost is Rabin-Karp's own, not a law of string search — KMP reports the same k occurrences in O(n + m) total however large k is, because it never re-reads a character it has already matched. Spurious hits: with the hash spread evenly over P residues each window collides with probability about 1/P, so the expected count is n / P and the expected wasted work is n · m / P characters. On a million-character text with a 100-character pattern and P = 10^9, that is 10^6 × 100 / 10^9 — a tenth of one character comparison.

Best case: O(n + m). There is no early exit, so the loop always performs all n - m + 1 rolling updates. The floor is therefore the O(m) preprocessing plus one constant-time step per window, and it is reached exactly when no window's hash ever equals the pattern's — no verification runs at all.

Average case: O(n + m). Add the three: O(m) preprocessing, O(n) rolling, O(k · m + n · m / P) verification. With few occurrences and a large prime the last term rounds away, leaving a bound linear in the total input size.

Worst case: O(n·m). This happens whenever every window produces a hit. The easy trigger is a text and pattern of one repeated character, where every window is a genuine match and each verification costs the full m:

worst_text = "a" * 2000
worst_pattern = "a" * 100
found = rabin_karp_search(worst_text, worst_pattern)
print(len(found), "windows matched;", len(found) * len(worst_pattern), "characters compared")

1901 windows matched; 190100 characters compared

KMP finds those same 1,901 matches in about n + m = 2,100 steps, which is the concrete version of "KMP has the better worst case".

The dangerous trigger is adversarial. Base and modulus are usually constants in the source, so anyone who can read the code can build a text whose every window collides with the pattern hash, turning a linear search quadratic on demand. If the input comes from strangers, pick the modulus or base randomly at process start.

Space: O(1) working memory. A handful of integers, however long the text is. Two things sit outside that count: the returned list of positions, which is output rather than working space, and the text[start:start + m] slice, which Python materialises as a throwaway m characters wide each time a hash agrees. The multi-pattern version below needs O(k) for a dictionary of k pattern hashes.

One caveat about "constant time" arithmetic: it holds while the numbers fit in a machine word. Keep the modulus below 2^31 if intermediate products must fit in a signed 64-bit integer, since (window - leaving * high_power) * base can reach roughly P × B² — about 2^47 when P is near 2^31 and B is 256. Python integers are arbitrary precision so nothing overflows, but arithmetic on numbers far larger than a word stops being constant time.

Searching for many patterns at once

This is where Rabin-Karp stops being an academic exercise. A window's hash does not know what it is being compared against, so instead of comparing it to one number, look it up in a dictionary of many.

Put every pattern's hash in a dict keyed by hash and run the same single pass, checking each window's hash against the dict. That lookup is O(1) whether you are searching for 4 patterns or 40,000. KMP cannot do this — its skip table is built from one specific pattern, so k patterns means k passes.

def rabin_karp_multi(text: str, patterns: list[str], base: int = 256,
                     modulus: int = 1_000_000_007) -> list[tuple[int, str]]:
    """Find every occurrence of any of the patterns, which must share one length.

    The rolling hash is unchanged; only the comparison changes, from one number
    to a dictionary lookup. That lookup costs the same whether there is one
    pattern or a hundred thousand.
    """
    if not patterns:
        return []
    m = len(patterns[0])
    if any(len(item) != m for item in patterns):
        raise ValueError("every pattern must have the same length")

    n = len(text)
    if m == 0 or m > n:
        return []

    by_hash: dict[int, list[str]] = {}
    for item in patterns:
        digest = 0
        for character in item:
            digest = (digest * base + ord(character)) % modulus
        by_hash.setdefault(digest, []).append(item)

    high_power = pow(base, m - 1, modulus)
    window = 0
    for index in range(m):
        window = (window * base + ord(text[index])) % modulus

    found: list[tuple[int, str]] = []
    for start in range(n - m + 1):
        # Two different patterns can collide too, so still verify each candidate.
        for candidate in by_hash.get(window, ()):
            if text[start:start + m] == candidate:
                found.append((start, candidate))
        if start < n - m:
            window = ((window - ord(text[start]) * high_power) * base
                      + ord(text[start + m])) % modulus

    return found


sentence = "the quick brown fox jumps over the lazy dog"
print(rabin_karp_multi(sentence, ["quick", "brown", "jumps", "zebra"]))

needles = [f"id{index:03d}" for index in range(1000)]  # 1,000 patterns, all length 5
log_line = "user id042 opened id777 and closed id042 again"
print(rabin_karp_multi(log_line, needles))

[(4, 'quick'), (10, 'brown'), (20, 'jumps')]
[(5, 'id042'), (18, 'id777'), (35, 'id042')]

The second search checks a thousand patterns against 46 characters in one pass, at the same scanning cost as searching for one — only the dictionary build is proportional to the number of patterns. The dict values are lists rather than single strings because two patterns can collide with each other, which needs the same verification treatment.

The restriction is real, though: every pattern must be the same length, because one rolling hash tracks one window size. Patterns of d distinct lengths need d passes. If d is large, use Aho-Corasick instead — it handles mixed lengths in one pass with an automaton built from a trie.

When to use it, and when not to

Use it for many patterns of one length. Blocklists of fixed-length tokens, a set of known hashes, a dictionary of 5-grams: one pass, one dict lookup per position.

Use it when the rolling hash is the product, not the search. Content-defined chunking, fingerprinting every substring of a document, finding duplicate blocks between two files — these want a hash of every window and never do exact matching at all. The rolling update is the only reason they are affordable.

Use it when you want a short, obviously-correct implementation. It is about 20 lines you can re-derive from the base-B idea. KMP's failure function is easier to get subtly wrong.

Do not use it for a single pattern in Python. The built-ins are written in C and are strictly better:

import re

print("banana".find("ana"))
print([match.start() for match in re.finditer("(?=ana)", "banana")])

1
[1, 3]

str.find and the in operator use CPython's fastsearch: a Boyer-Moore-Horspool variant with a bitmask skip table, plus the two-way algorithm for long needles since Python 3.10 to avoid quadratic blowups. The lookahead in re.finditer("(?=ana)", …) is how you get overlapping matches, which a find loop misses unless it advances by one character rather than by the pattern length.

Do not use it when the worst case must be bounded. Rabin-Karp is O(n·m) in the worst case, and an attacker can construct that case. KMP is O(n + m) always, with no probabilistic argument anywhere in it.

Do not use it for approximate matching. A rolling hash tells you two windows are equal, never how similar they are — change one character and the hash is unrelated. For "how different are these two strings", reach for edit distance.

Where it shows up in the real world

rsync is the clearest example of the shape, even though it does not use a polynomial hash. To work out which blocks of a file the far end already has, it rolls a weak checksum (a variant of Adler-32) over every byte offset, looks each value up in a hash table of the remote block checksums, and only then confirms a hit with a strong hash — MD4 originally, MD5 in current versions. Cheap rolling hash, expensive verification.

Content-defined chunking is Rabin's fingerprint doing what it was designed for. The Low-Bandwidth Network File System (LBFS, 2001) introduced it: roll a hash over a small sliding window and cut a chunk boundary wherever the hash's low bits are all zero. Because boundaries depend on content rather than offset, inserting a byte at the front of a file shifts only the chunk containing it — every other chunk keeps its identity and does not need re-uploading. restic does exactly this with a Rabin fingerprint in its chunker package; borgbackup does the same job with a different rolling hash (buzhash), which tells you the structure matters more than the specific function.

Plagiarism and near-duplicate detection. The winnowing algorithm behind MOSS (Schleimer, Wilkerson and Aiken, 2003) hashes every k-gram of a document and keeps a sample as its fingerprint; Karp-Rabin hashing is what makes hashing every k-gram affordable. Broder's shingling work on near-duplicate web pages at AltaVista used Rabin fingerprints of overlapping word sequences for the same reason.

Interview and contest problems. Rolling hashes solve a family of questions that look unrelated to string search: counting distinct substrings of a given length, finding the longest substring common to two strings by binary searching the length, and comparing two long substrings in O(1) after O(n) preprocessing.

Common mistakes

Not verifying the characters. The serious one. Equal hashes are evidence, not proof. An implementation that reports a match on hash equality alone is wrong, and wrong rarely enough that testing will not find it.

Recomputing each window hash from scratch. If the inner loop walks m characters you have written the naive algorithm with extra arithmetic. The update must be O(1).

Using B^m instead of B^(m-1). The leading character sits in place m - 1. Every rolled hash after the first window is then wrong, so the search loses every match except one that starts at index 0 — and loses them silently. Test on a text where the pattern occurs late, not only at the start.

Letting the intermediate go negative in a language without floor modulo. window - leaving * high_power is often negative. Python's % returns a non-negative result, so this just works; in C, C++, Java, Go or Rust it does not, and the search silently misses matches unless you add modulus back.

A modulus that is too small, or shares factors with the base. Modulus 31 gave a spurious hit every few windows above, and a power-of-two modulus with base 256 is worse than small — it discards everything but the last few characters.

Reaching for Python's built-in hash(). It cannot roll, and it is randomised per process by default, so results are not reproducible across runs.

Forgetting overlapping matches. After a match at index i the next window starts at i + 1, not i + m. Skipping ahead by the pattern length misses the second "ana" in "banana".

Practice

Add an early return to rabin_karp_search so it stops at the first match and returns -1 when there is none, mirroring str.find.
Count the distinct substrings of length k in a string by rolling one hash across it and collecting the values in a set, then explain why the answer can be slightly too low.
Extend rabin_karp_multi to accept patterns of different lengths by grouping them by length and running one pass per group, and count how many passes a realistic blocklist needs.
Implement double hashing: carry two rolling hashes with different primes and treat a hit as a hit only when both agree. Measure how many spurious hits survive on the "banana" * 200 text with moduli 31 and 37.
Write a content-defined chunker: roll a hash over a 48-byte window and cut a boundary wherever hash % 4096 == 0. Insert a character at the front of the input and show that all chunks after the first are unchanged.

Summary

Rabin-Karp is a good algorithm attached to a great primitive. The search itself is rarely the right tool for one pattern in one text — the standard library is faster and KMP has a stronger guarantee. The rolling hash underneath it is what matters: constant-time hashing of every window is the foundation of backup deduplication, rsync's block matching, and every plagiarism detector that fingerprints k-grams. Learn the update arithmetic, and remember that a hash match is a hint that must always be checked.


Difficulty	Medium
Preprocessing	O(m) — hash the pattern and the first window
Best case	O(n + m) — no window's hash ever collides
Average case	O(n + m) — one O(1) roll per window, about n/P spurious hits
Worst case	O(n·m) — every window hits, so every window is verified
Space	O(1) working memory for one pattern, O(k) for k patterns
Finds overlapping matches	Yes — the window advances one character at a time
Multi-pattern	Yes — any number of patterns of one length, in a single pass
Deterministic	Yes on output, probabilistic only on running time
Data structure	String or byte sequence, plus a dict of pattern hashes
Use it when	Many same-length patterns, or you need a hash of every window
Avoid it when	One pattern and a bounded worst case matters — use KMP
Real-world use	rsync block matching, restic chunking, MOSS plagiarism fingerprints
Python equivalent	`text.find(pattern)` / `re.finditer` — CPython's C `fastsearch`

Keep reading

The KMP Algorithm — the other classic string search, with a guaranteed linear worst case and no hashing at all.
Hash Tables — where hash functions, collisions and buckets are explained from first principles.
The Sliding Window Technique — the same add-one-remove-one idea applied to sums and counts instead of hashes.
Fast Exponentiation — how pow(base, m - 1, modulus) computes a huge modular power in O(log m) steps.
Big O Notation — the counting arguments used to justify O(n + m) above.

The KMP Algorithm in Python: String Search Without Ever Going Backwards

Bimal Kshetri — Sun, 16 Aug 2026 16:37:03 +0000

Searching a text for a substring looks like it should be free. Line the pattern up at the start, compare characters until one disagrees, slide one place right, try again. Python spells it "needle" in haystack and you never think about it.

The obvious implementation has a failure mode that is easy to trigger. Looking for aaaaab inside ten a characters followed by a b — an eleven character text — costs 36 character comparisons. Scale both sides up and the cost is the product of the two lengths, not the sum.

The waste has a precise shape. Every time an alignment fails, the naive search slides the pattern one place right and starts from scratch, even though it just read five characters that matched and is about to read four of them again. It learned something and threw it away. Knuth-Morris-Pratt, published in 1977 by Donald Knuth, James Morris and Vaughan Pratt, keeps it: one small table computed from the pattern alone, and with it the pointer into the text never moves backwards, not once, on any input. The search costs O(n + m), the text length plus the pattern length, with no bad case hiding behind the average.

The idea

The naive version, and where it wastes work

Start with the version you would write without thinking: at every possible starting position, compare the pattern character by character and give up at the first disagreement.

def naive_search(text: str, pattern: str) -> list[int]:
    """Return every start index where pattern occurs in text, by brute force.

    For each possible alignment of the pattern against the text, compare
    characters left to right and abandon the alignment at the first mismatch.
    """
    n, m = len(text), len(pattern)
    if m == 0:
        return list(range(n + 1))

    matches: list[int] = []
    for start in range(n - m + 1):
        offset = 0
        while offset < m and text[start + offset] == pattern[offset]:
            offset += 1
        if offset == m:
            matches.append(start)
    return matches


print(naive_search("ABABDABABABCABAB", "ABABCABAB"))
print(naive_search("aaaaaaaaaab", "aaaaab"))
print(naive_search("banana", "ana"))

[7]
[5]
[1, 3]

That is correct, including the overlapping matches of ana in banana at 1 and 3. The cost only bites when the pattern nearly matches in many places, so add a counter and measure.

def naive_search_counted(text: str, pattern: str) -> tuple[list[int], int]:
    """Brute force again, but also report how many characters it compared."""
    n, m = len(text), len(pattern)
    matches: list[int] = []
    comparisons = 0

    for start in range(n - m + 1):
        offset = 0
        while offset < m:
            comparisons += 1
            if text[start + offset] != pattern[offset]:
                break
            offset += 1
        if offset == m:
            matches.append(start)
    return matches, comparisons


for sample_text, sample_pattern in [
    ("ABABDABABABCABAB", "ABABCABAB"),
    ("aaaaaaaaaab", "aaaaab"),
]:
    found, count = naive_search_counted(sample_text, sample_pattern)
    print(f"{sample_pattern} in {sample_text}: {found}, {count} comparisons")

ABABCABAB in ABABDABABABCABAB: [7], 26 comparisons
aaaaab in aaaaaaaaaab: [5], 36 comparisons

Thirty-six comparisons to search eleven characters. There are n - m + 1 alignments and each can run the full length of the pattern before failing, so the worst case is (n - m + 1) * m comparisons: O(n × m). On a 1 MB text with a 1,000 character pattern of that shape, a billion comparisons.

Repetitive input is not exotic: DNA has a four-letter alphabet, log lines repeat the same prefixes, binary formats are full of zero bytes. And if an attacker picks the input, they can make this quadratic on purpose.

What a partial match already tells you

Here is the observation the whole algorithm rests on. Suppose you aligned the pattern at text position s, matched j characters, then hit a mismatch. You now know something exact and free: text[s : s + j] equals pattern[: j].

The naive version deletes that. But you can work out in advance, without looking at the text at all, which shifts are worth trying. A match beginning at s + k, for a shift k between 1 and j, needs its first j - k characters to line up with text[s + k : s + j] — which you already know, because those characters are pattern[k : j]. So the shift is worth trying only if

pattern[: j - k] equals pattern[k : j]

which says: some prefix of the pattern is also a suffix of pattern[: j]. A string that is both a proper prefix and a suffix of another string is called a border. "Proper" means not the whole string — otherwise everything borders itself and nothing ever shifts.

To shift as little as possible, so you never jump over a match, you want the longest border. The longest border of ABABCABAB is ABAB: those four characters open the pattern and also close it.

So after matching all nine characters, the next alignment to try is not one place right and not nine — it is five places right, because four characters are already verified. The longest border of every prefix of the pattern is the entire table KMP needs. It is called the failure function, or the LPS array, for longest proper prefix which is also a suffix.

Watching it work

Building the table for ABABCABAB

lps[index] is the length of the longest border of pattern[: index + 1]. Build it left to right carrying one variable, length, the border length of the previous prefix.

Index 0, A. A one-character string has no proper prefix, so lps[0] = 0. Always.
Index 1, B. length is 0, so try to extend the empty border: is pattern[1] equal to pattern[0]? B against A, no. lps[1] = 0.
Index 2, A. length is 0. Is pattern[2] equal to pattern[0]? A against A, yes. length becomes 1, lps[2] = 1. The border of ABA is A.
Index 3, B. length is 1. Is pattern[3] equal to pattern[1]? Yes. length becomes 2, lps[3] = 2. The border of ABAB is AB.
Index 4, C. length is 2. Is pattern[4] equal to pattern[2]? C against A, no — AB cannot be extended. Fall back to the border of AB, which is lps[1] = 0, and try again: C against pattern[0], still no. Nothing is left, so lps[4] = 0. The C occurs nowhere else in the pattern, so it destroys every border.
Index 5, A. length is 0 and A matches pattern[0], so lps[5] = 1.
Index 6, B. Matches pattern[1]. lps[6] = 2.
Index 7, A. Matches pattern[2]. lps[7] = 3.
Index 8, B. Matches pattern[3]. lps[8] = 4 — the ABAB from the diagram above.

Index 4 is the step worth staring at. When a border cannot be extended, the next candidate is the border of that border. That is why the table can be built out of itself, and it is the same move the search makes on a mismatch.

Searching ABABDABABABCABAB

Now search the text ABABDABABABCABAB with i as the text pointer and j as the pattern pointer. j doubles as "how many characters currently match", so the pattern sits at text position i - j.

Alignment at 0. A, B, A, B match, so i and j both reach 4. Then text[4] is D and pattern[4] is C. Mismatch, four characters matched. Consult lps[3] = 2: two of them survive. Set j = 2 and leave i at 4. The pattern has slid right by two.

Alignment at 2. Compare text[4], still D, against pattern[2], A. Mismatch again with j = 2, so consult lps[1] = 0. Set j = 0, i unchanged.

Alignment at 4. With j at 0 there is nothing left to fall back on. D against A fails, and only now does i move, to 5.

Alignment at 5. A, B, A, B match again, taking i to 9 and j to 4. text[9] is A, pattern[4] is C. Mismatch, lps[3] = 2 again, j = 2, i stays at 9.

Alignment at 7. The payoff. Compare text[9] against pattern[2]: A against A, a match. The characters text[7:9] were never re-read, and the search rolls straight on to a full match starting at index 7.

Five alignments, 19 character comparisons, and the text pointer moved right 16 times and left zero times. Naive search needed 26 comparisons on the same input.

The code

Two functions. The table is the harder one.

def build_lps(pattern: str) -> list[int]:
    """Longest proper prefix that is also a suffix, for every prefix of pattern.

    lps[index] is the length of the longest string that is both a proper
    prefix and a suffix of pattern[:index + 1]. "Proper" means it is not the
    whole prefix, so the value is always at most index.
    """
    lps = [0] * len(pattern)
    length = 0  # length of the border of the prefix ending at index - 1

    for index in range(1, len(pattern)):
        # The border cannot be extended by pattern[index], so try the next
        # shortest border. lps[length - 1] is exactly that border's length.
        while length > 0 and pattern[index] != pattern[length]:
            length = lps[length - 1]

        if pattern[index] == pattern[length]:
            length += 1

        lps[index] = length

    return lps


demo_pattern = "ABABCABAB"
print(" ".join(demo_pattern))
print(" ".join(str(value) for value in build_lps(demo_pattern)))
print(build_lps("aaaaab"))
print(build_lps("abcdef"))

A B A B C A B A B
0 0 1 2 0 1 2 3 4
[0, 1, 2, 3, 4, 0]
[0, 0, 0, 0, 0, 0]

Those are the two extremes and one real case. A pattern of distinct characters has no borders, so every entry is 0 and KMP degenerates to naive search — which costs nothing, because naive search is already linear when nothing ever partially matches. A run of identical characters has the maximum possible borders, so the table counts up.

The search loop is short once the table exists.

def kmp_search(text: str, pattern: str) -> list[int]:
    """Every start index of pattern in text, found in O(len(text) + len(pattern)).

    The text index never decreases: a mismatch rewinds the pattern index to a
    shorter border instead of re-reading characters of the text.
    """
    if not pattern:
        return list(range(len(text) + 1))

    lps = build_lps(pattern)
    matches: list[int] = []
    text_index = 0
    pattern_index = 0

    while text_index < len(text):
        if text[text_index] == pattern[pattern_index]:
            text_index += 1
            pattern_index += 1
            if pattern_index == len(pattern):
                matches.append(text_index - pattern_index)
                # Overlapping matches count, so slide to the longest border
                # rather than starting the pattern from scratch.
                pattern_index = lps[pattern_index - 1]
        elif pattern_index > 0:
            pattern_index = lps[pattern_index - 1]  # text_index stays put
        else:
            text_index += 1

    return matches


print(kmp_search("ABABDABABABCABAB", "ABABCABAB"))
print(kmp_search("banana", "ana"))
print(kmp_search("aaaa", "aa"))
print(kmp_search("hello", "xyz"))

[7]
[1, 3]
[0, 1, 2]
[]

How the code maps to the idea

lps[length - 1] is the fallback, in both functions. If the current border has length length, its own longest border has length lps[length - 1]. Following that chain from a border to the border of the border enumerates every border of the current prefix, longest first — which is the complete list of shifts worth trying, in the order that shifts least.

The off-by-one is j - 1, not j. On a mismatch you have matched pattern[: j], a string of length j, whose table entry lives at index j - 1. Writing lps[pattern_index] is the most common KMP bug: it consults the border of a prefix one character longer than the one you actually matched, so the pattern shifts by the wrong amount in whichever direction that entry happens to point. On ABABCABAB it shifts too far and the match at index 7 disappears; on aa the entry equals the current j, so nothing shrinks and the loop spins forever.

The three branches cover every case. Characters agree: advance both pointers. Characters disagree with something matched: rewind j, hold i. Characters disagree with nothing matched: there is no partial match to salvage, so advance i. Note that i increases in two branches and decreases in none — the whole "never goes backwards" claim is visible in three lines.

After a full match, j becomes lps[m - 1] rather than 0. Resetting to 0 works for non-overlapping searches but misses aa at index 1 of aaaa. Using the border treats a completed match exactly like a mismatch one character past the end, which is what it is.

Edge cases fall out of the arithmetic. An empty text never enters the loop. A pattern longer than the text can never drive j to len(pattern). A pattern that never matches returns an empty list with no special casing. The one case needing a guard is the empty pattern, since pattern[0] raises IndexError; the early return matches Python's convention, where "abc".find("") is 0.

Here is every comparison the search makes, printed by the algorithm itself, so you can check the hand trace above line by line.

def kmp_trace(text: str, pattern: str) -> None:
    """Print one line per character comparison the search performs."""
    lps = build_lps(pattern)
    text_index = pattern_index = 0
    step = 0

    while text_index < len(text):
        step += 1
        i_before, j_before = text_index, pattern_index
        text_char, pattern_char = text[text_index], pattern[pattern_index]

        if text_char == pattern_char:
            text_index += 1
            pattern_index += 1
            action = "match"
            if pattern_index == len(pattern):
                action = f"MATCH at {text_index - pattern_index}"
                pattern_index = lps[pattern_index - 1]
        elif pattern_index > 0:
            action = f"mismatch: j = lps[{j_before - 1}] = {lps[j_before - 1]}"
            pattern_index = lps[j_before - 1]
        else:
            action = "mismatch: j is 0, so i moves"
            text_index += 1

        print(f"{step:>2}  i={i_before:>2} j={j_before}  {text_char} vs {pattern_char}  {action}")


kmp_trace("ABABDABABABCABAB", "ABABCABAB")

 1  i= 0 j=0  A vs A  match
 2  i= 1 j=1  B vs B  match
 3  i= 2 j=2  A vs A  match
 4  i= 3 j=3  B vs B  match
 5  i= 4 j=4  D vs C  mismatch: j = lps[3] = 2
 6  i= 4 j=2  D vs A  mismatch: j = lps[1] = 0
 7  i= 4 j=0  D vs A  mismatch: j is 0, so i moves
 8  i= 5 j=0  A vs A  match
 9  i= 6 j=1  B vs B  match
10  i= 7 j=2  A vs A  match
11  i= 8 j=3  B vs B  match
12  i= 9 j=4  A vs C  mismatch: j = lps[3] = 2
13  i= 9 j=2  A vs A  match
14  i=10 j=3  B vs B  match
15  i=11 j=4  C vs C  match
16  i=12 j=5  A vs A  match
17  i=13 j=6  B vs B  match
18  i=14 j=7  A vs A  match
19  i=15 j=8  B vs B  MATCH at 7

Read the i column downwards: 0, 1, 2, 3, 4, 4, 4, 5, 6, … It repeats, but it never goes down.

Complexity

A single loop iteration is not bounded — one mismatch can trigger a whole chain of fallbacks. The bound comes from an amortised argument: those fallbacks have to be paid for by earlier work.

Building the table costs at most 2m comparisons. Watch length. It increases by 1 at most once per value of index, so across the whole build it increases at most m - 1 times. Every iteration of the inner while loop strictly decreases it, and it never drops below 0. A quantity that goes up at most m - 1 times in total can come down at most m - 1 times in total. So the inner loop runs at most m - 1 times across the entire build, not per index, and the total comparison count stays under 2m. Preprocessing is O(m).

Searching costs at most 2n comparisons. The same argument one level up. pattern_index increases only in the branch where the characters agree, and that branch also increases text_index, so it goes up at most n times in total — and therefore comes down at most n times in total. Every loop iteration performs exactly one comparison and then either increases text_index or decreases pattern_index, so there are at most n + n iterations. The scan is O(n).

Together: O(n + m) for the best case, the average case and the worst case alike. There is no adversarial input. Measure it on the pathological shape from the opening, scaled up.

def kmp_search_counted(text: str, pattern: str) -> tuple[list[int], int, int]:
    """KMP that reports (matches, comparisons building lps, comparisons scanning)."""
    lps = [0] * len(pattern)
    length = 0
    build_comparisons = 0

    for index in range(1, len(pattern)):
        while length > 0 and pattern[index] != pattern[length]:
            build_comparisons += 1
            length = lps[length - 1]
        build_comparisons += 1
        if pattern[index] == pattern[length]:
            length += 1
        lps[index] = length

    matches: list[int] = []
    scan_comparisons = 0
    text_index = pattern_index = 0

    while text_index < len(text):
        scan_comparisons += 1
        if text[text_index] == pattern[pattern_index]:
            text_index += 1
            pattern_index += 1
            if pattern_index == len(pattern):
                matches.append(text_index - pattern_index)
                pattern_index = lps[pattern_index - 1]
        elif pattern_index > 0:
            pattern_index = lps[pattern_index - 1]
        else:
            text_index += 1

    return matches, build_comparisons, scan_comparisons


haystack = "a" * 5000 + "b"
needle = "a" * 100 + "b"
naive_matches, naive_count = naive_search_counted(haystack, needle)
kmp_matches, build_count, scan_count = kmp_search_counted(haystack, needle)

print(f"text {len(haystack)} chars, pattern {len(needle)} chars")
print(f"naive: {naive_count} comparisons, matches {naive_matches}")
print(f"kmp:   {build_count + scan_count} comparisons "
      f"({build_count} to build the table, {scan_count} to scan), matches {kmp_matches}")
print(f"kmp bound 2n + 2m = {2 * len(haystack) + 2 * len(needle)}")

text 5001 chars, pattern 101 chars
naive: 495001 comparisons, matches [4900]
kmp:   10100 comparisons (199 to build the table, 9901 to scan), matches [4900]
kmp bound 2n + 2m = 10204

Forty-nine times fewer comparisons, and the measured 10,100 sits just under the predicted ceiling of 10,204. Naive search pays the full 101 comparisons at every one of its 4,901 alignments; KMP pays just under two per text character, however repetitive the input is. Double the text and the pattern together and naive search quadruples, to 1,970,001; KMP only doubles, to 20,200.

Space is O(m): one integer per pattern character, plus the two indices. The text is read, never copied — no slicing, no buffer. The one thing that does grow with the text is the returned list of match positions, and that is output rather than working memory: yield each hit instead of collecting it and the scan itself holds O(m) state no matter how long the text is.

Correctness deserves checking rather than trusting, so here is an exhaustive comparison against brute force and against Python's own str.find, over every text of up to 8 characters and every pattern of up to 3 characters from a two-letter alphabet. Small alphabets maximise partial matches, which is precisely where wrong border logic shows up.

from itertools import product

pairs = 0
for text_length in range(9):
    for text_letters in product("ab", repeat=text_length):
        sample = "".join(text_letters)
        for pattern_length in range(1, 4):
            for pattern_letters in product("ab", repeat=pattern_length):
                needle_text = "".join(pattern_letters)
                expected = naive_search(sample, needle_text)
                assert kmp_search(sample, needle_text) == expected
                first = expected[0] if expected else -1
                assert sample.find(needle_text) == first
                pairs += 1

print(f"{pairs} text/pattern pairs agree with brute force and with str.find")

7154 text/pattern pairs agree with brute force and with str.find

When to use it, and when not to

In application Python, reach for the built-in. pattern in text, text.find(pattern), text.count(pattern) and re.finditer are implemented in C and will beat a hand-written KMP loop by roughly two orders of magnitude, because the constant factor of a Python-level loop dwarfs any saving in comparisons. CPython already protects you from the quadratic blow-up too: its string search is a Boyer-Moore-Horspool style skip loop with a Bloom filter over the pattern's characters, and since Python 3.10 it switches to the Crochemore-Perrin "two-way" algorithm for long needles, which carries the same O(n + m) worst case guarantee. Not KMP, same protection.

KMP is not the fastest algorithm in practice either. It reads every character of the text at least once. Boyer-Moore-Horspool does not: it compares from the right of the pattern, and when the mismatching text character occurs nowhere in the pattern it jumps forward by the pattern's full length. On English prose with a 20 character pattern it touches roughly one character in ten. KMP cannot beat linear, so on ordinary text it loses. What KMP sells is the guarantee — Horspool's worst case is still O(n × m).

Write KMP when one of these is true:

You cannot rewind the input. A socket, a pipe, or a multi-gigabyte file you refuse to buffer. Because the text pointer never goes backwards, KMP consumes each character exactly once on arrival, with O(m) state and no window at all.
The input is adversarial or genuinely repetitive. Attacker-supplied patterns, DNA, run-length-ish binary formats.
You need the failure function for something else. It generalises directly to Aho-Corasick for many patterns at once, and m - lps[m - 1] gives the smallest period of a string for free.
Your language's built-in search is naive. Java's String.indexOf is a straightforward O(n × m) scan.

The streaming case is where KMP is genuinely the right tool rather than a teaching exercise:

from typing import Optional


class StreamMatcher:
    """Report matches while consuming a stream one character at a time.

    Only the pattern, its LPS table and two integers are kept, so a 4 GB log
    is scanned with state proportional to the pattern and never to the stream.
    """

    def __init__(self, pattern: str) -> None:
        if not pattern:
            raise ValueError("pattern must not be empty")
        self.pattern = pattern
        self.lps = build_lps(pattern)
        self.pattern_index = 0
        self.consumed = 0

    def feed(self, character: str) -> Optional[int]:
        """Return the start index of a match ending at this character, else None."""
        while self.pattern_index > 0 and character != self.pattern[self.pattern_index]:
            self.pattern_index = self.lps[self.pattern_index - 1]

        if character == self.pattern[self.pattern_index]:
            self.pattern_index += 1

        self.consumed += 1

        if self.pattern_index == len(self.pattern):
            self.pattern_index = self.lps[self.pattern_index - 1]
            return self.consumed - len(self.pattern)
        return None


matcher = StreamMatcher("ana")
hits = [start for start in (matcher.feed(char) for char in "bananana") if start is not None]
print(hits)

[1, 3, 5]

Three overlapping matches in bananana, found while holding nothing but the pattern, its three-entry table and two integers.

Do not use KMP for approximate matching. The border logic assumes exact character equality; allow typos or wildcards and the argument collapses, so use edit distance. For many patterns at once, build an Aho-Corasick automaton rather than running KMP once per pattern. If you search the same text repeatedly, do not scan at all — index it once with a suffix array or an FM-index.

One more honest number: on the 16 character example above, KMP costs 9 comparisons to build the table plus 19 to scan, against naive search's 26. Preprocessing is not free, and on short texts it does not pay for itself. KMP wins when n is large relative to m, or when one table is reused across many texts.

Where it shows up in the real world

Intrusion detection. The failure function generalises from one pattern to a trie of thousands; the result is the Aho-Corasick automaton, and it is a multi-pattern matching engine in both the Snort and Suricata intrusion detection systems — Suricata will pick Intel's Hyperscan instead when its build includes it. They test each packet against tens of thousands of attack signatures in a single pass over the payload, which is only possible because the automaton never re-reads a byte.

Small-alphabet sequence scanning. DNA has four letters and is full of tandem repeats, so partial matches are constant and long — the input shape that pushes naive search towards its quadratic case, and the one KMP is immune to. It is a reasonable tool for scanning reads for a fixed motif such as a restriction site or a sequencing adapter. Note the limit: the large read aligners such as BWA and Bowtie do not use KMP, because they answer millions of queries against one fixed reference genome, and it pays to build a Burrows-Wheeler/FM index of that genome once instead of streaming it per query.

Period detection. For a pattern of length m, the value m - lps[m - 1] is its shortest period, and the string is that period repeated exactly when the period divides m. abcabcabc has a final LPS entry of 6, and 9 − 6 = 3 divides 9, so it is abc three times.

Mostly, though, it is a building block. Aho-Corasick, the Z-algorithm and the general theory of string periodicity all descend from the failure function. Learning KMP is how you get access to them.

Common mistakes

Using lps[j] instead of lps[j - 1]. The table is indexed by last-character position, so the entry for a matched prefix of length j lives at j - 1. Getting it wrong shifts by the wrong amount, and the symptom depends on the pattern: searching ABABCABAB silently loses its match, searching aba reports a phantom hit at index 1 of abba, and searching aa never terminates.

Advancing the text pointer on a mismatch when j is above 0. The mismatched text character has not been consumed yet — the whole point is to re-test it against the shorter alignment. Advancing past it skips matches starting inside the region you already read.

Resetting j to 0 after a full match. It looks harmless and breaks overlapping searches: aa in aaaa returns two matches instead of three. Use lps[m - 1].

Letting the whole prefix count as its own border. If lps[index] could equal index + 1, the fallback length = lps[length - 1] would never shrink and the build loop would spin forever. "Proper" is doing real work in the definition.

Rebuilding the table inside the loop. Scanning 10,000 files for one pattern should call build_lps once, outside the file loop. Calling kmp_search per file is still correct and still linear, but it throws away the thing preprocessing buys you.

Forgetting the empty pattern. pattern[0] raises IndexError on an empty string. Decide what an empty pattern means, then handle it explicitly.

Practice

Add a first_only: bool = False parameter to kmp_search that returns as soon as it finds one match, and confirm the comparison count drops on a text with many occurrences.
Count the overlapping occurrences of aa in a string of 1,000 a characters, and explain from the LPS table why the answer is 999.
Write smallest_period(text) using len(text) - lps[-1], and use it to decide whether a string is some shorter string repeated a whole number of times.
Decide whether one string is a rotation of another with a single KMP search: check the lengths agree, then search the first inside the second concatenated with itself.
Extend the matcher so that ? in the pattern matches any character, then explain why the LPS table can no longer be built the same way.

Summary

KMP replaces "slide by one and start over" with "slide by exactly as much as the pattern's own structure permits". The table that says how much is the longest border of every prefix, built in O(m) by the same fallback move the search uses. The result is a search whose text pointer only ever moves right, and a hard O(n + m) bound with no worst case to fear — but also an algorithm that reads every character, which is why skip-based searches beat it on ordinary prose and why Python's own str.find is built from different parts.


Difficulty	Hard
Best case	O(n + m) — the bound is the same in every case
Average case	O(n + m)
Worst case	O(n + m) — at most 2m comparisons building, 2n scanning
Preprocessing	O(m) time, on the pattern only
Space	O(m) — one integer per pattern character
Text pointer	Never decreases, so it runs on a stream
Amortised	Yes — fallbacks are paid for by earlier matches
Data structure	String plus the LPS array
Use it when	Input is repetitive or adversarial, cannot be rewound, or feeds Aho-Corasick
Avoid it when	A C-level `str.find` will do, or the text is ordinary prose
Real-world use	Aho-Corasick in Snort and Suricata; motif scanning; period detection
Python equivalent	`pattern in text`, `text.find(pattern)`, `re.finditer`

Keep reading

Rabin-Karp — the other linear-average substring search, built from rolling hashes instead of borders.
Tries — the structure Aho-Corasick hangs its failure links on, and the natural next step from here.
Big O Notation — the full treatment of the amortised argument used to prove the 2n bound.
The Two Pointers Technique — the same "each pointer only moves one way" discipline, applied to arrays.
Edit Distance — what to reach for when the match does not have to be exact.

The N-Queens Problem in Python: Backtracking at Its Clearest

Bimal Kshetri — Sun, 16 Aug 2026 16:36:56 +0000

Put eight queens on a chessboard so that none of them can capture another. That is the entire problem, posed by the chess composer Max Bezzel in 1848, and it has outlived almost every other puzzle of its era for one reason: it is the cleanest demonstration of backtracking anyone has found. The general version swaps eight for n — place n queens on an n by n board with no two sharing a row, a column or a diagonal.

The reason it is worth your time is the gap between the obvious approach and the good one. Choosing 8 squares out of 64 gives 4,426,165,368 boards to test. The good version examines 2,057 partial boards and finds all 92 solutions. Nothing clever was added to get there — the entire improvement comes from two things: describing the board in a way that makes illegal positions unrepresentable, and abandoning a branch the instant it becomes hopeless.

Be clear about what this is, though. N-queens is a benchmark and a teaching problem, not a production algorithm. If you only need one solution for large n there are closed-form constructions that place the queens directly in O(n) time with no search at all. Backtracking earns its place when you need every solution, an exact count, or the freedom to bolt extra constraints onto the search.

The idea

A queen attacks along four lines at once: her whole row, her whole column, and both diagonals running through her square. Every square on any of those lines is off limits to every other queen.

Look at how much one queen destroys. On a 4 by 4 board a queen placed near the middle attacks 11 of the other 15 squares, leaving 4. That destructiveness is why the naive search space is such a lie: almost none of it is reachable.

Start counting anyway, because the numbers are what justify everything that follows.

from math import comb, factorial

for n in (4, 6, 8):
    squares = n * n
    print(f"n = {n}:  {comb(squares, n):>13,}  boards with n queens dropped anywhere")
    print(f"        {n ** n:>13,}  boards with exactly one queen per row")
    print(f"        {factorial(n):>13,}  boards with one per row and one per column")

n = 4:          1,820  boards with n queens dropped anywhere
                  256  boards with exactly one queen per row
                   24  boards with one per row and one per column
n = 6:      1,947,792  boards with n queens dropped anywhere
               46,656  boards with exactly one queen per row
                  720  boards with one per row and one per column
n = 8:  4,426,165,368  boards with n queens dropped anywhere
           16,777,216  boards with exactly one queen per row
               40,320  boards with one per row and one per column

Two observations collapse that first column into the third, and both are worth stating precisely because they do all the work.

Exactly one queen goes in every row. There are n queens and n rows, and two queens in the same row would attack each other, so no row can hold two. n queens spread over n rows with at most one each means every row holds exactly one. So a candidate board is fully described by n numbers: for each row, which column its queen sits in. That single change takes 4.4 billion boards down to 16.7 million for n = 8.

No column is used twice either, by exactly the same argument. So those n numbers are a permutation of 0 through n − 1. That takes 16.7 million down to 40,320.

Everything left is diagonals. And rather than generate all 40,320 permutations and filter them, build the permutation one row at a time and reject as early as possible: place a queen in row 0, then row 1, then row 2, and the moment a row has no legal square, throw away the last queen and try the next column for it instead.

That removal step is what makes it backtracking rather than plain recursion. The board must be restored to exactly its previous state before the next column is tried, or later branches inherit queens that are no longer there.

Naming the diagonals

The whole algorithm now rests on one question, asked once per candidate square: is this square attacked by any queen already placed? If answering it costs a scan over the placed queens, that is O(n) per square and O(n²) per row. There is a much better way.

Give every diagonal a name. Walk one step down and to the right, and both the row and the column go up by one, so row - col does not change. That means every square on a \ diagonal shares one value of row - col, and squares on different \ diagonals never do.

The other direction works the same way with the other sign. Walk one step down and to the left and the row goes up by one while the column goes down by one, so row + col does not change. Every / diagonal has its own constant value of row + col.

On an n by n board row - col ranges from −(n − 1) to n − 1 and row + col ranges from 0 to 2n − 2 — 2n − 1 diagonals in each direction, 30 in total for a chessboard. So keep three sets:

used_columns — the columns already holding a queen.
used_down_diagonals — the row - col values already taken.
used_up_diagonals — the row + col values already taken.

A square at (row, col) is safe exactly when col, row - col and row + col are all absent from their respective sets. Three set lookups, each O(1) on average because Python hashes small integers to themselves. No scan, no matter how many queens are already down. That is the trick the whole implementation is built on.

Watching it work

Take n = 4 and follow the search exactly as the code will run it, always trying columns left to right.

Row 0, column 0. The sets become columns {0}, down {0}, up {0}.

Row 1. Column 0 is a used column. Column 1 has row - col = 0, already taken by the queen at (0, 0) — they are on the same \ diagonal. Column 2 is clear on all three counts, so the queen goes there. The sets are now columns {0, 2}, down {0, -1}, up {0, 3}.

Row 2. Every square fails:

Column 0 — the column is used.
Column 1 — row + col = 3, which the queen at (1, 2) already owns.
Column 2 — the column is used.
Column 3 — row - col = -1, which the queen at (1, 2) already owns.

That is the first backtrack. The search returns to row 1, removes the queen from column 2 — restoring the sets to columns {0}, down {0}, up {0} — and tries column 3 instead.

From there it gets one row further. Row 2 accepts column 1, then row 3 has nothing left: columns 0, 3 and 1 are used, and column 2 sits on row - col = 1, which the queen at (2, 1) owns. Back up to row 2 — column 2 lies on the same \ diagonal as the queen at (0, 0) and column 3 is taken, so row 2 is finished too. Row 1 then runs out of columns. The entire subtree under "row 0, column 0" contains no solution, and the search finally moves the first queen to column 1.

That branch works out immediately: row 1 takes column 3, row 2 takes column 0, row 3 takes column 2. Four queens, no conflicts, first solution found.

The remaining two branches are mirror images of the first two: column 2 mirrors column 1 and yields the second solution, column 3 mirrors column 0 and yields nothing. Seventeen partial boards examined in total, against 24 permutations and 1,820 raw placements.

The code

Start with the version that follows straight from the counting above: generate every permutation of the columns, then check the diagonals. It is short, obviously correct, and a useful reference to test the fast version against.

from itertools import permutations


def is_safe_arrangement(columns: tuple[int, ...]) -> bool:
    """True when no two queens in a one-per-row, one-per-column board share a diagonal."""
    for upper in range(len(columns)):
        for lower in range(upper + 1, len(columns)):
            # Two queens share a diagonal when the row gap equals the column gap.
            if lower - upper == abs(columns[lower] - columns[upper]):
                return False
    return True


def count_by_permutation(n: int) -> int:
    """Count solutions by testing every column permutation. Correct, and wasteful."""
    return sum(1 for columns in permutations(range(n)) if is_safe_arrangement(columns))


for n in range(4, 9):
    print(f"n = {n}: {count_by_permutation(n):>3} solutions from {factorial(n):>6,} permutations")

n = 4:   2 solutions from     24 permutations
n = 5:  10 solutions from    120 permutations
n = 6:   4 solutions from    720 permutations
n = 7:  40 solutions from  5,040 permutations
n = 8:  92 solutions from 40,320 permutations

Those are the right answers — 92 for the classic eight-queens board — and note that the counts are not monotonic: six queens have only 4 solutions where five queens have 10. The problem is spikier than it looks.

This version dies quickly, though. permutations is a generator, so it hands the boards over one at a time rather than materialising all of them — but it still hands over every one of the n!, and nothing is judged until all n queens are already placed. n = 12 means 479 million full boards, and n = 15 means over a trillion. Worse, it learns nothing from failure: for n = 8 a permutation starting 0, 1 already has two queens on the same diagonal, but the loop still produces all 720 completions of that prefix and tests every one of them from scratch.

Backtracking fixes precisely that. Build the permutation left to right and check each queen against the ones already placed, so a doomed prefix is abandoned once instead of re-derived thousands of times.

def solve_n_queens(n: int) -> list[list[int]]:
    """Every placement of n queens on an n by n board with no two attacking.

    A solution is a list of n column indices: entry r holds the column of the
    queen in row r. One queen per row is baked into the shape of the answer,
    so the search never considers anything else.
    """
    solutions: list[list[int]] = []
    used_columns: set[int] = set()
    used_down_diagonals: set[int] = set()   # row - col is constant along a "\" line
    used_up_diagonals: set[int] = set()     # row + col is constant along a "/" line
    placement: list[int] = []

    def place_queen_in(row: int) -> None:
        if row == n:
            # Every row is filled, so this partial board is a complete solution.
            solutions.append(placement.copy())
            return

        for col in range(n):
            if (col in used_columns
                    or row - col in used_down_diagonals
                    or row + col in used_up_diagonals):
                continue

            used_columns.add(col)
            used_down_diagonals.add(row - col)
            used_up_diagonals.add(row + col)
            placement.append(col)

            place_queen_in(row + 1)

            # Undo before trying the next column, so the three sets always
            # describe exactly the queens still standing on the board.
            placement.pop()
            used_up_diagonals.remove(row + col)
            used_down_diagonals.remove(row - col)
            used_columns.remove(col)

    place_queen_in(0)
    return solutions


print(solve_n_queens(4))
print(f"n = 1: {solve_n_queens(1)}   n = 2: {solve_n_queens(2)}   n = 3: {solve_n_queens(3)}")

[[1, 3, 0, 2], [2, 0, 3, 1]]
n = 1: [[0]]   n = 2: []   n = 3: []

The two solutions for n = 4 are the ones traced by hand above, in the order the search finds them. The small boards behave correctly with no special cases: one queen on a 1 by 1 board is trivially fine, and 2 and 3 genuinely have no solutions, which the code discovers by exhausting the search rather than by being told.

A list of column indices is a compact answer but a poor picture, so draw it.

def render(solution: list[int]) -> str:
    """Draw one solution: Q for a queen, . for an empty square."""
    return "\n".join(
        " ".join("Q" if col == queen_col else "." for col in range(len(solution)))
        for queen_col in solution
    )


for number, solution in enumerate(solve_n_queens(4), start=1):
    print(f"n = 4, solution {number}: columns {solution}")
    print(render(solution))

first_eight = solve_n_queens(8)[0]
print(f"n = 8, first solution: columns {first_eight}")
print(render(first_eight))

n = 4, solution 1: columns [1, 3, 0, 2]
. Q . .
. . . Q
Q . . .
. . Q .
n = 4, solution 2: columns [2, 0, 3, 1]
. . Q .
Q . . .
. . . Q
. Q . .
n = 8, first solution: columns [0, 4, 7, 5, 2, 6, 1, 3]
Q . . . . . . .
. . . . Q . . .
. . . . . . . Q
. . . . . Q . .
. . Q . . . . .
. . . . . . Q .
. Q . . . . . .
. . . Q . . . .

How the code maps to the idea

The row parameter is the recursion depth. There is no loop over rows anywhere, because "one queen per row" is expressed by the structure of the recursion: each call handles exactly one row, and place_queen_in(row + 1) moves to the next. Nothing in the code can produce a board with two queens in a row, so no check for it exists.

The base case is row == n. Reaching it means all n rows were filled without a conflict, so the current placement is a solution. There is no validity test at the bottom — every queen was checked before it was placed, so a complete board is correct by construction.

placement.copy() is load-bearing. placement is one list that the search mutates for the whole run. Appending it directly would store one more reference to that same list every time a solution is found — 92 of them for n = 8 — and the list is emptied again on the way out, so every one of those references ends up pointing at []. Copy it, or the function returns a pile of aliases to nothing.

The continue is the pruning. When a square fails the three-set test, no queen is placed and no recursion happens, so every board that would have grown out of that square is skipped. That is the entire performance story — the whole subtree disappears from one if.

The three remove calls are the backtrack. They mirror the three add calls, undoing them in reverse order. The order does not matter functionally — the three sets are independent — but mirroring makes it obvious that nothing was missed. Forget one and the sets slowly fill with ghosts, and the search silently reports too few solutions.

The nested function closes over the sets, so every recursive call shares one copy of each. Passing copies down would also be correct, and would cost O(n) per call for nothing. Depth is only n + 1 frames, so Python's 1,000-frame recursion limit never binds.

Complexity

The honest headline: the upper bound is loose, the real cost is much smaller, and nobody can state the real cost exactly.

The upper bound. Because the search fixes one queen per row and never repeats a column, the nodes at depth k are at most n × (n − 1) × … × (n − k + 1), which is n! / (n − k)!. Summing over all depths gives a total node count of at most

n! × (1/0! + 1/1! + … + 1/n!), which is less than e × n! ≈ 2.72 n!.

Each node loops over n columns and does three O(1) set lookups per column, so the work per node is O(n). Multiply: O(n × n!) in the worst case, and that is genuinely an upper bound rather than a description of what happens.

What actually happens. The pruning removes most of that tree, and there is no known closed form for how much. So measure it.

def count_n_queens(n: int) -> tuple[int, int]:
    """Return (number of solutions, number of partial boards the search examined)."""
    used_columns: set[int] = set()
    used_down_diagonals: set[int] = set()
    used_up_diagonals: set[int] = set()
    boards_examined = 0

    def place_queen_in(row: int) -> int:
        nonlocal boards_examined
        boards_examined += 1
        if row == n:
            return 1

        found = 0
        for col in range(n):
            if (col in used_columns
                    or row - col in used_down_diagonals
                    or row + col in used_up_diagonals):
                continue

            used_columns.add(col)
            used_down_diagonals.add(row - col)
            used_up_diagonals.add(row + col)
            found += place_queen_in(row + 1)
            used_up_diagonals.remove(row + col)
            used_down_diagonals.remove(row - col)
            used_columns.remove(col)

        return found

    return place_queen_in(0), boards_examined


print(f"{'n':>3} {'solutions':>10} {'boards examined':>16} {'n!':>14}")
for n in range(4, 13):
    solutions, boards = count_n_queens(n)
    print(f"{n:>3} {solutions:>10,} {boards:>16,} {factorial(n):>14,}")

  n  solutions  boards examined             n!
  4          2               17             24
  5         10               54            120
  6          4              153            720
  7         40              552          5,040
  8         92            2,057         40,320
  9        352            8,394        362,880
 10        724           35,539      3,628,800
 11      2,680          166,926     39,916,800
 12     14,200          856,189    479,001,600

At n = 12 the bound allows about 1.3 billion nodes and the search examines 856,189 — roughly 1,500 times fewer. But watch the ratios down that middle column: 2,057, then 8,394, then 35,539, then 166,926, then 856,189. Each step multiplies the work by more than the step before it: 4.1×, then 4.2×, then 4.7×, then 5.1×. Pruning changed the constant and the base, not the shape: the growth is still worse than exponential. On this code n = 14 already takes minutes and n = 15 tens of minutes, and n = 20 is out of reach no matter how long you wait.

Space: O(n). The recursion is n + 1 frames deep and the three sets hold at most n entries each. count_n_queens therefore runs in linear space no matter how astronomically many solutions it counts. solve_n_queens is different — it stores every solution, so its memory is O(n × number of solutions), which for n = 12 is 14,200 lists of 12 integers. If you only want the count, count.

One more honest note: counting solutions is far harder than finding one. The exact totals are known only up to n = 27, and that value was produced in 2016 by a distributed FPGA project that ran for about a year. No formula for the sequence is known.

When to use it, and when not to

Use backtracking for n-queens when you want all solutions, an exact count, or a solution satisfying extra constraints you invented ("no queen on the main diagonal", "these three squares must be occupied"). It is also the right shape for the whole family of constraint puzzles — Sudoku, graph colouring, crosswords, the knight's tour — where the three-set pattern generalises directly.

Do not use it to find a single solution for large n. Explicit constructions exist that place n non-attacking queens for every n ≥ 4 in O(n) time by writing the column indices from a formula, with no search whatsoever. If all you need is one valid board, searching for it is the wrong tool by an enormous margin.

Do not use it for very large n even with the constraints relaxed. Local search beats it badly here: the min-conflicts heuristic starts from a random full board, picks a queen that is currently under attack, and moves it within its own column to the square that conflicts with the fewest others. It solves the million-queens problem in around fifty moves on average. Backtracking cannot get past a few dozen.

Do not hand-roll it when the constraints get complicated. Once the rules stop being three simple sets, a real constraint solver — Google's OR-Tools CP-SAT, or a SAT solver — will out-prune anything you write by hand, because it learns from conflicts instead of just backing up one level.

And if the search tree has overlapping subproblems rather than merely infeasible ones, backtracking is the wrong family entirely; that is what dynamic programming is for. N-queens has a trace of overlap: on an 8 by 8 board the prefixes 1, 3, 0, 2 and 2, 0, 3, 1 block exactly the same columns and the same diagonals, so rows 4 through 7 face an identical subproblem. There is far too little of it to pay for a memo table, though. The key would have to be a set of columns plus two sets of diagonals, and that key repeats so rarely that the cache needs an entry for very nearly every node — for which it cuts n = 8 from 2,057 nodes to 1,999, and n = 12 from 856,189 to 821,299. Under five per cent of the work saved, for a table almost as large as the search itself. Plain backtracking stays the right family here.

Where it shows up in the real world

Not as itself. Nobody ships software that places queens on chessboards.

What it has instead is a real career as a proving ground. Niklaus Wirth built his 1971 paper Program Development by Stepwise Refinement — one of the founding documents of structured programming — around the eight queens problem, and the modern shape of this algorithm is essentially his. Constraint toolkits including Google's OR-Tools and MiniZinc ship it as a standard example, because everyone already knows the answers, so a new solver can be checked against them. Russell and Norvig's Artificial Intelligence: A Modern Approach uses it as the running example for local search, which is where the million-queens figure comes from.

The transferable part is the pattern, and that does ship. A Sudoku solver is this algorithm with different bookkeeping: three collections of used digits, one per row, one per column and one per 3 by 3 box, with the same place-recurse-undo cycle around them. Once you can write n-queens without thinking, you can write that in twenty minutes.

Common mistakes

Forgetting to undo. Leave out one of the three remove calls and the sets fill with queens that are no longer on the board. Nothing crashes; the search just reports too few solutions. Drop the down-diagonal remove and n = 8 comes back as 0 rather than 92, because the ghosts block every branch the search tries after its first retreat.

Storing the live list. solutions.append(placement) instead of solutions.append(placement.copy()) stores the same object every time. You get one reference per solution — 92 of them at n = 8 — all pointing at the same list, and that list is empty by the time the search returns, so the answer prints as ninety-two copies of [].

Using the same expression for both diagonals. If used_down is also fed row + col, the code compiles, runs, and reports 2,113 solutions for n = 8 — it is now only checking columns and one diagonal direction. Always check n = 8 against 92; it is the cheapest regression test in the world.

Sizing the diagonal arrays by n. Swapping the sets for boolean lists is a reasonable optimisation, but there are 2n − 1 diagonals in each direction, not n. Allocating [False] * n and indexing it by row - col raises nothing, because Python reads negative indices from the end — diagonal v and diagonal v − n end up sharing a slot, so the search blocks squares it should not and under-counts. That does not fail loudly and it does not fail uniformly: n = 8 comes back as 0 instead of 92 and n = 9 as 116 instead of 352, but n = 5 comes back as 10, which is the correct answer. Test the bug on a small odd board and it passes. A list of length exactly 2 * n - 1 does work, since the negative wrap covers exactly one period, but used_down[row - col + n - 1] on the same list is the version you can reread in six months.

Scanning the placed queens for every candidate square. Looping over placement and testing abs(row - other_row) == abs(col - other_col) is correct and is what the permutation version does. It also turns an O(1) check into an O(n) one, which multiplies the whole runtime by n for no benefit.

Checking rows. Some implementations carry a used_rows set. It can never fire — the recursion visits each row exactly once — so it is pure overhead, and its presence usually signals a misunderstanding of why the one-queen-per-row rule was free.

Practice

Add an early return so the search stops at the first solution instead of enumerating all of them, and confirm it returns [1, 3, 0, 2] for n = 4.
Replace the three sets with three boolean lists, remembering the + n - 1 shift for the down diagonals, and check that every count from n = 4 to n = 10 is unchanged.
Print how many of the 92 eight-queens solutions have a queen in the corner square, using the rendering function to spot-check a few by eye.
For even n, restrict the queen in row 0 to the left half of the board and double the resulting count; verify it matches the full search for n = 8, 10 and 12, and work out why the trick needs care for odd n.
Reuse the place-recurse-undo skeleton to solve a 9 by 9 Sudoku, tracking used digits per row, per column and per 3 by 3 box.

Summary

N-queens is the problem to reach for when you want to understand backtracking, because every part of it is visible. The search space shrinks from 4.4 billion to 40,320 through two sentences of reasoning, the conflict test collapses to three O(1) lookups because row - col and row + col name the diagonals, and the pruning turns an intractable bound into 2,057 examined boards. Learn the three-set trick properly; you will use it again the next time a puzzle has "no two of these may share a line" in its rules.


Difficulty	Hard
Time, upper bound	O(n × n!) — fewer than e × n! partial boards, O(n) work at each
Time, measured	2,057 boards at n = 8; 856,189 at n = 12, about 1,500× below the bound
Space	O(n) — n + 1 recursion frames plus three sets of at most n entries
Space to list all solutions	O(n × number of solutions) — 14,200 lists of 12 ints at n = 12
Conflict check	O(1) — three set lookups, never a scan of placed queens
Key insight	One queen per row; `row - col` and `row + col` identify the two diagonals
Data structure	Three sets plus one list of column indices
Use it when	You need all solutions, an exact count, or custom extra constraints
Avoid it when	You need one solution for large n — use the O(n) construction or min-conflicts
Real-world use	A benchmark for constraint and SAT solvers; the pattern behind Sudoku solvers
Python equivalent	None in the standard library; `itertools.permutations` gives the brute force

Keep reading

Recursion and Backtracking in Python — the mental model this post applies, built from scratch.
Depth-First Search — the same traversal, run over an explicit graph instead of an implicit search tree.
Hash Tables in Python — why those three set lookups really are O(1), and when they are not.
Dynamic Programming Explained — what to do instead when the branches of the search overlap.
Big O Notation — the counting argument behind the n! bound above.

Recursion and Backtracking in Python: Building the Mental Model

Bimal Kshetri — Sun, 16 Aug 2026 16:36:25 +0000

A recursive function is one that calls itself. That sentence is the whole definition, and it is also why recursion confuses people: it sounds circular, like a dictionary that defines "recursion" as "see recursion". It is not circular, because every call works on a smaller problem than the one that made it, and because there is always a smallest problem the function answers outright instead of asking again.

Once that clicks, a large part of this series stops being a list of tricks. Merge sort, quick sort, tree traversal, depth-first search, dynamic programming, divide and conquer, and every puzzle solver ever written are the same shape: solve the trivial case directly, break everything else into smaller versions of itself, combine the answers.

Two halves follow. First plain recursion: the base case, the call stack drawn frame by frame, why a missing base case crashes instead of hanging, and where Python's limits sit. Then backtracking, which is recursion plus one move — undo the choice you just made and try the next one. That single addition is how you enumerate every subset, every permutation, every legal Sudoku grid.

The idea

A recursive function has exactly two parts, and if either one is missing the function is broken.

The base case. A version of the problem so small that you answer it directly, with no further calls. factorial(1) is 1. An empty list has length 0. A leaf node has no children to visit.
The recursive case. Everything else. You express the answer in terms of the same function applied to something strictly smaller, then do the small amount of work that turns that answer into yours.

Factorial is the standard first example because it is defined recursively in mathematics before anyone thinks about code. 4! means 4 × 3 × 2 × 1, which is the same as 4 × 3!. In general n! = n × (n - 1)!, and 1! = 1 stops it.

The part beginners find genuinely hard is not writing that down. It is believing that the machine can keep track of four half-finished multiplications at once. It can, because of the call stack.

Every time a Python function is called, the interpreter allocates a frame: a small block of memory holding that call's arguments, its local variables, and the place to resume once it returns. Frames pile up. The frame at the top is the one currently running; everything below it is paused, waiting for the value it asked for.

That picture is the entire mystery. factorial(4) cannot finish its multiplication until it knows factorial(3), so it sits there holding the number 4 while a new frame runs. Four frames exist at once, each with its own private n. Nothing is shared, nothing is overwritten.

Two rules follow, and every recursion bug is a violation of one of them:

The base case must exist and must be reachable. Not just present in the source — actually reached by the arguments you pass.
Every recursive call must move strictly closer to it. factorial(n - 1) shrinks. factorial(n) does not. Neither does factorial(n - 1) when n is negative and the base test is n == 1.

Watching it work

Trace factorial(4) by hand. Going down, each call defers its multiplication:

factorial(4) = 4 * factorial(3)      frame 1 opens, holds 4
factorial(3) = 3 * factorial(2)      frame 2 opens, holds 3
factorial(2) = 2 * factorial(1)      frame 3 opens, holds 2
factorial(1) = 1                     frame 4 hits the base case

Nothing has been multiplied yet. Four frames are open and three of them are stuck on a half-written expression. Then the base case returns and the stack unwinds, in exactly the reverse order:

factorial(1) hands 1 back to factorial(2), which finally computes 2 × 1 = 2 and hands that to factorial(3), which computes 3 × 2 = 6, which lets factorial(4) compute 4 × 6 = 24. The work happens on the way up, not on the way down. That asymmetry is worth remembering — it is why a recursive function can build an answer in an order you never explicitly wrote.

Fibonacci changes the picture in one important way: the recursive case calls itself twice. fib(n) is fib(n - 1) + fib(n - 2), with fib(0) = 0 and fib(1) = 1. One call per level becomes a branching tree of calls.

Fifteen calls in total to compute fib(5). But look at the highlighted chain: at the instant the leftmost fib(1) is running, only five frames exist. The other ten calls have not started yet. That distinction matters:

Total calls determine the running time. For fib they explode.
Maximum depth determines the memory, because that is the tallest the stack ever gets. For fib it is only n.

Count them separately, because they diverge. fib(30) makes 2,692,537 calls in a stack only 30 frames tall and runs fine; factorial(2000) makes 2,000 calls in a stack 2,000 frames tall and crashes.

The code

Factorial first, exactly as described:

def factorial(n: int) -> int:
    """n! computed from its own definition: n times the factorial below it."""
    if n <= 1:  # base case: 0! and 1! are both 1, and neither needs recursion
        return 1
    return n * factorial(n - 1)  # recursive case: the same job, one size smaller


print(factorial(4))
print([factorial(k) for k in range(7)])

24
[1, 1, 2, 6, 24, 120, 720]

If you do not yet believe the stack diagram, make the function narrate itself. Indentation here is literally the stack depth:

def factorial_traced(n: int, depth: int = 0) -> int:
    """factorial again, announcing each frame as it opens and as it returns."""
    pad = "  " * depth
    print(f"{pad}call factorial({n})")
    if n <= 1:
        print(f"{pad}base case, return 1")
        return 1
    result = n * factorial_traced(n - 1, depth + 1)
    print(f"{pad}return {n} * factorial({n - 1}) = {result}")
    return result


factorial_traced(4)

call factorial(4)
  call factorial(3)
    call factorial(2)
      call factorial(1)
      base case, return 1
    return 2 * factorial(1) = 2
  return 3 * factorial(2) = 6
return 4 * factorial(3) = 24

Four calls down, four returns up, and every return line is a frame finishing the multiplication it started. Printing a traced version of a recursion you do not understand is the fastest debugging technique in this entire series.

Now Fibonacci, with a second function that counts the calls instead of guessing at them:

def fib(n: int) -> int:
    """The nth Fibonacci number, written straight from the definition."""
    if n < 2:  # fib(0) = 0 and fib(1) = 1 are given, not computed
        return n
    return fib(n - 1) + fib(n - 2)


def fib_call_count(n: int) -> int:
    """Total calls the function above makes while computing fib(n)."""
    if n < 2:
        return 1
    return 1 + fib_call_count(n - 1) + fib_call_count(n - 2)


for k in (5, 10, 20, 30):
    print(f"fib({k:>2}) = {fib(k):>6}   calls: {fib_call_count(k):>9}   deepest stack: {k}")

fib( 5) =      5   calls:        15   deepest stack: 5
fib(10) =     55   calls:       177   deepest stack: 10
fib(20) =   6765   calls:     21891   deepest stack: 20
fib(30) = 832040   calls:   2692537   deepest stack: 30

Adding ten to n multiplies the calls by about 123 while the stack grows by ten frames. Correct, and unusable past about n = 35. The fix is not to abandon recursion — it is to stop recomputing answers you already have, which is dynamic programming, and in Python it is one decorator: functools.lru_cache.

How the code maps to the idea

The base test comes first. Always. If the recursive call is written above the base test, the call happens before the test can stop it, and nothing stops it.

n - 1 is the contract. Each call receives an argument strictly closer to the base case, and the base case catches everything at or below the boundary. Writing if n == 1 instead of if n <= 1 looks equivalent and is not: factorial(0) then recurses to −1, −2, −3 and never matches.

The return value is the only channel. Each frame gets one value back from the frame above it and is responsible for turning that into its own answer. factorial multiplies. fib adds. Merge sort merges. That final combining step is where a divide-and-conquer algorithm actually does its work.

fib is exponential because the branches overlap, not because it recurses. fib(5) computes fib(3) twice, fib(2) three times, fib(1) five times. Recursion is not slow. Recomputation is slow.

Where recursion breaks

Delete the base case and something interesting happens: the program does not hang. It stops, quickly, with an exception.

import sys


def countdown(n: int) -> int:
    """No base case: every call makes another call, one frame deeper."""
    return countdown(n - 1)


try:
    countdown(5)
except RecursionError:
    print("countdown(5) stopped with RecursionError - it never hung")

print("frames allowed by default:", sys.getrecursionlimit())

countdown(5) stopped with RecursionError - it never hung
frames allowed by default: 1000

An infinite loop spins forever using no extra memory. An infinite recursion consumes a frame per call, so CPython counts frames and raises RecursionError at 1000 of them. That limit is a guard rail, not a law of nature — sys.setrecursionlimit(20000) raises it — but raise it only when the depth is genuinely bounded and merely larger than 1000. Frames are real memory, and on older CPython builds each Python call also consumed C stack, so a limit set high enough crashes the interpreter outright instead of raising a catchable exception.

The ceiling has one blunt consequence: a recursion whose depth grows with input size is a bug waiting for a big input. Walking a linked list of 5,000 nodes recursively passes your tests and dies in production.

Functional languages solve this with tail-call optimisation. A call is in tail position when it is the entire return expression — nothing is left to do after it comes back — so the compiler can reuse the current frame instead of stacking a new one. Scheme guarantees it. Python does not do it at all:

def sum_to(n: int, total: int = 0) -> int:
    """Tail recursive: the recursive call is the entire return expression."""
    if n == 0:
        return total
    return sum_to(n - 1, total + n)


print(sum_to(800))

try:
    print(sum_to(5000))
except RecursionError:
    print("sum_to(5000) overflows anyway - Python keeps every frame")


def sum_to_loop(n: int) -> int:
    """The same arithmetic as a loop: one frame, any n."""
    total = 0
    for value in range(1, n + 1):
        total += value
    return total


print(sum_to_loop(5000))
print(sum_to_loop(1_000_000))

320400
sum_to(5000) overflows anyway - Python keeps every frame
12502500
500000500000

sum_to is textbook tail recursion and it still blows the stack, while the loop handles a million in one frame. Guido van Rossum, Python's creator, wrote about this on his Neopythonic blog in 2009 and rejected the feature deliberately, not for lack of time. His main objection was debuggability: eliminating frames destroys the traceback, and a stack trace that has silently dropped its middle is a much worse tool than a slightly slower program. He also argued that Python is not a functional language, that iteration is the idiomatic way to loop, and that programmers should not have to reason about whether a call is in tail position to know whether their code will crash.

So in Python, the rule is blunt. Use recursion when the depth is naturally small — the height of a balanced tree, the number of digits in a number, log n levels of divide and conquer. Use a loop or an explicit stack when the depth scales with the data — which is exactly why the iterative version of depth-first search exists.

Backtracking: choose, explore, un-choose

Everything above computes a single value. Backtracking answers a different kind of question: enumerate every arrangement that satisfies some constraints, or find one that does.

The mechanism is three lines, repeated at every level:

Choose. Commit to one option at the current position.
Explore. Recurse to fill the next position, given that commitment.
Un-choose. When the recursion returns, take the commitment back, so the next option starts from a clean state.

The un-choose step is the whole trick, and it is what makes this backtracking rather than plain recursion. The recursion explores one branch to its very end, then rewinds to the last decision point and takes the other road — depth-first search over a tree of decisions that you never build in memory. The tree exists only as the sequence of calls.

Start with subsets. To list every subset of [1, 2, 3] you decide, for each element in turn, take it or leave it. Three elements, two choices each, 2³ = 8 subsets.

def subsets(items: list[int]) -> list[list[int]]:
    """Every subset of items, by taking or leaving each element in turn."""
    result: list[list[int]] = []
    path: list[int] = []

    def explore(index: int) -> None:
        if index == len(items):
            result.append(list(path))  # copy, because path keeps changing
            return
        path.append(items[index])  # choose items[index]
        explore(index + 1)         # explore every ending that includes it
        path.pop()                 # un-choose, putting path back as it was
        explore(index + 1)         # explore every ending that leaves it out

    explore(0)
    return result


print(subsets([1, 2, 3]))
print(len(subsets([1, 2, 3, 4, 5])))

[[1, 2, 3], [1, 2], [1, 3], [1], [2, 3], [2], [3], []]
32

Read the output against the diagram: the leaves come out left to right, take-branch first. index == len(items) is the base case — every element has been decided, so the path is a finished subset. path.append is the choice, path.pop is the undo, and the two explore calls are the two branches out of every node.

Permutations, and one list rewound

Permutations are the same pattern with a different constraint. Instead of two choices per position, every unused element is a candidate, and each element must be used exactly once.

from itertools import permutations as stdlib_permutations


def permutations(items: list[int]) -> list[list[int]]:
    """Every ordering of items, choosing one unused element per position."""
    result: list[list[int]] = []
    path: list[int] = []
    used = [False] * len(items)

    def explore() -> None:
        if len(path) == len(items):
            result.append(list(path))
            return
        for index, value in enumerate(items):
            if used[index]:  # each element fills exactly one position
                continue
            used[index] = True   # choose
            path.append(value)
            explore()            # explore
            path.pop()           # un-choose, both halves of the state
            used[index] = False

    explore()
    return result


print(permutations([1, 2, 3]))
print(len(permutations([1, 2, 3, 4, 5, 6])))
print(permutations("abc") == [list(p) for p in stdlib_permutations("abc")])

[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
720
True

That last True says the hand-written function produces the same orderings in the same order as itertools.permutations. In real code use the standard library version: it is written in C and yields results lazily instead of building a list of 720. Write your own when the constraints are custom, which is exactly when backtracking earns its place.

Notice that two things get undone in the loop: path.pop() and used[index] = False. Every piece of state you mutate on the way down has to be restored on the way up. Miss one and later branches inherit a corrupted world.

That diagram answers the question everyone asks eventually: when do you copy, and when do you mutate and undo?

Mutate and undo the working state. There is one path list for the entire search, and one used array. Appending is amortised O(1), popping is O(1), and no new list is allocated per node. This is why backtracking is fast despite visiting a huge number of nodes.
Copy when you record an answer. result.append(list(path)) makes a snapshot. Writing result.append(path) instead stores a reference to the one shared list, which keeps changing — so at the end every entry in result is the same object, and subsets([1, 2, 3]) returns eight empty lists. This is the single most common backtracking bug.
Copy when you pass state down only if you cannot undo it cleanly. Passing path + [value] into the recursive call is correct and much easier to reason about, but it allocates a new list at every node, which turns O(1) per node into O(n).

Pruning: not exploring what cannot work

So far the search visits everything. Backtracking becomes powerful when you check the constraints during the descent and abandon a branch the moment it cannot lead to a solution. Every node you skip takes its whole subtree with it.

Take a concrete problem: given twelve parcel weights, find every subset weighing exactly 60 kg. The plain version decides take-or-leave for each parcel and checks the total at the leaf. The pruned version stops at any node where the running total has already passed 60 (weights are positive, so it can only grow) or where everything still undecided adds up to less than what is still needed.

The diagram uses four weights so it fits on a page. [5, 4, 3, 2] with target 6 has a full tree of 31 nodes; the five abandoned branches shown remove 16 of them, leaving 15 visited, and the single solution [4, 2] is still found. Here are both versions on the twelve-parcel problem, counting nodes:

PARCELS = [31, 27, 15, 42, 8, 19, 23, 11, 37, 5, 14, 29]


def subsets_summing_to(weights: list[int], target: int) -> tuple[list[list[int]], int]:
    """Every subset that sums to target, exploring all 2**n leaves."""
    found: list[list[int]] = []
    nodes = 0

    def explore(index: int, path: list[int], total: int) -> None:
        nonlocal nodes
        nodes += 1
        if index == len(weights):
            if total == target:
                found.append(list(path))
            return
        path.append(weights[index])
        explore(index + 1, path, total + weights[index])
        path.pop()
        explore(index + 1, path, total)

    explore(0, [], 0)
    return found, nodes


def subsets_summing_to_pruned(weights: list[int], target: int) -> tuple[list[list[int]], int]:
    """The same answers, abandoning branches that provably cannot work."""
    # remaining[index] is the sum of everything from index onwards: the most
    # that can still be added once the first index parcels are decided.
    remaining = [0] * (len(weights) + 1)
    for index in range(len(weights) - 1, -1, -1):
        remaining[index] = remaining[index + 1] + weights[index]

    found: list[list[int]] = []
    nodes = 0

    def explore(index: int, path: list[int], total: int) -> None:
        nonlocal nodes
        nodes += 1
        if total > target:  # weights are positive, so the total only grows
            return
        if total + remaining[index] < target:  # not enough left to reach target
            return
        if index == len(weights):
            if total == target:
                found.append(list(path))
            return
        path.append(weights[index])
        explore(index + 1, path, total + weights[index])
        path.pop()
        explore(index + 1, path, total)

    explore(0, [], 0)
    return found, nodes


for target in (60, 25):
    plain, plain_nodes = subsets_summing_to(PARCELS, target)
    pruned, pruned_nodes = subsets_summing_to_pruned(PARCELS, target)
    print(f"target {target}: nodes {plain_nodes} -> {pruned_nodes} with pruning, "
          f"solutions {len(plain)}, identical answers: {plain == pruned}")

found, nodes = subsets_summing_to_pruned([5, 4, 3, 2], 6)
print(f"the diagram's tree: {found} found in {nodes} nodes instead of 31")

target 60: nodes 8191 -> 929 with pruning, solutions 9, identical answers: True
target 25: nodes 8191 -> 141 with pruning, solutions 1, identical answers: True
the diagram's tree: [[4, 2]] found in 15 nodes instead of 31

Same answers, 8,191 nodes down to 929. Target 25 does better still, at 141 nodes, because the first parcel weighs 31: taking it is already over target, so the take-it half of the tree collapses to the single node that gets rejected at depth 1 and its 4,094 descendants are never visited. Pruning pays most when it fires near the root, where the subtrees it discards are largest.

Be honest about what pruning buys, though. It does not change the growth rate. There are still 2ⁿ subsets, and an input where nothing can be ruled out early walks the whole tree. Pruning turns "impossible" into "fast enough for the inputs I actually have", which is usually what you need and never a guarantee.

Complexity

For any recursion, two independent quantities:

Time = number of calls × work per call. Count the calls with the recurrence the function itself describes.

factorial(n) makes exactly n calls, each doing one multiplication: O(n).
fib(n) makes 2 × fib(n + 1) − 1 calls. Fibonacci numbers grow like φⁿ where φ ≈ 1.618, so this is O(1.618ⁿ), loosely quoted as O(2ⁿ). The measured numbers above confirm it: 21,891 calls at n = 20 and 2,692,537 at n = 30, a factor of 123 for ten extra terms, and 1.618¹⁰ ≈ 123.
subsets(items) visits 2ⁿ⁺¹ − 1 nodes — a full binary tree of n + 1 levels — and copies an average of n/2 elements at each of its 2ⁿ leaves. The copying dominates: O(n · 2ⁿ). You cannot do better, because the output itself contains that many numbers.
permutations(items) produces n! results and scans all n candidates at each of the n levels: O(n · n!).

Space = maximum stack depth × frame size, plus the shared state. This is where beginners overestimate the cost.

factorial(n) and sum_to(n) reach depth n: O(n) stack. That is the reason both hit RecursionError.
fib(n) reaches depth n despite making exponentially many calls: O(n) stack.
subsets and permutations reach depth n and hold one path of at most n elements: O(n) working space, plus whatever the collected results occupy.

For divide-and-conquer recursions that halve the input — binary search, merge sort — the depth is log₂ n, which is 20 for a million items and 30 for a billion. That is why those algorithms recurse safely in Python and a linked-list walk does not.

When to use it, and when not to

Use recursion when the problem is defined in terms of itself and the depth is bounded. Trees, nested structures like JSON, grammars, divide and conquer. A recursive tree traversal is six lines and obviously correct; the iterative one needs a stack, a visited marker and careful ordering.

Use backtracking when you must search combinations under constraints and there is no formula, no greedy rule and no polynomial algorithm. Sudoku, N-Queens, crossword filling, timetabling, exact set cover.

Do not use recursion when the depth scales with input size past a few hundred; convert it to a loop or keep an explicit list as your stack. Do not use plain recursion when subproblems repeat — add functools.lru_cache or rewrite it as dynamic programming, which is the difference between 2,692,537 calls and 31 distinct subproblems.

Do not use backtracking when a cheaper structure exists. Subset-sum over small integers is a dynamic programming table. Shortest paths are Dijkstra. Permutations are itertools. Backtracking should be a considered choice, not a reflex.

Where it shows up in the real world

Parsers. CPython's own parser is recursive descent: PEP 617 replaced the old LL(1) grammar with a PEG parser in Python 3.9, and PEG parsing backtracks by design, trying alternatives in order and rewinding the input position when one fails. The standard library's JSON decoder is recursive too, in both its C and its pure-Python scanner, which is why json.loads on a deeply nested document raises RecursionError.

Regular expressions. Python's re module is a backtracking engine. a* first grabs as much as it can, and if the rest of the pattern then fails, it gives a character back and retries — choose, explore, un-choose, over the characters of a string. The same idea in twelve lines, matching glob patterns:

def matches(pattern: str, text: str) -> bool:
    """Glob-style matching where '?' is any one character and '*' is any run."""
    if not pattern:
        return not text
    if pattern[0] == "*":
        # Let the star match nothing; if the rest of the pattern then fails,
        # backtrack and hand the star one more character.
        return matches(pattern[1:], text) or (bool(text) and matches(pattern, text[1:]))
    if not text:
        return False
    if pattern[0] in ("?", text[0]):
        return matches(pattern[1:], text[1:])
    return False


for pattern, text in [("*.py", "main.py"), ("*.py", "main.txt"),
                      ("a*c", "abbbc"), ("a?c", "ac"), ("*x*y*", "axolotly")]:
    print(f"{pattern:6} vs {text:9} -> {matches(pattern, text)}")

*.py   vs main.py   -> True
*.py   vs main.txt  -> False
a*c    vs abbbc     -> True
a?c    vs ac        -> False
*x*y*  vs axolotly  -> True

The or is the backtrack: try the shorter match, and only if the whole rest of the pattern fails does the star swallow another character. The standard library equivalent is fnmatch.fnmatch, which works by translating the glob into a regex. The cost of that same mechanism is real: in 2016 Stack Overflow was taken offline for 34 minutes by catastrophic backtracking in a single regular expression scanning a post with a very long run of whitespace. Backtracking engines can go exponential on innocent-looking patterns.

Constraint satisfaction and SAT solvers. The DPLL algorithm from 1962, still the skeleton inside modern SAT solvers, is backtracking search over variable assignments with propagation as its pruning rule. Prolog's entire execution model is backtracking. Sudoku solvers, exam timetablers and puzzle generators are all the same shape: choose a value for the next empty slot, check the constraints, explore, undo.

Nested data. A directory tree, an HTML document and a JSON object are recursive structures, so the code that walks them is recursive too — shutil.rmtree descends into each subdirectory by calling itself, and every tree algorithm in this series does the same.

Common mistakes

No base case, or an unreachable one. if n == 1 misses n = 0 and every negative number. Test the boundary: factorial(0), subsets([]), an empty tree.

Recursing on the same size. factorial(n) inside factorial compiles fine and dies with RecursionError. Every call must strictly shrink the problem.

Storing the shared path instead of a copy. result.append(path) gives you a list of identical, mutated lists. result.append(list(path)) is the fix.

Forgetting to undo part of the state. Popping path but leaving used[index] = True marks every element consumed forever, so the search never backs up past its first descent — permutations([1, 2, 3]) returns [[1, 2, 3]] and nothing else, one ordering instead of six. Undo everything you did, in the reverse order you did it.

Assuming recursion is why your code is slow. Usually it is recomputation. Add functools.lru_cache before rewriting anything as a loop, and measure.

Returning nothing from the recursive branch. Writing factorial(n - 1) instead of return n * factorial(n - 1) makes the function return None — a bug that Python reports far away from where it happened.

Practice

Write a recursive function that reverses a string, with the empty string as the base case.
Write a recursive count_leaves for a nested list such as [1, [2, [3, 4]], 5], then state its maximum stack depth.
Convert factorial into a loop, and confirm both agree for every n from 0 to 20.
Generate all subsets of size exactly k by pruning any branch whose path is already longer than k.
Write a Sudoku solver: find the first empty cell, try digits 1 to 9, keep any digit that breaks no row, column or box constraint, recurse, and undo it if the recursion fails.

Summary

Recursion is two parts and one data structure: a base case, a strictly smaller recursive case, and the call stack that remembers everything half-finished. Backtracking adds one line — undo the last choice — and turns that stack into a systematic search over every arrangement, with pruning as the lever that makes it practical. Count calls for time, count depth for space, and remember that Python gives you 1000 frames and no tail-call optimisation, so depth is a design constraint rather than a detail.


Difficulty	Medium
Time, linear recursion	O(n) — one call per level, as in `factorial`
Time, branching recursion	O(bᵈ) — b calls per level, d levels deep, so `fib` is O(1.618ⁿ)
Time, full enumeration	O(n · 2ⁿ) for all subsets, O(n · n!) for all permutations
Space	O(depth) — one frame per open call, plus one shared path
Python depth limit	1000 frames by default, from `sys.getrecursionlimit()`
Tail-call optimisation	None, and deliberately so — it would destroy tracebacks
Data structure	The call stack, plus one mutable `path` list you undo
Use it when	The problem is self-similar, or you must enumerate under constraints
Avoid it when	Depth grows with input size, or subproblems repeat — loop or memoize
Real-world use	Recursive descent and PEG parsers, regex engines, SAT and CSP solvers
Python equivalents	`itertools.permutations`, `functools.lru_cache`, `re`, `fnmatch`

Keep reading

Dynamic Programming — the fix for the exponential fib above, in one decorator or one table.
The N-Queens Problem — backtracking with real constraints, and the clearest pruning you will see.
Depth-First Search — the same descent applied to a graph, both recursively and with an explicit stack.
Stacks — what the call stack actually is, and how to build your own when recursion runs out of frames.
Merge Sort — divide and conquer, where recursion depth is log n and the combining step is the algorithm.

Huffman Coding in Python: How Compression Actually Compresses

Bimal Kshetri — Sun, 16 Aug 2026 16:36:18 +0000

A plain text file spends eight bits on the letter e and eight bits on the letter q, even though English uses e about a hundred times more often. That is the waste Huffman coding removes. It hands every symbol a bit pattern whose length depends on how often the symbol appears — one bit for the workhorses, ten or twelve for the oddities — and the result is provably the best you can do with whole numbers of bits.

"Provably the best" is a strong claim and Huffman earns it. Among all codes that give each symbol a whole number of bits and can be decoded without separators between symbols, none produces a shorter output than Huffman's on the same frequency table. The construction that achieves this is a loop with two heap pops and one heap push in it.

It is also still running on your machine right now. Huffman codes are the entropy stage of DEFLATE, so they are inside nearly every .zip archive, nearly every Content-Encoding: gzip response and nearly every PNG image — nearly, because DEFLATE can also emit uncompressed blocks and ZIP permits other methods entirely. They are inside baseline JPEG and inside MP3 too. This post builds a working coder on abracadabra, decodes a bit string by hand, measures the real compression ratio including the table you have to ship alongside the bits, and is honest about where modern codecs have moved past it.

The idea

Why variable-length codes need a rule

Start with the boring option. abracadabra uses five distinct characters, so you could number them 0 to 4 and spend three bits on each: 11 characters, 33 bits. Decoding is trivial — chop the stream into three-bit pieces.

Now try to do better. a appears five times out of eleven, so give it a short code and the rare c and d long ones. The moment code lengths differ, decoding stops being obvious. Assign E the code 0, T the code 1 and A the code 01, and the received bits 01 mean either E then T, or a single A. No amount of cleverness in the decoder fixes that; the information is genuinely not there.

The rule that saves you is prefix-free: no symbol's code may be a prefix of another symbol's code. 01 is broken because 0 is a prefix of it. With that rule, decoding becomes a single left-to-right scan with no lookahead. Read bits until what you hold matches some code, and that match is forced — a longer code that started the same way would have your bits as a prefix, which the rule forbids.

Codes are paths down a binary tree

There is a neat way to guarantee prefix-freeness by construction rather than by checking. Build a binary tree. Label every left edge 0 and every right edge 1. Put the symbols only at the leaves. A symbol's code is the sequence of edge labels from the root down to its leaf.

Prefix-freeness now comes for free. Code X is a prefix of code Y exactly when the path to X is the start of the path to Y — which would mean X sits on the way to Y, and therefore is not a leaf. Symbols live at leaves, so it cannot happen.

Decoding is now a walk. Put a finger on the root. For each bit, step left on 0 and right on 1. When the finger lands on a leaf, emit that symbol and jump back to the root. Here is the bit string 0101101110 decoded against the tree above, one line per symbol:

bits     path from the root                       lands on   emit
0        left                                     leaf       E
1 0      right, left                              leaf       T
1 1 0    right, right, left                       leaf       A
1 1 1    right, right, right                      leaf       O
0        left                                     leaf       E

Ten bits in, ETAOE out. The decoder never looked ahead, never backtracked and never needed a separator or a length field. That is the entire payoff of the prefix-free property.

Which tree, though?

Any tree with the symbols at the leaves gives a valid code. The one you want is the cheapest. The total number of bits for a message is

cost(tree) = sum over symbols of (frequency of symbol) x (depth of its leaf)

because a symbol at depth d costs d bits every time it appears. Minimising that means pushing frequent symbols towards the root and rare ones towards the bottom.

Huffman's rule for building that tree is one sentence: take the two lowest-weight nodes, join them under a new parent whose weight is their sum, and put the parent back in the pool. Repeat until one node is left. Start with one leaf node per symbol, weighted by its frequency. Each merge removes two nodes and adds one, so n symbols take exactly n - 1 merges.

The reason this greedy rule is correct is easier to see with a second way of writing the same cost. Every merge you perform adds one edge above everything already inside the two nodes you merged, so it adds one bit to every occurrence of every symbol underneath. That means each merged node's weight gets billed exactly once:

cost(tree) = sum of the weights of all the internal nodes

Cheap merges early, expensive merges late — and that is what "always merge the two lightest" does.

Watching it work

Take the string abracadabra. Count the characters first:

Symbol	Frequency
`a`	5
`b`	2
`c`	1
`d`	1
`r`	2

Eleven characters, five distinct symbols. Fixed-width would cost 33 bits and ASCII costs 88.

Two of those frequencies tie at 2 and two tie at 1, so before any merging you have to decide what "the two lowest" means when several nodes are equal. Pick a tie-break and stick to it, or two runs of your own code can produce two different (equally good) trees. The rule used here: nodes enter the pool in alphabetical order and carry a strictly increasing counter, and ties are broken by the lower counter. So the initial pool, in the order it will be drained, is c=1, d=1, b=2, r=2, a=5.

Four merges, since there are five symbols:

Merge 1. The two lightest are c=1 and d=1. Join them under a node of weight 2. The pool is now b=2, r=2, (cd)=2, a=5 — the new node sorts last among the twos because it received the newest counter.

Merge 2. The two lightest are b=2 and r=2. Join them under a node of weight 4. Pool: (cd)=2, (br)=4, a=5.

Merge 3. The two lightest are (cd)=2 and (br)=4. Join them under a node of weight 6. Pool: a=5, (cdbr)=6.

Merge 4. Only two nodes are left, so they merge into the root, weight 11 — which is the length of the string, as it must be, because every character is counted exactly once somewhere below the root.

Now read the tree. Left edges are 0, right edges are 1:

Symbol	Frequency	Code	Bits used
`a`	5	`0`	5
`b`	2	`110`	6
`c`	1	`100`	3
`d`	1	`101`	3
`r`	2	`111`	6

Total: 23 bits, against 33 for fixed width and 88 for ASCII. Check it against the other cost formula — the internal nodes weigh 2, 4, 6 and 11, and 2 + 4 + 6 + 11 = 23. The two ways of counting agree, as they must.

Encoding abracadabra is now a table lookup per character:

a   b    r    a   c    a   d    a   b    r    a
0   110  111  0   100  0   101  0   110  111  0

which runs together as 01101110100010101101110. Feed that back into the decoding walk and you get abracadabra — the first 0 lands immediately on the a leaf, then 110 steps right, right, left to b, and so on.

Why greedy is right here

Greedy algorithms are usually wrong, so the argument matters. It has two halves.

First, there is an optimal tree in which the two rarest symbols are siblings at the deepest level. Take any optimal tree. No internal node in it has only one child, because deleting such a node would shorten a code and cut the cost. So the deepest leaf has a sibling, and that sibling must be a leaf at the same depth — anything hanging below it would be deeper still. Call those two slots p and q, and swap the rarest symbol into p and the second-rarest into q. Each swap moves a lower frequency to a greater-or-equal depth and a higher frequency to a lesser-or-equal depth, so the total cost cannot go up. The tree was optimal, so the swapped tree is optimal too, and in it the two rarest symbols are siblings at the bottom.

Second, merging them leaves a smaller version of the same problem. Replace those two siblings with a single symbol whose frequency is their sum. Any tree for the smaller alphabet turns into a tree for the original by splitting that leaf in two, and the cost goes up by exactly the merged weight either way. So an optimal tree for the smaller problem gives an optimal tree for the bigger one.

Put them together and induct on the number of symbols: the first merge is safe, and what remains is the same problem one symbol shorter. That is a textbook exchange argument, and it is what separates Huffman from the greedy rules that merely look sensible.

The code

Counting first. collections.Counter is the standard-library tool for this and there is no reason to write your own tally loop.

from __future__ import annotations

from collections import Counter

TEXT = "abracadabra"
counts = Counter(TEXT)

for symbol, weight in sorted(counts.items()):
    print(f"{symbol}  {weight}")

print(f"{len(TEXT)} characters, {len(counts)} distinct symbols")
print(f"three bits each (fixed width): {3 * len(TEXT)} bits")
print(f"eight bits each (ASCII):       {8 * len(TEXT)} bits")

a  5
b  2
c  1
d  1
r  2
11 characters, 5 distinct symbols
three bits each (fixed width): 33 bits
eight bits each (ASCII):       88 bits

Now the tree. The pool of nodes waiting to be merged needs one operation — "give me the lightest" — repeated 2(n - 1) times, which is exactly what a min-heap is for. Python's is heapq, and it works on a plain list.

import heapq
from dataclasses import dataclass
from typing import Optional


@dataclass
class Node:
    """One node of the code tree. A leaf carries a symbol; an internal node
    carries None and exists only to hold two children together."""

    symbol: Optional[str]
    weight: int
    left: Optional["Node"] = None
    right: Optional["Node"] = None

    @property
    def is_leaf(self) -> bool:
        return self.left is None and self.right is None


def build_tree(counts: dict[str, int], trace: bool = False) -> Optional[Node]:
    """Merge the two lightest nodes until one node is left, and return it.

    The heap holds `(weight, tiebreak, node)`. The tiebreak counter is not
    decoration: it fixes the order of equal weights so the output is
    reproducible, and it guarantees heapq never has to compare two Nodes.
    """
    if not counts:
        return None

    tiebreak = 0
    heap: list[tuple[int, int, Node]] = []
    for symbol, weight in sorted(counts.items()):
        heap.append((weight, tiebreak, Node(symbol, weight)))
        tiebreak += 1
    heapq.heapify(heap)

    while len(heap) > 1:
        left_weight, _, left = heapq.heappop(heap)
        right_weight, _, right = heapq.heappop(heap)
        parent = Node(None, left_weight + right_weight, left, right)
        if trace:
            print(f"merge {name(left)} + {name(right)} -> {parent.weight}"
                  f"   pool: {pool(heap, parent, tiebreak)}")
        heapq.heappush(heap, (parent.weight, tiebreak, parent))
        tiebreak += 1

    return heap[0][2]


def name(node: Node) -> str:
    """A readable label: the symbol for a leaf, its leaves for a merged node."""
    if node.is_leaf:
        return f"{node.symbol}={node.weight}"
    return f"({''.join(symbols_of(node))})={node.weight}"


def symbols_of(node: Node) -> list[str]:
    if node.is_leaf:
        return [node.symbol]
    return symbols_of(node.left) + symbols_of(node.right)


def pool(heap: list[tuple[int, int, Node]], pending: Node, tiebreak: int) -> str:
    """The nodes still waiting to be merged, in the order the heap will pop."""
    waiting = sorted(heap + [(pending.weight, tiebreak, pending)],
                     key=lambda item: item[:2])
    return "  ".join(name(item[2]) for item in waiting)


root = build_tree(dict(counts), trace=True)
print(f"root weight: {root.weight}")

merge c=1 + d=1 -> 2   pool: b=2  r=2  (cd)=2  a=5
merge b=2 + r=2 -> 4   pool: (cd)=2  (br)=4  a=5
merge (cd)=2 + (br)=4 -> 6   pool: a=5  (cdbr)=6
merge a=5 + (cdbr)=6 -> 11   pool: (acdbr)=11
root weight: 11

Four merges, matching the hand trace line for line.

Turning the tree into a lookup table is one depth-first walk that carries the path so far:

def build_codes(root: Optional[Node]) -> dict[str, str]:
    """Map each symbol to its bit string. A left edge is 0, a right edge is 1."""
    if root is None:
        return {}
    if root.is_leaf:
        # A one-symbol alphabet has no edges at all, so the walk below would
        # hand it the empty code and nothing could ever be decoded.
        return {root.symbol: "0"}

    codes: dict[str, str] = {}

    def walk(node: Node, prefix: str) -> None:
        if node.is_leaf:
            codes[node.symbol] = prefix
            return
        walk(node.left, prefix + "0")
        walk(node.right, prefix + "1")

    walk(root, "")
    return codes


codes = build_codes(root)
total_bits = 0
for symbol, code in sorted(codes.items(), key=lambda pair: (len(pair[1]), pair[0])):
    cost = len(code) * counts[symbol]
    total_bits += cost
    print(f"{symbol}  {code:<4} x{counts[symbol]}  = {cost:>2} bits")
print(f"total: {total_bits} bits")

a  0    x5  =  5 bits
b  110  x2  =  6 bits
c  100  x1  =  3 bits
d  101  x1  =  3 bits
r  111  x2  =  6 bits
total: 23 bits

Encoding is a join over table lookups. Decoding is the finger-on-the-tree walk, written as a loop:

def encode(text: str, codes: dict[str, str]) -> str:
    return "".join(codes[symbol] for symbol in text)


def decode(bits: str, root: Optional[Node]) -> str:
    """Take one edge per bit. Landing on a leaf emits a symbol and restarts."""
    if root is None or not bits:
        return ""
    if root.is_leaf:
        return root.symbol * len(bits)

    out: list[str] = []
    node = root
    for bit in bits:
        node = node.left if bit == "0" else node.right
        if node.is_leaf:
            out.append(node.symbol)
            node = root

    if node is not root:
        raise ValueError("bit string ran out part-way down the tree")
    return "".join(out)


bits = encode(TEXT, codes)
print(bits)
print(f"{len(bits)} bits, against {8 * len(TEXT)} bits of ASCII")
print(decode(bits, root))
print(decode(bits, root) == TEXT)

01101110100010101101110
23 bits, against 88 bits of ASCII
abracadabra
True

How the code maps to the idea

The heap tuple is (weight, tiebreak, node) and every part earns its place. weight is the key the algorithm cares about. tiebreak increments on every push, so no two tuples are ever equal in their first two fields — which makes the output reproducible across runs and, less obviously, stops Python from ever reaching the third element during a comparison. Drop it and the first frequency tie crashes with a TypeError, because Node has no ordering.

heapq.heapify seeds the heap in O(n), not O(n log n), which is why the leaves are appended to a plain list first rather than pushed one at a time. It does not change the overall bound, but it is free.

The loop condition is len(heap) > 1, not a fixed range. Each pass pops two and pushes one, so the heap shrinks by exactly one and ends holding the root. It also handles a single-symbol alphabet for nothing: the body never runs and the lone leaf is returned.

build_codes carries the path in an argument, not a shared list. prefix + "0" gives each branch its own string, so there is nothing to undo on the way back up — a shared list you append to and forget to pop is the classic way to get this function subtly wrong. The single-symbol case still needs its own line, though: a tree of one leaf and no edges would hand aaaa the empty code and encode it to zero bits, so it gets the one-bit code 0 instead.

decode ends by checking that the finger came back to the root. A bit string that stops half-way down a branch is truncated or corrupt, and saying so beats silently dropping the last symbol.

The header you cannot skip

Here is the part that tutorials tend to leave out. The 23 bits above are meaningless on their own. The decoder needs the same tree the encoder used, and it has no way to derive it — the frequencies came from the original message, which is precisely what the decoder does not have. So a real Huffman file is header plus payload, and the header is not free.

You have two reasonable ways to write one. Storing raw frequencies means a symbol plus a count per entry, and the count needs four bytes to hold a large file's tally — five bytes an entry. Storing code lengths is much cheaper: given only "a is 1 bit, b is 3 bits, c is 3 bits, ..." you can rebuild a code with exactly those lengths by handing out bit patterns in a fixed order — shortest length first, symbols in order within a length. That is canonical Huffman, it is what DEFLATE does, and it costs two bytes an entry rather than five. The codes come out different from the ones the tree gave you, but the lengths are the same, so the payload is the same size.

The measurement below uses that cheaper header: a two-byte count of entries, then one byte of symbol and one byte of code length each.

import math


def measure(text: str) -> tuple[int, int, int, int]:
    """Return (distinct symbols, payload bits, header bits, stored bytes).

    The header is the honest part. Storing one byte of symbol and one byte of
    code length per entry, plus a two-byte entry count, is enough for the
    decoder to rebuild an identical canonical tree.
    """
    counts = dict(Counter(text))
    table = build_codes(build_tree(counts))
    payload = sum(len(table[symbol]) * count for symbol, count in counts.items())
    header = 16 + 16 * len(counts)
    return len(counts), payload, header, math.ceil((payload + header) / 8)


SAMPLE = (
    "Compression works because real data is lumpy. The letter e turns up in "
    "English text roughly a hundred times more often than the letter q, and yet "
    "a plain text file spends exactly eight bits on each of them. Huffman coding "
    "removes that waste by giving common symbols short codes and rare symbols "
    "long ones, and it does so optimally: no other code that spends a whole "
    "number of bits per symbol can beat it on the same frequency table."
)

print(f"{'input':<22}{'syms':>5}{'payload':>9}{'header':>8}{'stored':>8}{'ratio':>8}")
for label, text in [
    ("abracadabra", TEXT),
    ("sample paragraph", SAMPLE),
    ("sample x 4", SAMPLE * 4),
    ("sample x 40", SAMPLE * 40),
]:
    distinct, payload, header, stored = measure(text)
    ratio = stored / len(text)
    print(f"{label:<22}{distinct:>5}{payload:>9}{header:>8}{stored:>8}{ratio:>7.2f}x")

input                  syms  payload  header  stored   ratio
abracadabra               5       23      96      15   1.36x
sample paragraph         32     1843     528     297   0.69x
sample x 4               32     7372     528     988   0.57x
sample x 40              32    73720     528    9281   0.54x

Read the first row carefully: abracadabra compresses from 11 bytes to 15 bytes. Huffman made it bigger. The payload did shrink, from 88 bits to 23, but 96 bits of header swamped the saving. That is arithmetic, not a bug — header cost is proportional to alphabet size while the saving is proportional to message length, so a message loses when its alphabet is wide relative to its length. The ratio decides, not the length alone: eleven characters over five symbols stores as 15 bytes, but eleven characters over two symbols carries a six-byte header and stores as 8, a real win. DEFLATE knows this and offers a "stored" block type that copies bytes verbatim, used whenever compression would expand the data.

Repeat the same distribution down the table and the header amortises away: 0.69x, then 0.57x, then 0.54x, converging on the payload-only ratio of 1843 / (432 × 8) = 0.53.

The other number worth knowing is the floor. Shannon's source coding theorem says no code that assigns bits to symbols independently can beat the entropy of the distribution, and Huffman always lands within one bit per symbol of it — usually far closer.

def entropy_bits(text: str) -> float:
    """Shannon's lower bound: no per-symbol code can beat this on this text."""
    counts = Counter(text)
    return -sum(c * math.log2(c / len(text)) for c in counts.values())


for label, text in [("abracadabra", TEXT), ("sample paragraph", SAMPLE)]:
    _, payload, _, _ = measure(text)
    floor = entropy_bits(text)
    print(f"{label:<20} huffman {payload / len(text):.3f} bits/symbol, "
          f"entropy {floor / len(text):.3f}, overhead "
          f"{(payload - floor) / len(text):.3f}")

abracadabra          huffman 2.091 bits/symbol, entropy 2.040, overhead 0.051
sample paragraph     huffman 4.266 bits/symbol, entropy 4.224, overhead 0.042

Five hundredths of a bit per symbol off the theoretical floor. The gap exists because Huffman must round every code to a whole number of bits, and it is the one thing arithmetic coding fixes.

Complexity

Let n be the number of distinct symbols and m the length of the message. Keeping those two separate is essential — for English text n is a few dozen while m can be gigabytes.

Counting frequencies: O(m). One pass, one dictionary update per character.

Building the tree: O(n log n). Seeding costs O(n log n) as written, because the code sorts the symbols before it heapifies — that sort buys the alphabetical tie-break, not the heap, and heapify itself is O(n). The loop then runs exactly n - 1 times, because each pass consumes two nodes and produces one, taking the pool from n down to 1. Each pass does two heappop calls and one heappush, and every heap operation walks one root-to-leaf path of a heap holding at most n items, which is at most log₂ n steps. So the merging is 3(n - 1) heap operations at O(log n) each — that is the n and that is the log n. For 256 byte values, n - 1 is 255 merges and log₂ 256 is 8, so the whole tree build is around 6,000 elementary steps regardless of file size.

If the frequencies arrive already sorted you can drop the heap entirely and do it in O(n) with two FIFO queues: one holding the sorted leaves, one holding the merged nodes in the order they were created, which is automatically non-decreasing. The next-lightest node is always at the front of one of the two queues. That is a real technique, not a curiosity — it is why sorting-based Huffman implementations exist.

Building the code table: O(total code length). The walk visits each of the 2n - 1 nodes once, but the strings it builds cost the sum of all code lengths. That sum is n × (average depth), and the worst case is worse than it looks: frequencies that grow like the Fibonacci sequence produce a completely lopsided tree with a leaf at depth n - 1, giving O(n²) characters of code strings. In practice depth stays small, and there is a hard ceiling — to force a code of length d you need each level to be at least as heavy as the sum of the two below it, so total frequency must grow at least as fast as the Fibonacci numbers. A message of fewer than 4.8 billion symbols can never produce a code longer than 44 bits. DEFLATE tightens that to 15 bits by construction.

Encoding: O(m + B), where B is the number of output bits. One table lookup per input character plus the cost of writing the bits.

Decoding: O(B). Exactly one tree step per bit — no searching, no backtracking. Real decoders speed this up with a lookup table that consumes several bits at once, trading memory for a shorter walk.

Space: O(n) for the tree, plus the code table, plus O(B) for the output. The tree has 2n - 1 nodes, since a binary tree with n leaves and no one-child nodes has n - 1 internal nodes. The code table is the part that is easy to under-count: it holds every code as a string, so it costs the same sum of code lengths that building it cost — O(n²) characters on the Fibonacci worst case above, not O(n). build_codes is also recursive, so it holds one stack frame per level, up to n - 1 of them, which is why a deep enough tree hits Python's recursion limit rather than merely running slowly.

The n log n here is cheap in absolute terms because n is bounded by your alphabet, not your data. Compressing a 1 GB file with a byte alphabet still builds the tree in 255 merges. The linear O(m) passes over the data are what actually dominate the clock.

When to use it, and when not to

Use it when you control a format and the symbol frequencies are genuinely skewed. Huffman decoding is fast, branch-light and easy to make correct, which is why it survives inside formats that need to decode at video frame rates.

Use it as the last stage of a pipeline, not as the whole pipeline. This is the important limitation. Huffman models only how often each symbol appears; it is blind to order. A file that is ab repeated a million times has two symbols at 50% each, so Huffman spends one bit on each and saves nothing over the fixed-width baseline — even though a human can describe that file in one sentence. Repetition is caught by a dictionary method like LZ77, which replaces repeats with back-references. DEFLATE is exactly that pairing: LZ77 first, Huffman second on whatever skew is left.

Do not reach for it when a general-purpose compressor will do. If you just want a file to be smaller, use zlib, gzip or lzma from the standard library — C implementations of complete formats that will beat a hand-rolled symbol coder on essentially any real input.

Do not use it on short messages. As the table above shows, an eleven-byte input with five distinct symbols comes back as fifteen bytes: the payload fell from 88 bits to 23, saving eight bytes, and the twelve-byte header more than swallowed them. If you must compress small records, share one pre-agreed table across all of them so the header is paid once — which is what HTTP/2's HPACK does with its static Huffman table.

Do not use it when you need the last few percent. A symbol with probability 0.9 deserves 0.15 bits; Huffman must spend a whole one. Arithmetic coding, range coding and the newer asymmetric numeral systems all encode at fractional bit cost and close that gap.

Where it shows up in the real world

DEFLATE (RFC 1951) is the big one, and it is everywhere: .zip archives, gzip, the zlib library, PNG's image data, and HTTP responses sent with Content-Encoding: gzip. Each DEFLATE block runs LZ77 and then Huffman-codes the resulting literals, match lengths and distances. A block either uses a fixed table baked into the spec or carries its own dynamic one — and that dynamic table is stored as canonical code lengths which are themselves Huffman-coded, a second layer of the same trick.

Baseline JPEG Huffman-codes the quantised DCT coefficients, with separate tables for the DC and AC coefficients of the luminance and chrominance channels. JPEG also defines an arithmetic-coding mode that compresses a few percent better, but patent worries in the 1990s meant almost nothing implemented it, and baseline Huffman became the version everyone actually ships.

MP3 (MPEG-1 Audio Layer III) picks from a set of predefined Huffman tables to code the quantised frequency-domain values in each granule. Brotli, the algorithm behind Content-Encoding: br, uses Huffman codes alongside a large built-in dictionary and context modelling.

Zstandard uses Huffman for literals but switches to Finite State Entropy, a tANS coder, for the sequence data — a deliberate split, because tANS handles skewed distributions at fractional bit cost while Huffman decodes literals faster.

Be accurate about the trend, though. The highest-compression modern codecs have largely moved past Huffman for their main entropy stage: H.264 and H.265 use CABAC, a context-adaptive binary arithmetic coder, and AV1 uses a multi-symbol arithmetic coder. Huffman persists where decode speed and simplicity matter more than the last few percent — which is still an enormous amount of the world's data.

Common mistakes

Pushing bare nodes into the heap. Without a tie-break value, the first frequency tie makes heapq compare two Node objects:

broken_heap = [(1, Node("c", 1))]
try:
    heapq.heappush(broken_heap, (1, Node("d", 1)))
except TypeError as error:
    print(f"TypeError: {error}")

TypeError: '<' not supported between instances of 'Node' and 'Node'

The fix is the counter in the middle of the tuple. @dataclass(order=True) looks like the shortcut and is not one: it compares fields in declaration order, so the first field it reaches is symbol — a string on a leaf, None on an internal node. Leaves and merged nodes tie all the time, and the moment they do you swap one crash for another:

@dataclass(order=True)
class OrderedNode:
    symbol: Optional[str]
    weight: int
    left: Optional["OrderedNode"] = None
    right: Optional["OrderedNode"] = None


leaf_b = OrderedNode("b", 2)
merged_cd = OrderedNode(None, 2, OrderedNode("c", 1), OrderedNode("d", 1))
try:
    heapq.heappush([(2, leaf_b)], (2, merged_cd))
except TypeError as error:
    print(f"TypeError: {error}")

TypeError: '<' not supported between instances of 'NoneType' and 'str'

That is the pool after merge 1 of abracadabra, where (cd)=2 meets b=2 — ordered nodes do not survive even the worked example. The counter sidesteps the question entirely: two tuples never tie on their first two fields, so the Node is never compared.

Forgetting to store the table. An encoder that returns only the bit string has produced something nobody can decode, including itself after a restart. Ship the header or agree the table in advance; there is no third option.

Non-deterministic tie-breaking. Iterating a dictionary you built in arbitrary order, or breaking ties by object identity, gives you a different-but-equally-optimal tree each run. Every such tree is fine on its own, but if any part of your system rebuilds the tree separately it will disagree with the encoder, and the output will decode to garbage.

Putting a symbol on an internal node. It is tempting when an alphabet feels small. It destroys the prefix-free property and makes decoding ambiguous — the whole reason for the tree.

Treating a bit string as bytes. 23 bits is not a whole number of bytes. You must pad the last byte and store how many padding bits you added, or the decoder will walk a few extra edges at the end and emit a phantom symbol. Do not lean on the if node is not root check in decode for this: it only fires when the padding happens to strand the walk part-way down a branch. Pad abracadabra's 23 bits with one 0 and the extra step lands cleanly on the a leaf, so the check passes and you get abracadabraa back with no complaint at all.

Compressing already-compressed data. Running Huffman over a JPEG or a ZIP gives you a nearly uniform byte distribution, near-maximal entropy, and an output slightly larger than the input thanks to the header.

Practice

Add byte packing: convert the bit string into a bytes object with the padding length recorded, and confirm the round trip still returns the original text.
Write a serialise_header and parse_header pair that stores the code lengths, then rebuild the code table from lengths alone using canonical ordering — shortest length first, alphabetical within a length.
Find the break-even point: for the sample paragraph's distribution, how many characters must a message have before the compressed form is smaller than the original?
Replace the recursive build_codes walk with an explicit stack, so it cannot hit Python's recursion limit on a pathologically deep tree.
Build the Fibonacci worst case — frequencies 1, 1, 2, 3, 5, 8, 13, ... for 20 symbols — and check that the longest code really is 19 bits.

Summary

Huffman coding is the clearest example in the whole greedy family of a local rule that provably reaches a global optimum, and the exchange argument that proves it is short enough to hold in your head. Build a min-heap of leaves, merge the two lightest n - 1 times, read the codes off the tree. Then remember the two things the theory does not tell you: the header has to travel with the bits, and symbol frequencies alone cannot see repetition, which is why every serious format wraps Huffman around something else.


Difficulty	Hard
Build time	O(n log n) — `n - 1` merges, three heap operations each
Build time, sorted input	O(n) — two FIFO queues instead of a heap
Encode	O(m + B) — one table lookup per input symbol
Decode	O(B) — exactly one tree step per bit, no backtracking
Space	O(n) for the `2n - 1` tree nodes, plus the code table (the sum of all code lengths, O(n²) at worst), plus O(n) of recursion stack in `build_codes`, plus O(B) for the output
Optimal	Yes, among prefix codes with whole-bit symbol codes
Within	1 bit per symbol of the entropy floor, usually far less
Data structure	Binary tree built with a min-heap
Use it when	Frequencies are skewed and known, and decode speed matters
Avoid it when	The data is repetitive rather than skewed, or messages are tiny
Real-world use	DEFLATE (ZIP, gzip, PNG), baseline JPEG, MP3, Brotli, zstd literals
Python equivalent	`heapq` for the queue; `zlib` / `gzip` / `lzma` for actual compression

Keep reading

Heaps and Priority Queues — the structure doing all the work in the merge loop, built from scratch.
Greedy Algorithms — the family Huffman belongs to, and how to prove a greedy rule before trusting it.
Tries (Prefix Trees) — the other tree where the path spells the answer, used for autocomplete instead of compression.
The Coin Change Problem — a greedy rule that looks just as reasonable and is wrong, with the DP that fixes it.
Big O Notation — the counting arguments behind every bound quoted above.