13 — Shortest paths, spanning trees and beyond

Chapter 12 covered how to walk a graph: representations, depth-first and breadth-first search, topological order, union-find and cycle detection. This chapter builds on those walks. Most of the algorithms here are a traversal with one extra idea attached: a priority queue that always expands the cheapest frontier node, a relaxation loop repeated until nothing changes, a DFS that records the lowest discovery time reachable from each subtree, or a bitmask that remembers which vertices a path has already used.

The chapter is in four parts. Shortest paths come first (Dijkstra, Bellman-Ford, Floyd-Warshall, 0-1 BFS, A*), then minimum spanning trees (Prim and Kruskal), then the structural DFS algorithms (strongly connected components, bridges and articulation points, Eulerian paths), and finally the hard problems (Hamiltonian paths, the Traveling Salesman Problem and graph coloring), where exact algorithms are exponential and heuristics matter.

Every implementation below was checked against brute force on hundreds of random graphs before being included. Timings were measured with Python 3.14 on one Mac using time.perf_counter, taking the best of several runs. Exact numbers will differ elsewhere; the ratios between them are the point.

Notation: $V$ is the number of vertices, $E$ the number of edges, $w(u, v)$ the weight of edge $u \to v$, and $\delta(s, v)$ the true shortest-path distance from $s$ to $v$. Graphs are stored as adjacency lists unless stated otherwise.

Shortest paths: the shared idea

Every single-source shortest-path algorithm keeps a tentative distance $d[v]$ for each vertex, starting at $d[s] = 0$ and $d[v] = \infty$ elsewhere, and improves it with one operation, relaxation:

relax(u, v):  if d[u] + w(u, v) < d[v]:  d[v] = d[u] + w(u, v);  parent[v] = u

Relaxation never makes $d[v]$ smaller than the true distance, as long as $d$ only ever holds lengths of real paths. So $d[v] \ge \delta(s, v)$ always holds, and the algorithms differ only in which edges they relax, in what order, and how many times:

The parent pointers form a shortest-path tree: following them backwards from any vertex gives a shortest path to it. This works because of the optimal substructure of shortest paths: any prefix of a shortest path is itself a shortest path. If the prefix $s \leadsto u$ of a shortest path $s \leadsto u \to v$ could be shortened, the whole path could be too.

Dijkstra’s algorithm

The problem and the invariant

Dijkstra’s algorithm (Edsger Dijkstra, published 1959) finds the shortest distance from one source to every vertex in a graph with non-negative edge weights. It is BFS generalized to weights: BFS processes vertices in order of hop count, and Dijkstra processes them in order of distance.

The invariant: the algorithm keeps a set of finalized vertices whose distance is known to be correct. At each step it takes the unfinalized vertex with the smallest tentative distance, declares it final, and relaxes its outgoing edges.

Why the vertex with the smallest tentative distance $d[u]$ must already be correct: take any path from $s$ to $u$. It starts inside the finalized set and at some point leaves it, through a first vertex $x$ outside the set. The prefix up to $x$ costs at least $d[x]$, because $x$‘s predecessor was finalized and relaxed its edge to $x$. And $d[x] \ge d[u]$, because $u$ was chosen as the minimum. The rest of the path, from $x$ to $u$, costs at least zero, because no edge is negative. So every path to $u$ costs at least $d[u]$, and $d[u] = \delta(s, u)$.

That last step is the only place non-negativity is used, and it is exactly where negative edges break the algorithm.

Walk-through

        4
   A ------> B
   |       ^ |
 1 |     2/  | 1
   v     /   v
   C -------> D ------> E
        5          3

Edges: A→B 4, A→C 1, C→B 2, C→D 5, B→D 1, D→E 3. The heap holds (distance, vertex) pairs.

Implementation with heapq and lazy deletion

Python’s heapq is a binary min-heap with no decrease-key operation. The standard workaround is lazy deletion: when a vertex’s distance improves, push a new entry and leave the old one in the heap. When an entry is popped, check whether its vertex was already finalized and skip it if so. The entries skipped in steps 4 and 6 above are these stale entries.

import heapq

INF = float("inf")


def dijkstra(graph, source):
    """graph: dict node -> list of (neighbor, weight), weights >= 0."""
    dist = {source: 0}
    parent = {source: None}
    heap = [(0, source)]
    done = set()
    while heap:
        d, u = heapq.heappop(heap)
        if u in done:            # stale entry: u was already finalized
            continue
        done.add(u)
        for v, w in graph.get(u, ()):
            nd = d + w
            if nd < dist.get(v, INF):
                dist[v] = nd
                parent[v] = u
                heapq.heappush(heap, (nd, v))
    return dist, parent


def path_to(parent, target):
    if target not in parent:
        return None
    path = []
    while target is not None:
        path.append(target)
        target = parent[target]
    return path[::-1]

On the graph above:

graph = {
    "A": [("B", 4), ("C", 1)],
    "C": [("B", 2), ("D", 5)],
    "B": [("D", 1)],
    "D": [("E", 3)],
}
dist, parent = dijkstra(graph, "A")
print(dist)
print(path_to(parent, "E"))

{'A': 0, 'B': 3, 'C': 1, 'D': 4, 'E': 7}
['A', 'C', 'B', 'D', 'E']

An alternative to the done set is to compare the popped distance with the stored one, if d > dist[u]: continue. Both work for non-negative weights.

How many stale entries appear in practice? On a random directed graph with 10,000 vertices and 50,000 edges (weights 1 to 100):

reachable nodes 9923, pushes 14277, pops 14277, stale pops 4354

About 30% of pops were stale. The heap never holds more than $E + 1$ entries, so the cost of laziness is a constant factor, not a change in complexity.

Complexity

Space is $O(V + E)$ for the lazy version. If only one target matters, return as soon as the target is popped: its distance is final at that moment.

Why Dijkstra fails with negative edges

The proof above used “the rest of the path costs at least zero”. With a negative edge, a path can leave the finalized set through an expensive vertex and then come back cheaper.

        1           1
   S ------> A ------> T
   |         ^
 3 |        / -3
   v       /
   B -----+

Edges: S→A 1, S→B 3, B→A −3, A→T 1. The true distances are A = 0 (via B) and T = 1. Dijkstra pops A first at distance 1, because 1 < 3, declares it final, and relaxes A→T to 2. When B is popped at 3, the edge B→A offers 0, but A is already final. In the textbook form, finalized vertices are never updated:

def dijkstra_textbook(graph, source):
    # The classic form: once a node is popped it is final and never updated.
    dist = {source: 0}
    heap = [(0, source)]
    done = set()
    while heap:
        d, u = heapq.heappop(heap)
        if u in done:
            continue
        done.add(u)
        for v, w in graph.get(u, ()):
            if v not in done and d + w < dist.get(v, INF):
                dist[v] = d + w
                heapq.heappush(heap, (d + w, v))
    return dist

Dijkstra:     {'S': 0, 'A': 1, 'B': 3, 'T': 2}
Bellman-Ford: {'S': 0, 'A': 0, 'B': 3, 'T': 1}

The lazy version shown earlier fares no better. It lowers dist["A"] to 0 when B is relaxed, but A is in done so it is never expanded again, and T stays wrong: it returns {'S': 0, 'A': 0, 'B': 3, 'T': 2}, which is internally inconsistent.

Two common fixes do not work:

• Adding a constant to every edge to make weights non-negative changes which path is shortest. A path with $k$ edges gains $k$ times the constant, so paths with fewer edges are favored. In the example, S→A→T gains 2c while S→B→A→T gains 3c. The correct reweighting is Johnson’s: compute potentials $p(v)$ with one Bellman-Ford run and use $w'(u, v) = w(u, v) + p(u) - p(v)$, which preserves shortest paths because every $s \leadsto t$ path changes by the same $p(s) - p(t)$.
• Removing the done check so that improved vertices are re-expanded does give correct answers when there is no negative cycle, but the number of re-expansions can grow exponentially in the worst case. At that point Bellman-Ford’s $O(VE)$ is the better guarantee.

When to use, and pitfalls

Use Dijkstra for single-source shortest paths with non-negative weights: road networks, network routing (OSPF computes shortest paths with Dijkstra), latency graphs, grids with varying terrain cost.

• Zero-weight edges are fine. Only negative ones break the proof.
• Do not mark a vertex visited when it is pushed, only when it is popped. Marking on push is correct for BFS but wrong for Dijkstra: a later, cheaper path to the same vertex would be ignored.
• Tuples compare element by element. If two entries tie on distance, Python compares the vertices next. That fails for vertices that are not comparable (dicts, custom objects). Insert a counter from itertools.count() as a tie-breaker.
• Floating-point weights can make exact comparisons fragile when paths of equal length are computed in different orders; this matters for tie handling, not for correctness of the minimum.
• Dijkstra solves the single-pair problem as a by-product. For repeated single-pair queries on a large fixed graph (routing), production systems add preprocessing such as contraction hierarchies or use A* with landmarks.

Classic problem: Network Delay Time

A signal is sent from node k through a network of directed, timed edges. How long until every node has received it? The answer is the largest shortest-path distance from k, or −1 if some node is unreachable.

def network_delay_time(times, n, k):
    graph = defaultdict(list)
    for u, v, w in times:
        graph[u].append((v, w))
    dist = {}
    heap = [(0, k)]
    while heap:
        d, u = heapq.heappop(heap)
        if u in dist:
            continue
        dist[u] = d
        for v, w in graph[u]:
            if v not in dist:
                heapq.heappush(heap, (d + w, v))
    return max(dist.values()) if len(dist) == n else -1

This version records a distance only when a vertex is popped, which is a compact way to write lazy Dijkstra. network_delay_time([(2, 1, 1), (2, 3, 1), (3, 4, 1)], 4, 2) returns 2.

Bellman-Ford

The problem and the invariant

The Bellman-Ford algorithm, published independently by Lester Ford (1956) and Richard Bellman (1958), with Edward Moore describing the same idea in 1959, handles negative edge weights and detects negative cycles.

The idea is to relax every edge, and repeat. The invariant: after round $i$, $d[v]$ is at most the length of the shortest path from $s$ to $v$ that uses at most $i$ edges. This holds by induction. The shortest path with at most $i + 1$ edges is a shortest path with at most $i$ edges to some $u$, followed by the edge $u \to v$, and round $i + 1$ relaxes that edge.

If there is no negative cycle, every shortest path is simple and has at most $V - 1$ edges, so $V - 1$ rounds are enough. If a $V$-th round still improves some distance, a path with $V$ or more edges is shorter than every simple path, which is only possible if it goes around a negative cycle.

Implementation

def bellman_ford(n, edges, source):
    """n nodes 0..n-1, edges: list of (u, v, w). Returns (dist, has_neg_cycle)."""
    dist = [INF] * n
    dist[source] = 0
    for _ in range(n - 1):
        changed = False
        for u, v, w in edges:
            if dist[u] != INF and dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
                changed = True
        if not changed:
            break
    for u, v, w in edges:
        if dist[u] != INF and dist[u] + w < dist[v]:
            return dist, True
    return dist, False

The dist[u] != INF check matters. Without it, $\infty + (-3)$ is still $\infty$ for floats, but with integer sentinels such as 10**18 an unreachable vertex could “improve” its neighbors. The changed flag stops early once a full round makes no update.

This detects only negative cycles reachable from the source. To detect any negative cycle in the graph, start with every distance at 0, which is equivalent to adding a virtual source with a zero-weight edge to every vertex.

Time: the worst case versus typical inputs

Bellman-Ford is $O(VE)$ time and $O(V)$ space. With the early exit, it often finishes in far fewer than $V - 1$ rounds. On random sparse graphs with $E = 5V$:

n=1000 m=5000  dijkstra    0.9 ms   bellman-ford      2.1 ms
n=2000 m=10000  dijkstra    1.9 ms   bellman-ford      4.0 ms
n=4000 m=20000  dijkstra    4.1 ms   bellman-ford      7.8 ms

Random graphs have small diameter, so distances settle in a handful of rounds. The order in which edges are listed decides how far information travels per round. On a single path whose edges are listed backwards, each round advances the frontier by one edge:

path n=500: reverse order     5.8 ms   forward order  0.06 ms
path n=1000: reverse order    22.8 ms   forward order  0.12 ms
path n=2000: reverse order   103.3 ms   forward order  0.24 ms

Doubling $n$ multiplies the reverse-order time by about 4, which is the quadratic $V \cdot E$ behavior, while the forward order finishes in one productive round. The queue-based variant often called SPFA (only re-relax edges out of vertices whose distance just changed) is faster on typical inputs but has the same $O(VE)$ worst case.

Finding the cycle itself

To return the cycle rather than a flag, record parent pointers and run $V$ rounds. If the last round still relaxed some vertex $x$, then following parent pointers from $x$ for $V$ steps is guaranteed to land inside the cycle, from where one more walk traces it.

def find_negative_cycle(n, edges):
    """Return a list of nodes forming a negative cycle, or None."""
    dist = [0] * n               # virtual source to every node
    parent = [-1] * n
    x = -1
    for _ in range(n):
        x = -1
        for u, v, w in edges:
            if dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
                parent[v] = u
                x = v
    if x == -1:
        return None
    for _ in range(n):           # walk back n steps to land inside the cycle
        x = parent[x]
    cycle, y = [x], parent[x]
    while y != x:
        cycle.append(y)
        y = parent[y]
    return cycle[::-1]

Classic problem: currency arbitrage

Given exchange rates, is there a sequence of trades that ends with more money than it started with? A cycle with $r_1 r_2 \cdots r_k > 1$ is the same as $\sum -\log r_i < 0$, so weighting each edge by $-\log(\text{rate})$ turns arbitrage into negative-cycle detection. With made-up rates:

import math
names = ["USD", "EUR", "GBP", "JPY"]
rates = {
    ("USD", "EUR"): 0.92, ("EUR", "GBP"): 0.86, ("GBP", "USD"): 1.28,
    ("USD", "JPY"): 150.0, ("JPY", "EUR"): 0.0061,
}
idx = {x: i for i, x in enumerate(names)}
edges = [(idx[a], idx[b], -math.log(r)) for (a, b), r in rates.items()]
cycle = find_negative_cycle(4, edges)
print([names[i] for i in cycle])

['GBP', 'USD', 'JPY', 'EUR']
product of rates around the cycle: 1.0072

Classic problem: Cheapest Flights Within K Stops

Find the cheapest route from src to dst using at most k intermediate stops. Plain Dijkstra does not apply directly because the cheapest route might use too many stops. Bellman-Ford’s invariant is exactly what is needed: after round $i$, distances reflect paths of at most $i$ edges. Run $k + 1$ rounds, with one change: relax from a copy of the previous round’s distances, so that a single round cannot chain several edges together.

def cheapest_flights(n, flights, src, dst, k):
    dist = [INF] * n
    dist[src] = 0
    for _ in range(k + 1):              # at most k stops = k + 1 edges
        prev = dist[:]                  # relax from the previous round only
        for u, v, w in flights:
            if prev[u] + w < dist[v]:
                dist[v] = prev[u] + w
    return -1 if dist[dst] == INF else dist[dst]

With flights 0→1 100, 1→2 100, 2→0 100, 1→3 600, 2→3 200 and k = 1, the answer is 700 (0→1→3). The cheaper route 0→1→2→3 costs 400 but needs two stops. Without the copy, the function would return 400.

Floyd-Warshall: all pairs

The problem and the invariant

Floyd-Warshall computes the shortest distance between every pair of vertices. It was published by Robert Floyd in 1962, building on Stephen Warshall’s 1962 transitive-closure algorithm (Bernard Roy described the closure version in 1959).

It is a dynamic program over which vertices are allowed as intermediate stops. Let $D_k[i][j]$ be the shortest path from $i$ to $j$ whose intermediate vertices all come from $\{0, \ldots, k-1\}$. Allowing vertex $k$ as well, the best path either avoids $k$ or passes through it once:

$D_{k+1}[i][j] = \min\left(D_k[i][j],\ D_k[i][k] + D_k[k][j]\right)$

Row $k$ and column $k$ do not change during iteration $k$ (since $D[k][k] = 0$ without negative cycles), so the update can be done in place in a single matrix. The loop over $k$ must be the outermost. Putting it inside is a common bug that produces wrong answers without any error.

Implementation with path reconstruction

To reconstruct paths, keep nxt[i][j]: the first step on the best path from $i$ to $j$. When the path through $k$ wins, the first step from $i$ toward $j$ becomes the first step from $i$ toward $k$.

def floyd_warshall(n, edges):
    dist = [[INF] * n for _ in range(n)]
    nxt = [[None] * n for _ in range(n)]
    for i in range(n):
        dist[i][i] = 0
        nxt[i][i] = i
    for u, v, w in edges:
        if w < dist[u][v]:
            dist[u][v] = w
            nxt[u][v] = v
    for k in range(n):
        dk = dist[k]
        for i in range(n):
            di = dist[i]
            dik = di[k]
            if dik == INF:
                continue
            ni = nxt[i]
            for j in range(n):
                nd = dik + dk[j]
                if nd < di[j]:
                    di[j] = nd
                    ni[j] = ni[k]
    return dist, nxt


def fw_path(nxt, i, j):
    if nxt[i][j] is None:
        return None
    path = [i]
    while i != j:
        i = nxt[i][j]
        path.append(i)
    return path

Binding rows to local names (dk, di) and skipping rows where dist[i][k] is infinite are Python-specific speedups. They do not change the algorithm.

edges = [(0, 1, 3), (0, 2, 8), (1, 2, 2), (2, 3, 1), (3, 0, 4), (1, 3, 7)]
dist, nxt = floyd_warshall(4, edges)
for row in dist:
    print(row)
print(fw_path(nxt, 0, 3), fw_path(nxt, 3, 2))

[0, 3, 5, 6]
[7, 0, 2, 3]
[5, 8, 0, 1]
[4, 7, 9, 0]
[0, 1, 2, 3] [3, 0, 1, 2]

The direct edge 0→2 costs 8, but going through 1 costs 5. The direct edge 1→3 costs 7, but 1→2→3 costs 3.

Complexity and measurements

Time is $\Theta(V^3)$ regardless of the number of edges; space is $\Theta(V^2)$. On random graphs with $10V$ edges:

n=100      20.7 ms  
n=200     155.1 ms  7.5x
n=400    1416.4 ms  9.1x

Each doubling costs about 8 times as much, as expected for a cubic algorithm. In pure Python, a few hundred vertices is the practical limit for interactive use; vectorized code (NumPy updating one $k$ at a time) moves the limit to thousands.

When to use

Negative cycles: after the algorithm, $D[i][i] < 0$ for every vertex $i$ on a negative cycle. Any pair whose path can touch such a cycle has no shortest path. Numbers in the matrix can also grow exponentially negative in that case, which is harmless with Python ints but overflows fixed-width integers.

A classic use is Find the City With the Smallest Number of Neighbors at a Threshold Distance: run Floyd-Warshall once and count, for each city, how many entries in its row are within the threshold.

0-1 BFS

When every edge weighs 0 or 1, Dijkstra’s priority queue can be replaced by a deque. The key observation: the values of $d$ in Dijkstra’s queue are always within one of each other here. Popping a vertex at distance $d$ creates entries at $d$ (weight-0 edges) or $d + 1$ (weight-1 edges). So putting weight-0 neighbors on the front and weight-1 neighbors on the back keeps the deque sorted, with at most two distinct distance values in it at any time.

from collections import deque


def zero_one_bfs(graph, source):
    dist = {source: 0}
    dq = deque([source])
    while dq:
        u = dq.popleft()
        for v, w in graph.get(u, ()):
            nd = dist[u] + w
            if nd < dist.get(v, INF):
                dist[v] = nd
                if w == 0:
                    dq.appendleft(v)
                else:
                    dq.append(v)
    return dist

graph = {
    "a": [("b", 0), ("c", 1)],
    "b": [("c", 1), ("d", 1)],
    "c": [("d", 0), ("e", 1)],
    "d": [("e", 0)],
}
print(zero_one_bfs(graph, "a"))

{'a': 0, 'b': 0, 'c': 1, 'd': 1, 'e': 1}

Time is $O(V + E)$ and space $O(V)$. A vertex can be pushed twice: first at $d + 1$ through a weight-1 edge, then at $d$ through a weight-0 edge found later. It cannot improve a third time, because nothing in the deque is below $d$. Its second pop finds nothing to improve, so the extra work is bounded by a constant factor.

Where 0-1 weights come from: “free” moves versus “costly” ones. Examples: the minimum number of walls to remove to cross a grid (moving into an empty cell is free, into a wall costs 1), the minimum number of edge reversals to make a path exist in a directed graph (following an edge costs 0, traversing one backwards costs 1), and the minimum number of arrow changes to reach the bottom-right of a grid of arrows.

Classic problem: Minimum Obstacle Removal to Reach Corner

def min_obstacles(grid):
    """0 = empty, 1 = obstacle to remove. Moving into a cell costs grid value."""
    rows, cols = len(grid), len(grid[0])
    dist = [[INF] * cols for _ in range(rows)]
    dist[0][0] = 0
    dq = deque([(0, 0)])
    while dq:
        r, c = dq.popleft()
        for nr, nc in ((r + 1, c), (r - 1, c), (r, c + 1), (r, c - 1)):
            if 0 <= nr < rows and 0 <= nc < cols:
                w = grid[nr][nc]
                if dist[r][c] + w < dist[nr][nc]:
                    dist[nr][nc] = dist[r][c] + w
                    if w == 0:
                        dq.appendleft((nr, nc))
                    else:
                        dq.append((nr, nc))
    return dist[rows - 1][cols - 1]

For [[0, 1, 1], [1, 1, 0], [1, 1, 0]] the answer is 2. The generalization to weights $0, 1, \ldots, k$ is Dial’s algorithm, which uses $k + 1$ buckets in a circular array.

A* search

The problem and the idea

Dijkstra expands vertices in order of distance from the source, in every direction. When there is one specific goal and some estimate of how far each vertex is from it, the search can be steered. A* (Peter Hart, Nils Nilsson and Bertram Raphael, 1968) orders the frontier by

$f(v) = g(v) + h(v)$

where $g(v)$ is the best known cost from the start, as in Dijkstra, and $h(v)$ is a heuristic estimate of the remaining cost to the goal. Dijkstra is A* with $h = 0$.

Admissible and consistent heuristics

• Admissible: $h(v) \le \delta(v, \text{goal})$ for every $v$. The estimate never overshoots. This guarantees the first time the goal is popped, its cost is optimal, provided the search allows a vertex to be reopened when a cheaper path to it is found later.
• Consistent (also called monotone): $h(u) \le w(u, v) + h(v)$ for every edge, and $h(\text{goal}) = 0$. This is a triangle inequality. Consistency implies admissibility, and it gives the stronger guarantee that a vertex’s $g$ is final the first time it is popped, so a closed set works exactly as in Dijkstra.

The cleanest way to see why consistency works: A* with a consistent $h$ is Dijkstra on the reweighted graph $w'(u, v) = w(u, v) - h(u) + h(v)$. Consistency says $w' \ge 0$, so Dijkstra’s proof applies. Every path from start to goal changes by the same constant $h(\text{goal}) - h(\text{start})$, so the shortest path stays the same. (This is the same potential trick as Johnson’s algorithm, with $-h$ as the potential.)

An admissible but inconsistent heuristic combined with a closed set that never reopens vertices can return a suboptimal path. In practice, most heuristics built from a relaxed version of the problem (straight-line distance, ignoring walls) are consistent.

Euclidean distance is admissible on a 4-directional grid too, but weaker than Manhattan, so it expands more vertices. A heuristic that overestimates (for example $2 \times$ Manhattan, “weighted A*”) expands fewer vertices but gives up optimality; weighted A* with weight $\varepsilon$ returns a path at most $\varepsilon$ times the optimum.

Implementation on a grid

from itertools import count


def grid_neighbors(grid, r, c):
    rows, cols = len(grid), len(grid[0])
    for dr, dc in ((1, 0), (-1, 0), (0, 1), (0, -1)):
        nr, nc = r + dr, c + dc
        if 0 <= nr < rows and 0 <= nc < cols and grid[nr][nc] != "#":
            yield nr, nc


def manhattan(a, b):
    return abs(a[0] - b[0]) + abs(a[1] - b[1])


def astar_grid(grid, start, goal, h):
    """Unit-cost 4-directional grid. Returns (distance, nodes_expanded)."""
    g = {start: 0}
    tie = count()
    heap = [(h(start, goal), next(tie), start)]
    closed = set()
    expanded = 0
    while heap:
        f, _, u = heapq.heappop(heap)
        if u in closed:
            continue
        closed.add(u)
        expanded += 1
        if u == goal:
            return g[u], expanded
        for v in grid_neighbors(grid, *u):
            ng = g[u] + 1
            if ng < g.get(v, INF):
                g[v] = ng
                heapq.heappush(heap, (ng + h(v, goal), next(tie), v))
    return None, expanded

Passing h = lambda a, b: 0 turns this into Dijkstra, which makes a fair comparison: identical code, only the heuristic differs. A second variant breaks ties on equal $f$ by preferring the smaller $h$, that is, the vertex closer to the goal. The heap key becomes (ng + hv, hv, next(tie), v); nothing else changes.

Measured node expansions

On 300 × 300 grids, all three searches returned the same distance in every case (checked). “Expanded” counts vertices popped and closed.

grid                       dist  Dijkstra  A* FIFO ties  A* ties by h
open, corner to corner      598    90,000        90,000           599
                             ms     119.5         126.9           0.9
open, 200 apart in a row    200    52,800           201           201
                             ms      70.8           0.3           0.4
20% obstacles               598    71,662        40,099         6,043
                             ms      89.8          61.9          10.7
30% obstacles               602    61,938        15,929        12,197
                             ms      77.9          21.5          21.2
wall, gap at far end        879    89,710        87,010        43,941
                             ms     114.1         118.8          56.2

(The “wall” grid has a vertical wall down the middle with a gap only at the far end, start and goal on either side of it at the bottom; the true distance of 879 is far above the Manhattan estimate of 299.)

What these numbers show:

1. When the heuristic is exact, tie-breaking decides everything. On an open grid from corner to corner, every cell lies on some shortest path, so every cell has $f = 598$. With first-in-first-out ties, A* expands all 90,000 cells, exactly like Dijkstra. Preferring the smaller $h$ among equal $f$ makes it walk straight to the goal: 599 expansions, about 130 times faster.
2. When start and goal are in a straight line, cells off the line have larger $f$, so A* expands only the 201 cells on the path, against Dijkstra’s 52,800 cells within distance 200.
3. Obstacles reduce the advantage. The heuristic ignores walls, so the more the true path deviates from the straight line, the more the estimate undershoots and the more A* behaves like Dijkstra.
4. A heuristic that is badly wrong buys little. With the wall forcing a long detour, A* with FIFO ties expanded 97% as many cells as Dijkstra, and paid the extra cost of computing $h$.

A* is the standard algorithm for pathfinding in games and robotics, and for puzzle solving (the 15-puzzle with the Manhattan-distance-sum heuristic). Its memory use, which is all of the frontier plus the closed set, is the usual limit; IDA* (iterative deepening A*) trades time for memory when that matters.

Minimum spanning trees

The problem

Given a connected, undirected, weighted graph, a spanning tree is a set of $V - 1$ edges that connects all vertices without a cycle. A minimum spanning tree (MST) is one with the smallest total weight. Applications: designing the cheapest network that connects all sites (cable, pipes, roads), clustering (remove the $k - 1$ heaviest MST edges to get $k$ clusters, which is single-linkage clustering), and as a building block for approximation algorithms such as the 2-approximation for metric TSP.

If the graph is not connected, the same algorithms produce a minimum spanning forest, one tree per component.

An MST is not a shortest-path tree. The MST minimizes the total edge weight; a shortest-path tree minimizes each vertex’s distance from one root. On a cycle of four vertices with weights 1, 1, 1, 2.5, the MST drops the 2.5 edge, and the vertex at the far end of that edge is 3 away from its neighbor in the tree instead of 2.5.

The cut property

Both classic algorithms rest on one fact. A cut splits the vertices into two non-empty sets; an edge crosses the cut if its endpoints are on different sides.

Cut property. For any cut, a minimum-weight edge crossing it belongs to some MST. If it is the unique minimum, it belongs to every MST.

Proof by exchange: take an MST $T$ that does not contain the light edge $e = (u, v)$. Adding $e$ to $T$ creates a cycle. That cycle crosses the cut at $e$, so it must cross it again at some other edge $e'$. Remove $e'$: the result is still a spanning tree, and its weight changed by $w(e) - w(e') \le 0$. So a tree containing $e$ is also minimal.

The companion cycle property: the unique heaviest edge on any cycle is in no MST. And if all edge weights are distinct, the MST is unique.

Prim’s algorithm applies the cut property to the cut between “the tree so far” and “everything else”. Kruskal’s applies it to the cut between “one component of the forest so far” and “everything else”.

Kruskal’s algorithm

Kruskal’s algorithm (Joseph Kruskal, 1956): sort all edges by weight, and add each edge unless it would join two vertices already in the same component. Union-find, from chapter 12, answers “same component?” in near-constant time.

def kruskal(n, edges):
    parent = list(range(n))
    size = [1] * n

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    total, chosen = 0, []
    for w, u, v in sorted((w, u, v) for u, v, w in edges):
        ru, rv = find(u), find(v)
        if ru == rv:
            continue
        if size[ru] < size[rv]:
            ru, rv = rv, ru
        parent[rv] = ru
        size[ru] += size[rv]
        total += w
        chosen.append((u, v, w))
        if len(chosen) == n - 1:
            break
    return total, chosen

Why each added edge is safe: when Kruskal considers edge $(u, v)$ with $u$ and $v$ in different components, take the cut between $u$‘s component and the rest. Every lighter edge has already been considered, and none of them crosses this cut (each was either added, merging its endpoints into one component, or skipped because its endpoints were already together). So $(u, v)$ is a lightest crossing edge.

Time: $O(E \log E)$ for the sort, plus $O(E \, \alpha(V))$ for union-find. Space $O(V + E)$.

Prim’s algorithm

Prim’s algorithm (Robert Prim, 1957; Vojtěch Jarník described it in 1930, and Dijkstra again in 1959) grows a single tree from a start vertex, always adding the lightest edge that leaves the tree. It looks like Dijkstra with one change: the heap key is the edge weight, not the distance from the start.

def prim(n, adj):
    """adj: list of lists of (v, w), undirected. Returns total weight of MST/forest."""
    seen = [False] * n
    total = 0
    for s in range(n):
        if seen[s]:
            continue
        heap = [(0, s)]
        while heap:
            w, u = heapq.heappop(heap)
            if seen[u]:
                continue
            seen[u] = True
            total += w
            for v, wv in adj[u]:
                if not seen[v]:
                    heapq.heappush(heap, (wv, v))
    return total

The outer loop over s makes it handle disconnected graphs by starting a new tree in each component. With a binary heap and lazy deletion, time is $O(E \log E)$; with an array and a linear scan for the minimum, it is $O(V^2)$, which is optimal for dense graphs.

Walk-through and verification

         4
    0 ------- 1
     \       / \
    3 \   1 /   \ 2
       \   /     \
         2 ------ 3 ---2--- 4
             4     \        |
                  7 \       | 6
                     \      |
                      +---- 5

Edges: 0–1 4, 0–2 3, 1–2 1, 1–3 2, 2–3 4, 3–4 2, 4–5 6, 3–5 7.

Kruskal, in sorted order: 1–2 (1) add; 1–3 (2) add; 3–4 (2) add; 0–2 (3) add; 0–1 (4) skip, 0 and 1 are connected through 2; 2–3 (4) skip; 4–5 (6) add. Five edges for six vertices, and done.

(14, [(1, 2, 1), (1, 3, 2), (3, 4, 2), (0, 2, 3), (4, 5, 6)])
14

The second line is Prim’s total on the same graph. Prim and Kruskal may choose different edges when weights tie, but the total weight is always the same. Checked on 1,000 random graphs with weights from 1 to 10, so with many ties, and graphs that were sometimes disconnected:

Prim and Kruskal agreed on 1000/1000 random graphs (weights 1..10, many ties)

Both were also checked against an exhaustive search over all edge subsets on small graphs.

Which to use

V=10000 E=50000: kruskal   23.6 ms  prim   31.5 ms
V=10000 E=500000: kruskal  368.4 ms  prim  994.7 ms

In Python, Kruskal tends to win on edge lists because sorted runs in C, while Prim pushes and pops through heapq once per edge. In compiled languages the picture is closer.

Classic problem: Min Cost to Connect All Points

Connect points in the plane with the minimum total Manhattan distance. The graph is complete, with $V(V-1)/2$ implicit edges, so array-based Prim at $O(V^2)$ beats anything that builds and sorts the edge list.

def min_cost_connect_points(points):
    """Dense Prim, O(V^2), Manhattan distances."""
    n = len(points)
    best = [INF] * n
    in_tree = [False] * n
    best[0] = 0
    total = 0
    for _ in range(n):
        u = min((i for i in range(n) if not in_tree[i]), key=best.__getitem__)
        in_tree[u] = True
        total += best[u]
        ux, uy = points[u]
        for v in range(n):
            if not in_tree[v]:
                d = abs(ux - points[v][0]) + abs(uy - points[v][1])
                if d < best[v]:
                    best[v] = d
    return total

For [(0, 0), (2, 2), (3, 10), (5, 2), (7, 0)] it returns 20, matching Kruskal on the explicit edge list.

Strongly connected components

The problem

In a directed graph, two vertices are strongly connected if each can reach the other. This is an equivalence relation, so it partitions the vertices into strongly connected components (SCCs). Collapsing each SCC to a single vertex gives the condensation, which is always a DAG. Many problems on general directed graphs reduce to a problem on that DAG: 2-SAT (a formula is unsatisfiable exactly when some variable and its negation share an SCC), finding the minimum number of edges to add to make a graph strongly connected, deadlock analysis, and finding mutually dependent modules in a codebase.

   0 --> 1 --> 3 --> 4 --> 5
   ^     |           ^     |
   |     v           |     v
   +---- 2           +---- 6       4 --> 7

   SCCs: {0,1,2}  {3}  {4,5,6}  {7}
   Condensation: {0,1,2} -> {3} -> {4,5,6} -> {7}

Kosaraju’s algorithm

Kosaraju’s algorithm (attributed to S. Rao Kosaraju, 1978, unpublished; published by Micha Sharir in 1981) uses two DFS passes.

1. Run DFS on the graph and record vertices in order of finish time.
2. Reverse every edge. Run DFS on the reversed graph, taking start vertices in order of decreasing finish time. Each DFS tree in this pass is one SCC.

Why it works: the vertex that finishes last lies in a source component of the condensation (no edges come into that component from other components). Reversing the edges turns it into a sink, so a DFS from it in the reversed graph cannot escape its own component. Once that component is removed, the next vertex in finish order lies in a source of the remaining condensation, and so on.

def kosaraju(n, adj):
    order, seen = [], [False] * n
    for s in range(n):
        if seen[s]:
            continue
        seen[s] = True
        stack = [(s, iter(adj[s]))]
        while stack:
            u, it = stack[-1]
            for v in it:
                if not seen[v]:
                    seen[v] = True
                    stack.append((v, iter(adj[v])))
                    break
            else:
                stack.pop()
                order.append(u)          # post-order (finish time)
    radj = [[] for _ in range(n)]
    for u in range(n):
        for v in adj[u]:
            radj[v].append(u)
    comp = [-1] * n
    c = 0
    for s in reversed(order):
        if comp[s] != -1:
            continue
        comp[s] = c
        stack = [s]
        while stack:
            u = stack.pop()
            for v in radj[u]:
                if comp[v] == -1:
                    comp[v] = c
                    stack.append(v)
        c += 1
    return comp

The first pass must produce a true post-order, so it is written with an iterator per stack frame: a vertex is appended to order only after all of its neighbors are finished. The usual “pop, then push all neighbors” iterative DFS does not give correct finish times. The second pass only needs reachability, so any traversal works there.

Tarjan’s algorithm

Tarjan’s algorithm (Robert Tarjan, 1972) finds SCCs in a single DFS. Each vertex gets a discovery index and a low-link value: the smallest index reachable from its DFS subtree using tree edges plus at most one edge back to a vertex still on the stack. Vertices are pushed onto an explicit stack when discovered. A vertex whose low-link equals its own index is the root of an SCC: nothing in its subtree reaches above it, so the vertices on the stack above it, plus itself, form one component and are popped together.

def tarjan_scc(n, adj):
    index = [-1] * n
    low = [0] * n
    on_stack = [False] * n
    stack, comps = [], []
    counter = 0

    def visit(u):
        nonlocal counter
        index[u] = low[u] = counter
        counter += 1
        stack.append(u)
        on_stack[u] = True
        for v in adj[u]:
            if index[v] == -1:
                visit(v)
                low[u] = min(low[u], low[v])
            elif on_stack[v]:
                low[u] = min(low[u], index[v])
        if low[u] == index[u]:           # u is the root of an SCC
            comp = []
            while True:
                w = stack.pop()
                on_stack[w] = False
                comp.append(w)
                if w == u:
                    break
            comps.append(comp)

    for s in range(n):
        if index[s] == -1:
            visit(s)
    return comps

The on_stack check is essential. An edge to a vertex that is already finished and assigned to another SCC is a cross edge into a different component and must not lower the low-link.

On the example graph above (adjacency [[1], [2, 3], [0], [4], [5, 7], [6], [4], []]):

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

The first line is Kosaraju’s component label per vertex; its components come out in topological order of the condensation. The second is Tarjan’s list of components, which come out in reverse topological order (sinks first). Both properties were checked on random graphs, and they are useful: a DP over the condensation can consume Tarjan’s output directly.

Comparison

Recursion depth in Python is limited to 1,000 frames by default, so a recursive Tarjan fails on a path of a few thousand vertices. Either raise the limit with sys.setrecursionlimit (and possibly the thread stack size) or convert to an explicit stack, which takes some care because the low-link update after returning from a child must be preserved. An iterative version was tested alongside the recursive one shown.

Bridges and articulation points

The problem

In an undirected connected graph, a bridge is an edge whose removal disconnects the graph. An articulation point (cut vertex) is a vertex whose removal disconnects it. They identify single points of failure: a network link whose loss partitions a data center, a router that splits the network, a road whose closure isolates a town.

The brute-force approach removes each edge or vertex and re-runs a connectivity check, costing $O(E(V + E))$ or $O(V(V + E))$. Tarjan’s low-link technique finds all of them in one DFS, $O(V + E)$.

Low-link values

Run DFS and record disc[u], the time $u$ was discovered. In an undirected graph, DFS produces only tree edges and back edges (edges to an ancestor); there are no cross edges. Define low[u] as the smallest discovery time reachable from $u$‘s subtree using tree edges downward and then one back edge.

• Tree edge $(u, v)$ is a bridge if and only if low[v] > disc[u]. Nothing in $v$‘s subtree has a back edge to $u$ or above, so the only connection is $(u, v)$ itself.
• A non-root vertex $u$ is an articulation point if it has a child $v$ with low[v] >= disc[u]. $v$’s subtree can reach $u$ but not above it, so removing $u$ cuts it off. (Note $\ge$ for vertices, $>$ for edges.)
• The DFS root is an articulation point if and only if it has two or more DFS children. The general rule does not apply to it, because nothing is above it.

def bridges_and_cut_vertices(n, adj):
    """Undirected simple graph; adj[u] = list of neighbors. Recursive."""
    disc = [-1] * n
    low = [0] * n
    bridges, cuts = [], set()
    timer = 0

    def dfs(u, parent):
        nonlocal timer
        disc[u] = low[u] = timer
        timer += 1
        children = 0
        for v in adj[u]:
            if v == parent:
                continue
            if disc[v] != -1:                 # back edge
                low[u] = min(low[u], disc[v])
            else:                              # tree edge
                children += 1
                dfs(v, u)
                low[u] = min(low[u], low[v])
                if low[v] > disc[u]:
                    bridges.append((u, v))
                if parent != -1 and low[v] >= disc[u]:
                    cuts.add(u)
        if parent == -1 and children > 1:
            cuts.add(u)

    for s in range(n):
        if disc[s] == -1:
            dfs(s, -1)
    return bridges, cuts

   0           4
   | \       / |
   |  2 --- 3  |
   | /       \ |
   1           5 --- 6

Two triangles joined by the edge 2–3, plus a pendant edge 5–6:

bridges: [(5, 6), (2, 3)]
articulation points: [2, 3, 5]

Both outputs were checked against the brute-force “remove and count components” definition on 500 random graphs.

Pitfall: parallel edges. Skipping the edge back to the parent vertex is correct only for simple graphs. With two parallel edges between $u$ and $v$, neither is a bridge, but the code above would report one. For multigraphs, give each edge an id and skip only the specific edge used to arrive.

Classic problem: Critical Connections in a Network asks exactly for the bridges; for the input [[0, 1], [1, 2], [2, 0], [1, 3]] the only one is [1, 3].

A related concept is biconnected components: maximal subgraphs with no articulation point. They are found with the same DFS by keeping a stack of edges and popping a component whenever low[v] >= disc[u].

Eulerian paths and circuits

The problem

An Eulerian path uses every edge exactly once. An Eulerian circuit is one that starts and ends at the same vertex. Euler’s 1736 analysis of the Königsberg bridges, which showed no such walk existed there, is usually counted as the beginning of graph theory.

Unlike the Hamiltonian problem (every vertex once), this one is solvable in linear time: existence depends only on degrees and connectivity.

The intuition for the degree condition: every time a walk passes through a vertex, it uses one edge in and one edge out. Only the start and end can be unbalanced.

Hierholzer’s algorithm

Hierholzer’s algorithm (Carl Hierholzer, published 1873) builds the path in $O(E)$. Walk from the start, using up edges, until stuck. Being stuck can only happen at the end vertex. Then backtrack along the walk; at any vertex that still has unused edges, start a new sub-walk from there, which must return to that vertex, and splice it in. With a stack, the splicing is automatic: a vertex is appended to the output when it has no unused edges left, and the output, reversed, is the path.

from collections import defaultdict


def eulerian_path_directed(edges):
    """edges: list of (u, v). Returns node sequence or None."""
    adj = defaultdict(list)
    indeg, outdeg = defaultdict(int), defaultdict(int)
    for u, v in edges:
        adj[u].append(v)
        outdeg[u] += 1
        indeg[v] += 1
    nodes = set(indeg) | set(outdeg)
    start = next(iter(nodes))
    plus = minus = 0
    for x in nodes:
        d = outdeg[x] - indeg[x]
        if d == 1:
            plus += 1
            start = x
        elif d == -1:
            minus += 1
        elif d != 0:
            return None
    if not ((plus == minus == 0) or (plus == minus == 1)):
        return None
    if plus == 0:
        start = next(x for x in nodes if outdeg[x] > 0)
    stack, path = [start], []
    while stack:
        u = stack[-1]
        if adj[u]:
            stack.append(adj[u].pop())
        else:
            path.append(stack.pop())
    path.reverse()
    return path if len(path) == len(edges) + 1 else None   # connectivity check

The final length check replaces a separate connectivity test: if the edges are split across disconnected parts, the walk from start cannot use all of them.

Walk-through on edges 0→1, 1→2, 2→0, 0→3, 3→4. Vertex 0 has out − in = 1, so it is the start. adj[0] is [1, 3] and pop() takes 3 first:

stack [0]          take 0->3        stack [0, 3]
stack [0, 3]       take 3->4        stack [0, 3, 4]
4 has no edges     output 4         stack [0, 3]
3 has no edges     output 3         stack [0]
0 has 0->1 left    take it          stack [0, 1] -> [0, 1, 2] -> [0, 1, 2, 0]
then everything unwinds:            output 0, 2, 1, 0
output = [4, 3, 0, 2, 1, 0], reversed = [0, 1, 2, 0, 3, 4]

The first walk went straight to the dead end 4; the cycle 0→1→2→0 was spliced in before it automatically.

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

The second call, on edges 0→1 and 0→2, returns None: vertex 0 has out − in = 2.

For undirected graphs the same algorithm works, but each edge appears in two adjacency lists and must be marked used by id so it is not traversed twice.

Classic problem: Reconstruct Itinerary

Given airline tickets as [from, to] pairs, reconstruct the trip that starts at "JFK" and uses every ticket exactly once; if several exist, return the lexicographically smallest. This is an Eulerian path, with a greedy twist: always take the smallest available destination. Hierholzer’s splicing makes the greedy choice safe, because a smallest-first walk that hits a dead end too early just becomes the tail of the final route.

def find_itinerary(tickets):
    adj = defaultdict(list)
    for a, b in sorted(tickets, reverse=True):
        adj[a].append(b)                 # smallest destination ends up last
    stack, route = ["JFK"], []
    while stack:
        while adj[stack[-1]]:
            stack.append(adj[stack[-1]].pop())
        route.append(stack.pop())
    return route[::-1]

['JFK', 'MUC', 'LHR', 'SFO', 'SJC']
['JFK', 'ATL', 'JFK', 'SFO', 'ATL', 'SFO']

This was checked against brute force over all ticket orderings. Eulerian paths also appear in DNA sequence assembly (de Bruijn graphs, where reads become edges) and in the Cracking the Safe problem (the shortest string containing every $n$-digit code is an Eulerian circuit of a de Bruijn graph).

Hamiltonian paths and circuits

A Hamiltonian path visits every vertex exactly once; a Hamiltonian circuit (cycle) returns to the start. The definition differs from the Eulerian one by a single word, but the difficulty is completely different: deciding whether a Hamiltonian path or cycle exists is NP-complete, both for directed and undirected graphs (the directed and undirected Hamiltonian cycle problems are on Richard Karp’s 1972 list of 21 NP-complete problems). No polynomial-time algorithm is known, and none exists unless P = NP. No simple degree condition settles it either; there are only sufficient conditions, such as Dirac’s theorem: if every vertex in a graph with $n \ge 3$ vertices has degree at least $n/2$, a Hamiltonian cycle exists.

Approaches

• Backtracking extends a path one vertex at a time and undoes the choice when stuck. Worst case $O(n!)$, with pruning helping a lot on some inputs. Chapter 15 covers the technique.
• Bitmask DP (Bellman, and Held and Karp, 1962) solves it in $O(2^n n^2)$ time and $O(2^n n)$ space. The state is “the set of vertices visited, and which one the path currently ends at”. Since paths with the same visited set and endpoint are interchangeable for the future, only existence needs to be stored. That collapses $n!$ orderings into $2^n \cdot n$ states, which is practical up to about $n = 20$ in compiled languages.

The version below stores, for each subset, a bitset of possible end vertices, so the state is one int per subset:

def has_hamiltonian_path(n, adj_matrix):
    full = 1 << n
    dp = [0] * full          # dp[mask] = bitset of end nodes reachable
    for v in range(n):
        dp[1 << v] = 1 << v
    nbr = [sum(1 << j for j in range(n) if adj_matrix[i][j]) for i in range(n)]
    for mask in range(full):
        ends = dp[mask]
        while ends:
            low = ends & -ends
            v = low.bit_length() - 1
            ends ^= low
            nxt = nbr[v] & ~mask
            while nxt:
                b = nxt & -nxt
                nxt ^= b
                dp[mask | b] |= b
    return dp[full - 1] != 0

Processing masks in increasing numeric order is valid because adding a bit always produces a larger number, so every state is complete before it is read. Checked against all permutations on 300 random graphs with up to 7 vertices. The bit tricks used here (x & -x isolates the lowest set bit) are explained in chapter 17.

Classic problem: Shortest Path Visiting All Nodes asks for the shortest walk (vertices may repeat) that visits every vertex of an unweighted graph. It is BFS over states (current vertex, visited mask), $O(2^n \cdot n)$ states, a close relative of the DP above.

The Traveling Salesman Problem

The problem

Given $n$ cities and the distance between each pair, find the shortest tour that visits every city exactly once and returns to the start. It is the weighted version of the Hamiltonian cycle problem, and its decision version is NP-complete. Trying every tour costs $(n-1)!/2$ for symmetric distances; at $n = 20$ that is about $6 \times 10^{16}$ tours.

Held-Karp dynamic programming

Fix city 0 as the start. Let $C(S, j)$ be the length of the shortest path that starts at 0, visits exactly the cities in $S$ (which contains 0 and $j$), and ends at $j$. The base case is $C(\{0\}, 0) = 0$, with $C(S, 0) = \infty$ for any larger $S$, since the path may not return to 0 early. Then for $j \ne 0$

$C(S, j) = \min_{i \in S \setminus \{j\}} \Big( C(S \setminus \{j\}, i) + d(i, j) \Big)$

and the optimal tour is $\min_j \big( C(\text{all}, j) + d(j, 0) \big)$. There are $2^{n-1} \cdot n$ states and each takes $O(n)$ work: $O(2^n n^2)$ time and $O(2^n n)$ space. This is the algorithm of Michael Held and Richard Karp (1962), also found independently by Richard Bellman the same year. It remains the best known worst-case exact bound for general TSP.

The implementation below is written “push” style: each state extends itself to every unvisited city.

def held_karp(dist):
    n = len(dist)
    full = 1 << n
    dp = [[INF] * n for _ in range(full)]
    dp[1][0] = 0
    for mask in range(1, full, 2):        # every set contains city 0
        row = dp[mask]
        for last in range(n):
            cost = row[last]
            if cost == INF:
                continue
            dl = dist[last]
            for nxt in range(n):
                if mask & (1 << nxt):
                    continue
                nm = mask | (1 << nxt)
                c = cost + dl[nxt]
                if c < dp[nm][nxt]:
                    dp[nm][nxt] = c
    last_row = dp[full - 1]
    return min(last_row[j] + dist[j][0] for j in range(1, n)) if n > 1 else 0

Stepping mask by 2 skips all subsets without city 0, halving the work. Checked against all $(n-1)!$ tours on 200 random asymmetric instances with up to 8 cities.

Measured growth

Random asymmetric instances, one per size:

n=10        1.5 ms  
n=11        3.6 ms  2.34x
n=12        8.5 ms  2.38x
n=13       19.5 ms  2.29x
n=14       51.5 ms  2.64x
n=15      119.1 ms  2.31x
n=16      260.9 ms  2.19x
n=17      632.1 ms  2.42x
n=18     1388.8 ms  2.20x
n=19     3279.7 ms  2.36x

Each additional city multiplies the time by about 2.3. The formula predicts $2 \cdot (n+1)^2 / n^2$, which is about 2.2 at these sizes. Memory grows the same way: at $n = 19$ the table has $2^{19} \times 19 \approx 10^7$ entries. Every three or so additional cities cost a factor of ten, so pure Python reaches its practical limit in the low twenties, and even compiled code stops in the high twenties. Beyond that, exact solvers use branch and bound with linear-programming relaxations (the Concorde solver has solved instances with tens of thousands of cities), and everyone else uses heuristics.

Nearest-neighbor heuristic and its quality

The nearest-neighbor heuristic starts at a city and repeatedly moves to the closest unvisited one, then returns home. It runs in $O(n^2)$.

def nearest_neighbor_tour(dist, start=0):
    n = len(dist)
    unvisited = set(range(n)) - {start}
    tour, cur, total = [start], start, 0
    while unvisited:
        nxt = min(unvisited, key=lambda j: dist[cur][j])
        total += dist[cur][nxt]
        unvisited.remove(nxt)
        tour.append(nxt)
        cur = nxt
    total += dist[cur][start]
    return total, tour

How far from optimal is it? On 100 random instances per size, with cities placed uniformly in the unit square and Euclidean distances, comparing with Held-Karp’s exact optimum (ratio 1.000 means optimal):

n= 8  mean 1.082  median 1.072  worst 1.354  optimal in 20/100
n=10  mean 1.109  median 1.090  worst 1.415  optimal in 10/100
n=12  mean 1.110  median 1.091  worst 1.386  optimal in 3/100

Running it from every possible start city and keeping the best tour costs $n$ times as much and helps substantially:

n= 8  mean 1.022  worst 1.129  optimal in 52/100
n=10  mean 1.027  worst 1.173  optimal in 38/100
n=12  mean 1.036  worst 1.242  optimal in 24/100

On these random instances nearest neighbor is about 10% above optimal on average, with occasional tours 40% too long, and it finds the optimum less often as $n$ grows. These figures are specific to small random Euclidean instances. Nearest neighbor has no constant-factor guarantee: even for metric distances, the best known bound on its ratio is $O(\log n)$, and there are constructed inputs on which the ratio grows without limit as $n$ increases.

Better options, in increasing order of effort:

The metric guarantees require the triangle inequality. For general (non-metric) distances, no constant-factor approximation exists unless P = NP.

Graph coloring

The problem

A proper coloring assigns each vertex a color so that adjacent vertices differ. The chromatic number $\chi(G)$ is the minimum number of colors needed. Coloring models conflicts: exam timetabling (vertices are exams, edges join exams sharing a student, colors are time slots), register allocation in compilers (vertices are variables, edges join variables live at the same time, colors are registers), frequency assignment for radio transmitters, and Sudoku (a 9-coloring of a fixed graph with some colors given).

• 2-coloring is the bipartite check from chapter 12, solvable in $O(V + E)$ by BFS.
• Deciding whether a graph is $k$-colorable is NP-complete for every fixed $k \ge 3$.
• Every planar graph is 4-colorable (the Four Color Theorem, proved by Appel and Haken in 1976 with computer assistance).

Backtracking for k-coloring

Assign colors to vertices one at a time, trying each color that no already-colored neighbor uses, and backtrack when a vertex has no legal color. Coloring high-degree vertices first fails earlier on bad branches.

def k_color(n, adj, k):
    color = [-1] * n
    order = sorted(range(n), key=lambda v: -len(adj[v]))

    def place(i):
        if i == n:
            return True
        v = order[i]
        used = {color[u] for u in adj[v]}
        for c in range(k):
            if c not in used:
                color[v] = c
                if place(i + 1):
                    return True
        color[v] = -1
        return False

    return color if place(0) else None

The worst case is $O(k^n)$. To find $\chi(G)$, try $k = 1, 2, 3, \ldots$ until it succeeds. On the Petersen graph (10 vertices, 15 edges, contains 5-cycles, so not bipartite):

Petersen 2-colorable: False
Petersen 3-coloring: [0, 1, 0, 1, 2, 1, 0, 2, 2, 1]

A cheap optimization is symmetry breaking: colors are interchangeable, so the first vertex can be fixed to color 0, and in general a vertex never needs a color more than one above the largest used so far. This cuts the search by up to a factor of $k!$.

Greedy coloring and its dependence on order

Greedy coloring visits vertices in some order and gives each the smallest color not used by an already-colored neighbor. It runs in $O(V + E)$ and never uses more than $\Delta + 1$ colors, where $\Delta$ is the maximum degree (each vertex has at most $\Delta$ neighbors to avoid).

def greedy_color(adj, order):
    color = {}
    for v in order:
        used = {color[u] for u in adj[v] if u in color}
        c = 0
        while c in used:
            c += 1
        color[v] = c
    return color

The result depends heavily on the order. For every graph, some order makes greedy optimal (order the vertices by color class in an optimal coloring), but a bad order can be very bad. The crown graph shows how bad: $2k$ vertices $u_0, \ldots, u_{k-1}$ and $v_0, \ldots, v_{k-1}$, with $u_i$ adjacent to $v_j$ whenever $i \ne j$. It is bipartite, so 2 colors suffice.

4 interleaved: 4  side by side: 2
8 interleaved: 8  side by side: 2
16 interleaved: 16  side by side: 2

In the interleaved order $u_0, v_0, u_1, v_1, \ldots$, the pair $u_i, v_i$ is not adjacent but each is adjacent to all earlier vertices of the other side, which hold colors $0..i-1$. So greedy uses $k$ colors, $k/2$ times the optimum, and the gap grows without bound.

Welsh-Powell

The Welsh-Powell algorithm (Dominic Welsh and Martin Powell, 1967) is greedy coloring with a specific order: decreasing degree. Their original description colors one class at a time (go down the sorted list giving color 1 to every vertex not adjacent to a color-1 vertex, then repeat with color 2), which produces exactly the same coloring as first-fit greedy in that order. Coloring the most constrained vertices first tends to help; it also gives the bound $\chi(G) \le \max_i \min(d_i + 1, i)$ where $d_1 \ge d_2 \ge \cdots$ are the sorted degrees.

def welsh_powell(n, adj):
    order = sorted(range(n), key=lambda v: -len(adj[v]))
    return greedy_color(adj, order)

How do these compare with the true chromatic number? On 50 random graphs with 30 vertices and edge probability 0.3, exact $\chi$ was computed by backtracking, and the table counts how many colors each method used above $\chi$ (the dict maps “extra colors” to “number of graphs”):

G(30, 0.3), 50 graphs; colors used minus chromatic number:
greedy, index order:    {0: 2, 1: 23, 2: 23, 3: 2}
Welsh-Powell:           {0: 12, 1: 33, 2: 5}
greedy, worst of 200 random orders: {2: 5, 3: 42, 4: 3}
exact backtracking: median 3.0 ms, max 263.1 ms per graph

Welsh-Powell was optimal on 12 of 50 graphs and at most 2 colors over; greedy in arbitrary index order was optimal on only 2. The worst of 200 random orders was always at least 2 colors over. Exact backtracking was fast on most of these small graphs, but the maximum was nearly 100 times the median, which is typical of exponential search: the cost depends on how quickly the instance prunes.

DSatur (Daniel Brélaz, 1979) is a stronger dynamic ordering: at each step color the vertex whose neighbors already use the most distinct colors (its “saturation”), breaking ties by degree. It is exact on bipartite graphs and usually beats Welsh-Powell.

A pointer: maximum flow

One more family deserves a mention. A flow network has capacities on directed edges, and the maximum flow problem asks how much can be sent from a source to a sink. The Ford-Fulkerson method (1956) repeatedly finds a path with spare capacity in the residual graph (which includes reverse edges for undoing flow) and pushes flow along it. Choosing that augmenting path by BFS, the Edmonds-Karp algorithm (1972), bounds the work at $O(VE^2)$; Dinic’s algorithm improves this to $O(V^2 E)$. The max-flow min-cut theorem says the maximum flow equals the minimum total capacity of edges whose removal separates source from sink. Bipartite matching, edge-disjoint paths and many assignment problems reduce to max flow.

Choosing a shortest-path algorithm

The DAG row deserves a note: when the graph has no cycles, processing vertices in topological order and relaxing each one’s outgoing edges gives exact shortest (or longest) paths in linear time, even with negative weights. Critical-path scheduling in project management is this algorithm with longest paths.

Summary