Python Algorithms · Volume V

Deep Dives,
Proofs & Patterns

A comprehensive reference covering divide & conquer, greedy methods, backtracking, number theory, randomized algorithms, bit manipulation, computational geometry, and advanced data structures — with mathematical proofs, complexity analysis, and interactive demonstrations.

8 Major Topics Theorem Proofs 6 Interactive Demos Master Theorem 40+ Algorithms Complexity Analysis Real Python Code
§01

Divide & Conquer

T(n) = aT(n/b) + f(n)

The Master Theorem

The Master Theorem provides a cookbook solution for recurrences of the form T(n) = aT(n/b) + f(n) — where we split a problem of size n into a sub-problems each of size n/b, then combine in f(n) time.

⚡ Master Theorem (Cook-Levin-Akra-Bazzi)

Let T(n) = aT(n/b) + f(n) where a ≥ 1, b > 1. Define the watershed exponent c* = logb(a).

Case 1: f(n) = O(n^(c*−ε)) → T(n) = Θ(n^c*) [subproblems dominate] Case 2: f(n) = Θ(n^c*) → T(n) = Θ(n^c* · log n) [balanced] Case 3: f(n) = Ω(n^(c*+ε)) and regularity holds → T(n) = Θ(f(n)) [combine step dominates]

Worked examples:

Merge Sort: T(n) = 2T(n/2) + Θ(n) → c* = log₂(2) = 1, f = Θ(n¹) → Case 2 → Θ(n log n) Binary Search: T(n) = T(n/2) + Θ(1) → c* = log₂(1) = 0, f = Θ(n⁰) → Case 2 → Θ(log n) Karatsuba: T(n) = 3T(n/2) + Θ(n) → c* = log₂(3) ≈ 1.585 → Case 1 → Θ(n^1.585) Strassen: T(n) = 7T(n/2) + Θ(n²) → c* = log₂(7) ≈ 2.807 → Case 1 → Θ(n^2.807)
Proof sketch — recursion tree method

At depth k of the recursion tree there are a^k nodes, each handling a subproblem of size n/b^k. The work at depth k is:

a^k · f(n / b^k)

The tree has log_b(n) levels. Total work is the geometric series:

T(n) = Σ_{k=0}^{log_b n} a^k · f(n / b^k)

In Case 2, f(n) = Θ(n^c*) so each level contributes equal work n^c*, and with log_b(n) levels: T(n) = Θ(n^c* · log n). In Case 1, the root dominates each level and the geometric series sums to Θ(n^c*). In Case 3, the leaves dominate and f(n) wins. ∎

Binary Search — O(log n) Guaranteed

Python · Binary Search Suite
O(log n) O(1) space
def binary_search(arr, target):
    """Standard iterative binary search. Avoids integer overflow
    with the (lo + (hi-lo)//2) idiom instead of (lo+hi)//2."""
    lo, hi = 0, len(arr) - 1
    while lo <= hi:
        mid = lo + (hi - lo) // 2
        if arr[mid] == target:
            return mid
        elif arr[mid] < target:
            lo = mid + 1
        else:
            hi = mid - 1
    return -1


def binary_search_left(arr, target):
    """Find leftmost index where arr[i] >= target (lower bound).
    Returns len(arr) if target is greater than all elements."""
    lo, hi = 0, len(arr)
    while lo < hi:
        mid = (lo + hi) // 2
        if arr[mid] < target:
            lo = mid + 1
        else:
            hi = mid
    return lo


def binary_search_right(arr, target):
    """Find rightmost insertion point (upper bound)."""
    lo, hi = 0, len(arr)
    while lo < hi:
        mid = (lo + hi) // 2
        if arr[mid] <= target:
            lo = mid + 1
        else:
            hi = mid
    return lo


def search_rotated_sorted(arr, target):
    """Search a sorted array that has been rotated at some pivot.
    Key insight: one half is always sorted — identify which one
    and decide where to search in O(log n)."""
    lo, hi = 0, len(arr) - 1
    while lo <= hi:
        mid = (lo + hi) // 2
        if arr[mid] == target:
            return mid
        if arr[lo] <= arr[mid]:        # left half is sorted
            if arr[lo] <= target < arr[mid]:
                hi = mid - 1
            else:
                lo = mid + 1
        else:                          # right half is sorted
            if arr[mid] < target <= arr[hi]:
                lo = mid + 1
            else:
                hi = mid - 1
    return -1


def find_peak(arr):
    """Find any local maximum (a peak element) in O(log n).
    Works because there must be a peak in whichever half we
    go toward the larger neighbor."""
    lo, hi = 0, len(arr) - 1
    while lo < hi:
        mid = (lo + hi) // 2
        if arr[mid] < arr[mid + 1]:
            lo = mid + 1
        else:
            hi = mid
    return lo


# ── Binary Search on the Answer ─────────────────────────────────
def min_eating_speed(piles, h):
    """Koko eats bananas. What is the minimum speed k (bananas/hr)
    to eat all piles in h hours? Classic 'search on answer' pattern."""
    def can_finish(speed):
        return sum((p + speed - 1) // speed for p in piles) <= h

    lo, hi = 1, max(piles)
    while lo < hi:
        mid = (lo + hi) // 2
        if can_finish(mid):
            hi = mid
        else:
            lo = mid + 1
    return lo


# Example
arr = [1, 3, 5, 7, 9, 11, 13, 15, 17, 19]
print(binary_search(arr, 13))       # → 6
print(binary_search_left(arr, 8))   # → 4 (insertion point)
print(search_rotated_sorted([4,5,6,7,0,1,2], 0))  # → 4
print(find_peak([1,2,3,1]))         # → 2
Proof — Binary Search Correctness & O(log n)

Invariant: If target is in arr, it is always in arr[lo..hi]. Initially true (lo=0, hi=n-1). At each step, we check arr[mid] and discard the half that cannot contain target — preserving the invariant. When lo > hi the array is empty so target is absent.

Termination: The range [lo, hi] has size (hi − lo + 1). At each step, mid ≠ lo and mid ≠ hi when lo < hi, so the range strictly shrinks by at least 1 each iteration. Therefore it reaches 0 in finite steps.

Complexity: Each iteration halves the search space. Starting from n, after k iterations the space is n/2^k. We stop when n/2^k = 1, so k = log₂(n). Thus exactly ⌊log₂(n)⌋ + 1 iterations — O(log n).

Merge Sort — O(n log n) with Inversion Counting

Python · Merge Sort & Inversion Count
O(n log n) O(n) space
def merge_sort(arr):
    """Stable, O(n log n) sort via divide-and-conquer."""
    if len(arr) <= 1:
        return arr[:]
    mid = len(arr) // 2
    left  = merge_sort(arr[:mid])
    right = merge_sort(arr[mid:])
    return _merge(left, right)

def _merge(left, right):
    result = []
    i = j = 0
    while i < len(left) and j < len(right):
        if left[i] <= right[j]:   # <= keeps sort stable
            result.append(left[i]); i += 1
        else:
            result.append(right[j]); j += 1
    result.extend(left[i:])
    result.extend(right[j:])
    return result


def count_inversions(arr):
    """Count pairs (i,j) where i < j but arr[i] > arr[j].
    An inversion measures 'sortedness'. Exploits the merge step:
    whenever we pick right[j] over left[i], all remaining left
    elements (len(left)-i of them) form inversions with right[j].
    Total time: O(n log n) — same as merge sort."""
    if len(arr) <= 1:
        return arr[:], 0
    mid = len(arr) // 2
    left,  left_inv  = count_inversions(arr[:mid])
    right, right_inv = count_inversions(arr[mid:])

    merged = []
    inversions = left_inv + right_inv
    i = j = 0
    while i < len(left) and j < len(right):
        if left[i] <= right[j]:
            merged.append(left[i]); i += 1
        else:
            inversions += len(left) - i  # ← key insight
            merged.append(right[j]); j += 1
    merged.extend(left[i:])
    merged.extend(right[j:])
    return merged, inversions


def merge_sort_inplace(arr, lo=0, hi=None):
    """In-place merge sort — modifies arr directly."""
    if hi is None: hi = len(arr) - 1
    if lo >= hi: return
    mid = (lo + hi) // 2
    merge_sort_inplace(arr, lo, mid)
    merge_sort_inplace(arr, mid + 1, hi)
    _merge_inplace(arr, lo, mid, hi)

def _merge_inplace(arr, lo, mid, hi):
    left  = arr[lo:mid+1]
    right = arr[mid+1:hi+1]
    i = j = 0; k = lo
    while i < len(left) and j < len(right):
        if left[i] <= right[j]:
            arr[k] = left[i]; i += 1
        else:
            arr[k] = right[j]; j += 1
        k += 1
    while i < len(left):  arr[k] = left[i];  i += 1; k += 1
    while j < len(right): arr[k] = right[j]; j += 1; k += 1


# Examples
print(merge_sort([5,2,8,1,9,3]))      # [1,2,3,5,8,9]
_, inv = count_inversions([3,1,2,4])
print(f"Inversions: {inv}")            # 2 → (3,1),(3,2)

Karatsuba Multiplication — O(n1.585)

Schoolbook long multiplication takes O(n²) digit operations. Karatsuba reduces this to just 3 recursive multiplications instead of 4, achieving O(n^log₂3) ≈ O(n^1.585).

Python · Karatsuba Multiplication
O(n^1.585)
def karatsuba(x, y):
    """Multiply two integers using Karatsuba's algorithm.
    Splits each number at midpoint m:
      x = x_hi * 10^m + x_lo
      y = y_hi * 10^m + y_lo
    Normally needs 4 multiplications. Karatsuba's trick:
      z0 = x_lo * y_lo
      z2 = x_hi * y_hi
      z1 = (x_lo+x_hi)(y_lo+y_hi) - z0 - z2   ← 3 muls total!
    Result: z2*10^2m + z1*10^m + z0"""
    if x < 10 or y < 10:
        return x * y

    n = max(len(str(x)), len(str(y)))
    m = n // 2

    x_hi, x_lo = divmod(x, 10**m)
    y_hi, y_lo = divmod(y, 10**m)

    z0 = karatsuba(x_lo, y_lo)
    z2 = karatsuba(x_hi, y_hi)
    z1 = karatsuba(x_lo + x_hi, y_lo + y_hi) - z0 - z2

    return z2 * 10**(2*m) + z1 * 10**m + z0


# Verify correctness against Python's built-in *
a, b = 123456789, 987654321
assert karatsuba(a, b) == a * b
print(f"{a} × {b} = {karatsuba(a,b)}")
# 123456789 × 987654321 = 121932631112635269
🔑 Karatsuba Recurrence

T(n) = 3T(n/2) + Θ(n) — By Master Theorem Case 1 (c* = log₂3 ≈ 1.585 > 1):

T(n) = Θ(n^log₂(3)) ≈ Θ(n^1.585)

vs. schoolbook multiplication which is Θ(n²). For 1000-digit numbers: roughly 10^6 vs 10^4.8 ≈ 63,000 — a 16× speedup.

Interactive — Binary Search Step-Through
Click ▶ Search to begin.
§02

Greedy Algorithms

locally optimal → globally optimal

The Exchange Argument Proof Technique

Most greedy algorithm proofs use the exchange argument: assume an optimal solution OPT differs from the greedy solution G. Show you can "swap" the differing element so OPT becomes more like G without getting worse — eventually showing G ≥ OPT, so G is optimal.

Proof Template — Exchange Argument

1. Let G = greedy solution, OPT = any optimal solution.

2. Find the first position where they differ.

3. Show swapping OPT's choice for G's choice does not decrease the objective.

4. Repeat until OPT = G. Therefore G is also optimal.

Activity Selection Problem

Given n activities with start/finish times, select the maximum number of non-overlapping activities. Greedy: always pick the activity with the earliest finish time.

Python · Activity Selection + Weighted Variant
O(n log n)
def activity_selection(activities):
    """activities = [(start, finish, name), ...]
    Returns max set of non-overlapping activities.
    Proof: Greedy-choice property — always safe to pick
    the activity finishing earliest, as it leaves the most
    room for future activities."""
    activities = sorted(activities, key=lambda x: x[1])  # sort by finish
    selected = [activities[0]]
    last_finish = activities[0][1]

    for start, finish, name in activities[1:]:
        if start >= last_finish:        # compatible
            selected.append((start, finish, name))
            last_finish = finish

    return selected


def weighted_activity_selection(activities):
    """When activities have weights (profits), greedy fails.
    Use DP: dp[i] = max profit using activities 0..i.
    O(n log n) with binary search for latest non-overlapping."""
    import bisect
    n = len(activities)
    acts = sorted(enumerate(activities), key=lambda x: x[1][1])
    finish_times = [a[1] for _, a in acts]

    dp = [0] * (n + 1)
    for i in range(n):
        _, (start, finish, weight) = acts[i]
        # Find last activity that finishes <= start
        idx = bisect.bisect_right(finish_times, start, 0, i)
        dp[i+1] = max(dp[i], dp[idx] + weight)

    return dp[n]


def fractional_knapsack(capacity, items):
    """items = [(value, weight), ...]
    Take fractions of items. Greedy: highest value/weight ratio first.
    This WORKS for fractional but NOT for 0-1 knapsack."""
    items = sorted(items, key=lambda x: x[0]/x[1], reverse=True)
    total = 0.0
    for value, weight in items:
        if capacity >= weight:
            total    += value
            capacity -= weight
        else:
            total += value * (capacity / weight)
            break
    return total


def scheduling_minimize_lateness(jobs):
    """jobs = [(processing_time, deadline), ...]
    Schedule to minimize maximum lateness. Greedy: earliest deadline first.
    Proof by exchange: swapping any two adjacent jobs in wrong order
    can only decrease max lateness."""
    jobs = sorted(jobs, key=lambda x: x[1])  # earliest deadline first
    time = 0
    schedule = []
    max_lateness = 0
    for p, d in jobs:
        time += p
        lateness = max(0, time - d)
        max_lateness = max(max_lateness, lateness)
        schedule.append((time - p, time, lateness))
    return schedule, max_lateness


# Examples
acts = [(1,4,'A'),(3,5,'B'),(0,6,'C'),(5,7,'D'),(3,8,'E'),(5,9,'F'),(6,10,'G'),(8,11,'H'),(8,12,'I'),(2,13,'J'),(12,14,'K')]
sel = activity_selection(acts)
print("Selected:", [a[2] for a in sel])   # A B D H K

items = [(60,10),(100,20),(120,30)]
print(f"Fractional knapsack (50kg): {fractional_knapsack(50, items):.1f}")  # 240.0
Proof — Activity Selection Correctness

Claim: Always choosing the activity with the earliest finish time yields the maximum number of non-overlapping activities.

Proof by exchange argument: Let G = {g₁, g₂, …, gₖ} be the greedy solution sorted by finish time, and OPT = {o₁, o₂, …, oₘ} any optimal solution. We show k = m.

Induction: We claim that for each i, finish(gᵢ) ≤ finish(oᵢ). Base case: g₁ has the earliest finish of all activities, so finish(g₁) ≤ finish(o₁). Inductive step: given finish(gᵢ) ≤ finish(oᵢ), activity oᵢ₊₁ starts after finish(oᵢ) ≤ finish(gᵢ), so oᵢ₊₁ is available to the greedy at step i+1. Greedy picks gᵢ₊₁ with minimum finish, so finish(gᵢ₊₁) ≤ finish(oᵢ₊₁). By induction, greedy selects at least as many activities — k ≥ m. Since OPT is optimal m ≥ k, so k = m. ∎

Huffman Coding — Optimal Prefix-Free Codes

Python · Huffman Encoding (Full Implementation)
O(n log n)
import heapq
from collections import Counter

class HNode:
    __slots__ = ('char', 'freq', 'left', 'right')
    def __init__(self, char, freq, left=None, right=None):
        self.char  = char
        self.freq  = freq
        self.left  = left
        self.right = right
    def __lt__(self, other):
        return self.freq < other.freq  # min-heap by frequency


def build_huffman_tree(text):
    """Build optimal prefix-free code tree via greedy bottom-up.
    Always merge the two lowest-frequency trees — proven optimal."""
    freq = Counter(text)
    if len(freq) == 1:
        char = next(iter(freq))
        return HNode(char, freq[char])

    heap = [HNode(c, f) for c, f in freq.items()]
    heapq.heapify(heap)

    while len(heap) > 1:
        a = heapq.heappop(heap)
        b = heapq.heappop(heap)
        heapq.heappush(heap, HNode(None, a.freq + b.freq, a, b))

    return heap[0]


def _build_codes(node, prefix='', codes=None):
    if codes is None: codes = {}
    if node.char is not None:          # leaf
        codes[node.char] = prefix or '0'
    else:
        _build_codes(node.left,  prefix + '0', codes)
        _build_codes(node.right, prefix + '1', codes)
    return codes


def huffman_encode(text):
    root  = build_huffman_tree(text)
    codes = _build_codes(root)
    encoded = ''.join(codes[c] for c in text)
    return encoded, codes, root


def huffman_decode(encoded, root):
    result = []
    node = root
    for bit in encoded:
        node = node.left if bit == '0' else node.right
        if node.char is not None:
            result.append(node.char)
            node = root
    return ''.join(result)


def compression_stats(text, codes):
    original_bits = len(text) * 8
    encoded_bits  = sum(len(codes[c]) * text.count(c) for c in set(text))
    ratio = encoded_bits / original_bits
    return {
        'original_bits':  original_bits,
        'encoded_bits':   encoded_bits,
        'ratio':          ratio,
        'savings_pct':    (1 - ratio) * 100,
    }


# Example
text  = "huffman coding is a greedy algorithm for optimal prefix codes"
enc, codes, root = huffman_encode(text)
stats = compression_stats(text, codes)
print(f"Original: {stats['original_bits']} bits")
print(f"Encoded:  {stats['encoded_bits']} bits")
print(f"Savings:  {stats['savings_pct']:.1f}%")
print(f"Decoded matches: {huffman_decode(enc, root) == text}")

# Show codebook
for char in sorted(codes, key=lambda c: len(codes[c])):
    print(f"  {repr(char):4s} → {codes[char]}")

Kruskal's Minimum Spanning Tree

Python · Kruskal's MST + Prim's MST
O(E log E)
class UnionFind:
    """Disjoint Set Union with path compression + union by rank.
    Both find() and union() run in amortized O(α(n)) ≈ O(1)."""
    def __init__(self, n):
        self.parent = list(range(n))
        self.rank   = [0] * n

    def find(self, x):
        if self.parent[x] != x:
            self.parent[x] = self.find(self.parent[x])  # path compression
        return self.parent[x]

    def union(self, x, y):
        rx, ry = self.find(x), self.find(y)
        if rx == ry: return False
        if self.rank[rx] < self.rank[ry]: rx, ry = ry, rx
        self.parent[ry] = rx
        if self.rank[rx] == self.rank[ry]: self.rank[rx] += 1
        return True


def kruskal_mst(n, edges):
    """Kruskal's: sort edges by weight, add each if it doesn't
    form a cycle (checked by Union-Find). Greedy proof: the
    'cycle property' — the heaviest edge in any cycle is never
    in an MST. We safely skip those edges."""
    edges = sorted(edges)   # sort by (weight, u, v)
    uf    = UnionFind(n)
    mst   = []
    total = 0
    for weight, u, v in edges:
        if uf.union(u, v):
            mst.append((u, v, weight))
            total += weight
            if len(mst) == n - 1: break   # MST complete
    return mst, total


def prim_mst(n, adj):
    """Prim's: grow MST from vertex 0 using a min-heap.
    adj[u] = [(weight, v), ...]. O((V+E) log V)."""
    import heapq
    visited = [False] * n
    heap    = [(0, 0, -1)]   # (weight, vertex, from)
    mst     = []
    total   = 0

    while heap and len(mst) < n:
        w, u, parent = heapq.heappop(heap)
        if visited[u]: continue
        visited[u] = True
        if parent != -1:
            mst.append((parent, u, w))
            total += w
        for wt, v in adj[u]:
            if not visited[v]:
                heapq.heappush(heap, (wt, v, u))

    return mst, total


# Example: 5-node graph
edges = [(1,0,1),(3,0,2),(6,0,3),(5,1,2),(4,1,3),(2,2,4),(6,3,4)]
mst, cost = kruskal_mst(5, edges)
print(f"MST edges: {mst}, total weight: {cost}")
Interactive — Huffman Encoder
Enter text and click Encode.
§03

Backtracking

prune the state-space tree

General Backtracking Template

Backtracking systematically explores candidates for a solution and abandons (backtracks) as soon as it determines a partial candidate cannot lead to a valid solution. It is depth-first search on an implicit state-space tree, with pruning.

Python · Universal Backtracking Pattern
O(b^d)
def backtrack_template(state, choices, results):
    """
    state   : current partial solution
    choices : available candidates at this step
    results : accumulator for complete solutions

    Pattern:
      1. Base case:  if solution complete → record it
      2. Loop over candidates:
           a. Prune:  skip invalid candidates early
           b. Choose: extend state with candidate
           c. Recurse
           d. Undo:   restore state (backtrack)
    """
    if is_complete(state):
        results.append(state.copy())
        return
    for candidate in get_candidates(state, choices):
        if is_valid(state, candidate):  # prune
            make_choice(state, candidate)
            backtrack_template(state, choices, results)
            undo_choice(state, candidate)   # ← critical backtrack step


# ── Concrete example: generate all subsets ────────────────────
def subsets(nums):
    """Power set of nums. 2^n subsets total."""
    result  = []
    current = []
    def bt(start):
        result.append(current[:])           # every partial is valid
        for i in range(start, len(nums)):
            current.append(nums[i])
            bt(i + 1)
            current.pop()                   # backtrack
    bt(0)
    return result

# ── Combination sum: reach target with repeats ───────────────
def combination_sum(candidates, target):
    """Find all combinations that sum to target (reuse allowed)."""
    result = []
    candidates.sort()
    def bt(start, current, remaining):
        if remaining == 0:
            result.append(current[:])
            return
        for i in range(start, len(candidates)):
            c = candidates[i]
            if c > remaining: break         # pruning: sorted, so stop early
            current.append(c)
            bt(i, current, remaining - c)   # i not i+1: allow reuse
            current.pop()
    bt(0, [], target)
    return result

# ── Word Search in a grid ─────────────────────────────────────
def word_search(board, word):
    """DFS + backtracking to find word in grid. O(M*N*4^L)."""
    rows, cols = len(board), len(board[0])
    def dfs(r, c, idx):
        if idx == len(word): return True
        if r < 0 or r >= rows or c < 0 or c >= cols: return False
        if board[r][c] != word[idx]: return False
        tmp, board[r][c] = board[r][c], '#'   # mark visited
        found = any(dfs(r+dr, c+dc, idx+1)
                    for dr, dc in [(0,1),(0,-1),(1,0),(-1,0)])
        board[r][c] = tmp                      # restore (backtrack)
        return found
    return any(dfs(r, c, 0)
               for r in range(rows)
               for c in range(cols))

# Examples
print(subsets([1,2,3]))
# [[], [1], [1,2], [1,2,3], [1,3], [2], [2,3], [3]]

print(combination_sum([2,3,6,7], 7))
# [[2,2,3], [7]]

N-Queens — Constraint-Based Backtracking

Python · N-Queens (Bitmask + Sets)
O(n!)
def solve_n_queens(n):
    """Place n non-attacking queens on an n×n board.
    Track three sets of forbidden positions:
      cols  : column already occupied
      diag1 : (row - col) diagonal occupied  [↘ direction]
      diag2 : (row + col) diagonal occupied  [↙ direction]
    For any (r,c): if (r-c) same → same ↘ diagonal.
                   if (r+c) same → same ↙ diagonal."""
    solutions = []
    cols, diag1, diag2 = set(), set(), set()
    board = [0] * n    # board[row] = column of queen

    def bt(row):
        if row == n:
            solutions.append(board[:])
            return
        for col in range(n):
            if col in cols or (row-col) in diag1 or (row+col) in diag2:
                continue                       # attacked — prune
            cols.add(col); diag1.add(row-col); diag2.add(row+col)
            board[row] = col
            bt(row + 1)
            cols.discard(col); diag1.discard(row-col); diag2.discard(row+col)

    bt(0)
    return solutions


def format_board(solution):
    n = len(solution)
    lines = []
    for r in range(n):
        row = '·' * n
        row = row[:solution[r]] + 'Q' + row[solution[r]+1:]
        lines.append(row)
    return '\n'.join(lines)


def solve_n_queens_bitmask(n):
    """Fastest bitmask solution: uses integer bit-ops for O(1) conflict checks.
    cols, d1, d2 are n-bit integers. Each call tries valid columns in one pass."""
    solutions = []
    full = (1 << n) - 1  # all columns filled

    def bt(cols, d1, d2, board):
        if cols == full:
            solutions.append(board)
            return
        available = full & ~(cols | d1 | d2)
        while available:
            bit = available & (-available)   # lowest available bit
            col = bit.bit_length() - 1
            board.append(col)
            bt(cols | bit,
               (d1 | bit) << 1,
               (d2 | bit) >> 1,
               board)
            board.pop()
            available &= available - 1       # clear that bit

    bt(0, 0, 0, [])
    return solutions


# N solutions count:  N=4→2, N=8→92, N=12→14200
for n in range(1, 10):
    print(f"N={n}: {len(solve_n_queens_bitmask(n))} solutions")

Sudoku Solver — Constraint Propagation + Backtracking

Python · Sudoku Solver
O(9^81) worst case, fast in practice
def solve_sudoku(board):
    """board is a 9×9 list. 0 = empty. Modifies in place.
    Uses constraint sets for O(1) validity checks.
    Returns True if solved, False if no solution."""
    rows = [set() for _ in range(9)]
    cols = [set() for _ in range(9)]
    boxs = [set() for _ in range(9)]

    # Initialize constraints
    for r in range(9):
        for c in range(9):
            d = board[r][c]
            if d:
                b = (r//3)*3 + c//3
                rows[r].add(d); cols[c].add(d); boxs[b].add(d)

    def find_most_constrained():
        """MRV heuristic: pick empty cell with fewest legal values.
        Dramatically reduces backtracking in practice."""
        best = None; best_count = 10
        for r in range(9):
            for c in range(9):
                if board[r][c] == 0:
                    b = (r//3)*3 + c//3
                    used = rows[r] | cols[c] | boxs[b]
                    cnt  = sum(1 for d in range(1,10) if d not in used)
                    if cnt < best_count:
                        best_count = cnt
                        best = (r, c)
                        if cnt == 0: return None, None   # dead end
        return best

    def bt():
        pos = find_most_constrained()
        if pos == (None, None): return False  # dead end
        r, c = pos
        if r is None: return True             # all filled → solved

        b = (r//3)*3 + c//3
        used = rows[r] | cols[c] | boxs[b]

        for d in range(1, 10):
            if d in used: continue
            # place
            board[r][c] = d
            rows[r].add(d); cols[c].add(d); boxs[b].add(d)
            if bt(): return True
            # backtrack
            board[r][c] = 0
            rows[r].discard(d); cols[c].discard(d); boxs[b].discard(d)

        return False

    return bt()


# Example: world's hardest sudoku
hard = [
    [8,0,0,0,0,0,0,0,0],
    [0,0,3,6,0,0,0,0,0],
    [0,7,0,0,9,0,2,0,0],
    [0,5,0,0,0,7,0,0,0],
    [0,0,0,0,4,5,7,0,0],
    [0,0,0,1,0,0,0,3,0],
    [0,0,1,0,0,0,0,6,8],
    [0,0,8,5,0,0,0,1,0],
    [0,9,0,0,0,0,4,0,0],
]
solve_sudoku(hard)
for row in hard: print(row)
Interactive — N-Queens Visualizer
6
Pick N and click Solve.
§04

Number Theory Algorithms

gcd · primes · modular arithmetic

Euclidean Algorithm & GCD

Python · GCD Suite
O(log min(a,b))
def gcd(a, b):
    """Euclidean algorithm: gcd(a,b) = gcd(b, a mod b).
    Key insight: gcd doesn't change when we replace the larger
    number with the remainder — proved by Euclid c. 300 BCE."""
    while b:
        a, b = b, a % b
    return a

def lcm(a, b):
    """LCM via GCD: lcm(a,b) = a*b / gcd(a,b).
    Use a//gcd(a,b)*b to avoid integer overflow."""
    return a // gcd(a, b) * b

def extended_gcd(a, b):
    """Finds x, y such that a*x + b*y = gcd(a, b).
    This is Bézout's identity — x,y always exist.
    Returns (gcd, x, y). Useful for modular inverse."""
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return g, x, y

def mod_inverse(a, m):
    """Modular inverse of a mod m (exists iff gcd(a,m)=1).
    Finds x such that a*x ≡ 1 (mod m)."""
    g, x, _ = extended_gcd(a % m, m)
    if g != 1:
        raise ValueError(f"gcd({a},{m})={g} ≠ 1, no inverse")
    return x % m

def chinese_remainder(remainders, moduli):
    """Chinese Remainder Theorem: solve system x ≡ r_i (mod m_i)
    Requires moduli to be pairwise coprime.
    Returns unique solution x mod (product of moduli)."""
    M = 1
    for m in moduli: M *= m
    x = 0
    for r, m in zip(remainders, moduli):
        Mi = M // m
        x += r * Mi * mod_inverse(Mi, m)
    return x % M

# Examples
print(gcd(48, 18))                    # 6
print(extended_gcd(35, 64))           # (1, 11, -6) → 35*11 + 64*(-6) = 1
print(mod_inverse(3, 11))             # 4  (3*4 = 12 ≡ 1 mod 11)
print(chinese_remainder([2,3,2],[3,5,7]))  # 23
Proof — Euclidean Algorithm Correctness & O(log n)

Correctness: gcd(a, b) = gcd(b, a mod b). Proof: any common divisor d of a and b also divides a mod b (= a − ⌊a/b⌋·b), so it divides b and a mod b. Conversely, any common divisor of b and (a mod b) divides a. The set of common divisors is the same, so their GCDs are equal.

O(log min(a,b)) bound — Lamé's Theorem: After two iterations, the remainder strictly decreases: a mod b < a/2. Therefore every two steps halve the larger number, giving at most 2 log₂(min(a,b)) steps.

Sieve of Eratosthenes & Prime Utilities

Python · Prime Sieves & Factorization
O(n log log n)
def sieve(n):
    """Sieve of Eratosthenes: find all primes up to n.
    For each prime p, mark multiples starting from p² as composite.
    Use bytearray for 8× memory efficiency over a list of bools."""
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0] = is_prime[1] = 0
    for p in range(2, int(n**0.5) + 1):
        if is_prime[p]:
            is_prime[p*p::p] = bytearray(len(is_prime[p*p::p]))
    return [i for i, v in enumerate(is_prime) if v]


def segmented_sieve(lo, hi):
    """Find primes in [lo, hi] without sieving all of 0..hi.
    Uses small primes up to sqrt(hi) to sieve each segment.
    Memory: O(sqrt(hi)) instead of O(hi)."""
    limit = int(hi**0.5) + 1
    small_primes = sieve(limit)
    segment = bytearray([1]) * (hi - lo + 1)
    if lo == 1: segment[0] = 0
    if lo == 0: segment[0] = segment[1] = 0

    for p in small_primes:
        start = max(p*p, ((lo + p - 1) // p) * p)
        segment[start - lo::p] = bytearray(len(segment[start - lo::p]))

    return [lo + i for i, v in enumerate(segment) if v]


def prime_factorization(n):
    """Return list of (prime, exponent) pairs for n. O(sqrt(n))."""
    factors = []
    d = 2
    while d * d <= n:
        if n % d == 0:
            exp = 0
            while n % d == 0:
                n //= d; exp += 1
            factors.append((d, exp))
        d += 1
    if n > 1: factors.append((n, 1))
    return factors


def euler_totient(n):
    """φ(n) = count of integers 1..n coprime to n.
    Using inclusion-exclusion over prime factors."""
    result = n
    p = 2
    temp = n
    while p * p <= temp:
        if temp % p == 0:
            while temp % p == 0: temp //= p
            result -= result // p
        p += 1
    if temp > 1: result -= result // temp
    return result


# Fast modular exponentiation — O(log exp)
def pow_mod(base, exp, mod):
    """Right-to-left binary exponentiation (square-and-multiply).
    Proof: represents exp in binary: base^exp = base^(b_k·2^k + ... + b_0)
    = product of base^(2^i) for each set bit i."""
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:             # if lowest bit is set
            result = result * base % mod
        base = base * base % mod
        exp >>= 1               # shift to next bit
    return result


primes = sieve(50)
print(f"Primes ≤ 50: {primes}")
# [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]

print(prime_factorization(360))
# [(2,3),(3,2),(5,1)] → 360 = 2³·3²·5

print(f"φ(36) = {euler_totient(36)}")   # 12
print(f"2^31 mod 10^9+7 = {pow_mod(2,31,10**9+7)}")
📐 Sieve Complexity — O(n log log n)

The total work done marking composites is the sum over all primes p ≤ n of n/p:

W = Σ_{p prime, p ≤ n} n/p = n · Σ_{p prime, p ≤ n} 1/p ≈ n · ln(ln n)

By Mertens' second theorem, the sum of reciprocals of primes up to n is Θ(log log n), giving the total complexity of Θ(n log log n). In practice this is nearly linear — for n = 10⁷ it's about 3.9n, remarkably fast.

Miller-Rabin Probabilistic Primality Test

Python · Miller-Rabin (Deterministic for n < 3.3×10²⁴)
O(k log² n)
def miller_rabin(n, witnesses=None):
    """Miller-Rabin primality test. With specific witness sets,
    deterministic for all n up to certain bounds:
      witnesses {2,3,5,7,11,13}: n < 3,215,031,751
      witnesses {2,3,5,7,11,13,17,19,23,29,31,37}: n < 3.3 × 10^24

    Algorithm:
      Write n-1 = 2^r · d  (factor out 2s from n-1)
      For each witness a:
        Compute x = a^d mod n
        If x == 1 or x == n-1: continue (probably prime for this a)
        Square x up to r-1 times:
          If x becomes n-1: continue
          If x becomes 1 without hitting n-1: composite!
      If all witnesses pass: probably (or provably) prime
    """
    if n < 2:  return False
    if n == 2 or n == 3: return True
    if n % 2 == 0: return False

    # Write n-1 = 2^r * d
    r, d = 0, n - 1
    while d % 2 == 0:
        r += 1; d //= 2

    if witnesses is None:
        witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]

    for a in witnesses:
        if a >= n: continue
        x = pow_mod(a, d, n)
        if x == 1 or x == n - 1: continue
        for _ in range(r - 1):
            x = x * x % n
            if x == n - 1: break
        else:
            return False    # composite
    return True             # prime (deterministically for n < 3.3e24)


def next_prime(n):
    """Find next prime after n."""
    n += 1 if n % 2 == 0 else 2
    while not miller_rabin(n):
        n += 2
    return n


# Test
primes_to_check = [2, 17, 97, 104729, 15485863, 10**9+7, 2**31-1]
for p in primes_to_check:
    print(f"{p:15,d}  prime={miller_rabin(p)}")
Interactive — Sieve of Eratosthenes (2–120)
Click ▶ Animate Sieve or use Step.
§05

Randomized Algorithms

expected efficiency · Las Vegas · Monte Carlo

QuickSelect — Expected O(n) Selection

Python · QuickSelect & Variants
E[O(n)] O(1) space
import random

def quickselect(arr, k):
    """Find the kth smallest element (0-indexed) in expected O(n).
    Uses random pivot to avoid worst-case O(n²).
    Partitions array: elements ≤ pivot | pivot | elements > pivot
    Recurse only into the side containing rank k — no merging needed."""
    arr = arr[:]   # don't mutate input

    def _select(lo, hi, k):
        if lo == hi: return arr[lo]

        # Random pivot selection
        pivot_idx = random.randint(lo, hi)
        arr[pivot_idx], arr[hi] = arr[hi], arr[pivot_idx]
        pivot = arr[hi]

        # Three-way partition: [lo..i-1] ≤ pivot, i = pivot, [i+1..hi] > pivot
        i = lo
        for j in range(lo, hi):
            if arr[j] <= pivot:
                arr[i], arr[j] = arr[j], arr[i]
                i += 1
        arr[i], arr[hi] = arr[hi], arr[i]

        if i == k:        return arr[i]
        elif i < k:       return _select(i + 1, hi, k)
        else:             return _select(lo, i - 1, k)

    return _select(0, len(arr) - 1, k)


def nth_element(arr, k):
    """Partial sort: rearrange so arr[k] is the kth smallest,
    elements before it are ≤ arr[k], elements after ≥ arr[k].
    Useful for median-finding, top-k problems."""
    arr = arr[:]
    def _partition(lo, hi):
        pivot_idx = random.randint(lo, hi)
        arr[pivot_idx], arr[hi] = arr[hi], arr[pivot_idx]
        i = lo
        for j in range(lo, hi):
            if arr[j] <= arr[hi]:
                arr[i], arr[j] = arr[j], arr[i]; i += 1
        arr[i], arr[hi] = arr[hi], arr[i]
        return i
    lo, hi = 0, len(arr) - 1
    while lo < hi:
        pivot = _partition(lo, hi)
        if pivot < k:   lo = pivot + 1
        elif pivot > k: hi = pivot - 1
        else:           break
    return arr


def reservoir_sample(stream, k):
    """Sample k items uniformly from a stream of unknown size.
    After seeing n items, each has probability k/n of being chosen.
    Only requires O(k) memory regardless of stream size."""
    reservoir = []
    for i, item in enumerate(stream):
        if i < k:
            reservoir.append(item)
        else:
            j = random.randint(0, i)
            if j < k:
                reservoir[j] = item   # replace with probability k/(i+1)
    return reservoir


def randomized_quick_sort(arr):
    """Average O(n log n) via random pivot. Expected comparisons:
    2n ln n ≈ 1.386 n log₂n — within 39% of optimal."""
    if len(arr) <= 1: return arr
    pivot = random.choice(arr)
    less    = [x for x in arr if x <  pivot]
    equal   = [x for x in arr if x == pivot]
    greater = [x for x in arr if x >  pivot]
    return randomized_quick_sort(less) + equal + randomized_quick_sort(greater)


# Examples
arr = [7,2,1,6,3,9,4,8,5]
print(quickselect(arr, 4))              # 5 (5th smallest, 0-indexed)
print(quickselect(arr, 0))              # 1 (minimum)
print(quickselect(arr, len(arr)-1))     # 9 (maximum)
stream = range(1, 1001)
sample = reservoir_sample(stream, 5)
print(f"5 random from 1..1000: {sample}")
Proof — QuickSelect Expected O(n)

Let T(n) = expected comparisons on input of size n. With a random pivot landing at rank r (uniform in 0..n-1):

T(n) = (n - 1) + (1/n) · Σ_{r=0}^{n-1} T(max(r, n-1-r))

The max(r, n-1-r) is the size of the side we recurse into. Observe that half the time r is in [n/4, 3n/4] — a "good pivot" — which guarantees at most 3n/4 elements on each side. Let p_good = 1/2. Each "good" call reduces size by 1/4. The expected number of rounds until a good pivot is chosen is 2. So:

T(n) = n + T(3n/4) with probability 1/2 (geometric waiting for good pivot) ⟹ T(n) ≤ 4n → T(n) = O(n)

A more precise analysis by solving the recurrence yields T(n) ≤ 4n comparisons in expectation. ∎

Monte Carlo & Probabilistic Algorithms

Python · Monte Carlo Methods
O(n) samples
import random, math

def estimate_pi(n_samples, seed=None):
    """Monte Carlo Pi: sample (x,y) uniform in [0,1]².
    P(x²+y² ≤ 1) = π/4  (quarter circle area / unit square).
    By law of large numbers, frequency → π/4 as n → ∞.
    Error is O(1/√n) — need 10^6 samples for ~3 decimal places."""
    if seed: random.seed(seed)
    inside = sum(1 for _ in range(n_samples)
                 if random.random()**2 + random.random()**2 <= 1)
    return 4 * inside / n_samples


def monte_carlo_integration(f, a, b, n_samples=10**5):
    """Estimate ∫_a^b f(x) dx using Monte Carlo.
    E[f(U)] where U ~ Uniform[a,b] = (1/(b-a)) ∫_a^b f(x) dx
    So: integral ≈ (b-a) * mean(f(U_i))"""
    samples = [f(random.uniform(a, b)) for _ in range(n_samples)]
    return (b - a) * sum(samples) / n_samples


def randomized_load_balance(jobs, m):
    """Assign n jobs to m machines randomly. Expected makespan
    is O(max_job + Σjobs/m + log(n)/log(log(n))) with high probability."""
    machines = [0.0] * m
    assignments = []
    for job in jobs:
        machine = random.randrange(m)
        machines[machine] += job
        assignments.append(machine)
    return assignments, max(machines)


def karger_min_cut(graph, trials=50):
    """Karger's randomized min-cut algorithm. O(V² log V) total.
    Each single trial succeeds with probability ≥ 2/(V(V-1)).
    After O(V² log V) trials, error probability is 1/V."""
    import copy
    def contract(g):
        g = copy.deepcopy(g)
        vertices = list(g.keys())
        while len(vertices) > 2:
            u = random.choice(vertices)
            if not g[u]: continue
            v = random.choice(g[u])
            # Merge v into u
            for w in g[v]:
                if w != u:
                    g[u].append(w)
                    g[w].append(u)
                    g[w] = [x for x in g[w] if x != v]
            g[u] = [x for x in g[u] if x != v]
            del g[v]
            vertices = list(g.keys())
        u = vertices[0]
        return len(g[u])   # number of crossing edges

    graph_adj = {}
    for u, v in graph:
        graph_adj.setdefault(u, []).append(v)
        graph_adj.setdefault(v, []).append(u)

    return min(contract(graph_adj) for _ in range(trials))


# Examples
for n in [1000, 10000, 100000, 1000000]:
    pi_est = estimate_pi(n)
    print(f"n={n:8d}: π ≈ {pi_est:.6f}  error={abs(pi_est - math.pi):.6f}")

# ∫₀¹ x² dx = 1/3
est = monte_carlo_integration(lambda x: x**2, 0, 1)
print(f"∫₀¹ x² dx ≈ {est:.4f}  (true: 0.3333)")
Interactive — Monte Carlo π Estimator
π ≈ ?
Click ▶ Run to estimate π with random sampling.
§06

Bit Manipulation

O(1) tricks with binary representations

Essential Bit Operations Reference

OperationCodeResultUse Case
Get bit k(n >> k) & 10 or 1Check if bit k is set
Set bit kn | (1 << k)n with bit k = 1Enable a flag
Clear bit kn & ~(1 << k)n with bit k = 0Disable a flag
Toggle bit kn ^ (1 << k)n with bit k flippedSwitch state
Lowest set bitn & (-n)Isolated LSBFenwick tree, factor 2
Clear lowest bitn & (n-1)n with LSB removedCount bits (Kernighan)
Is power of 2n > 0 and (n & n-1) == 0boolPower detection
Round up to pow21 << n.bit_length()next power ≥ nBuffer sizing
XOR swapa^=b; b^=a; a^=ba,b swappedNo temp variable
Arithmetic meanlo + ((hi-lo)>>1)no overflowSafe binary search mid

Bit Counting, XOR Tricks & Bitmask DP

Python · Bit Manipulation Algorithms
O(1) or O(bits)
def count_bits_kernighan(n):
    """Brian Kernighan's trick: n & (n-1) clears the lowest set bit.
    Loop runs exactly k times where k = popcount(n). O(k)."""
    count = 0
    while n:
        n &= n - 1      # clear lowest set bit
        count += 1
    return count


def popcount_parallel(n):
    """SWAR (SIMD Within A Register) parallel bit count.
    Works in O(log bits) = O(5) steps for 32-bit integers.
    Used in hardware implementations."""
    n = n - ((n >> 1) & 0x55555555)
    n = (n & 0x33333333) + ((n >> 2) & 0x33333333)
    n = (n + (n >> 4)) & 0x0F0F0F0F
    return (n * 0x01010101) >> 24


def find_unique_xor(arr):
    """All elements appear twice except one. Find it in O(n), O(1) space.
    XOR is commutative, associative, and a^a=0, a^0=a.
    So all doubles cancel, leaving the unique element."""
    result = 0
    for x in arr: result ^= x
    return result


def find_two_unique_xor(arr):
    """Two elements appear once, rest appear twice.
    Step 1: XOR all → xor = a ^ b (both unique elements)
    Step 2: Find any set bit in xor — it differs between a and b
    Step 3: Split array by that bit → each group XORs to one answer."""
    xor = 0
    for x in arr: xor ^= x
    diff_bit = xor & (-xor)     # lowest differing bit
    a = b = 0
    for x in arr:
        if x & diff_bit: a ^= x
        else:            b ^= x
    return a, b


def subsets_bitmask(n):
    """Generate all subsets of {0,1,...,n-1} using bitmasks.
    Mask bit i set ↔ element i is in subset. 2^n total."""
    result = []
    for mask in range(1 << n):
        subset = [i for i in range(n) if mask >> i & 1]
        result.append(subset)
    return result


def tsp_dp_bitmask(dist):
    """Traveling Salesman Problem via Held-Karp algorithm.
    dp[mask][i] = min cost to visit exactly the cities in 'mask',
    ending at city i, starting from city 0.
    O(2^n · n²) time, O(2^n · n) space — exponential but tractable for n≤20."""
    n = len(dist)
    INF = float('inf')
    dp = [[INF] * n for _ in range(1 << n)]
    dp[1][0] = 0   # start at city 0, mask has only bit 0 set

    for mask in range(1 << n):
        for u in range(n):
            if not (mask >> u & 1) or dp[mask][u] == INF:
                continue
            for v in range(n):
                if mask >> v & 1: continue    # already visited
                new_mask = mask | (1 << v)
                cost = dp[mask][u] + dist[u][v]
                if cost < dp[new_mask][v]:
                    dp[new_mask][v] = cost

    full = (1 << n) - 1
    return min(dp[full][i] + dist[i][0] for i in range(1, n))


def gray_code(n):
    """Generate n-bit Gray code sequence (only 1 bit changes each step).
    Formula: gray(k) = k XOR (k >> 1)"""
    return [i ^ (i >> 1) for i in range(1 << n)]


# Examples
print(count_bits_kernighan(0b10110100))  # 4
print(find_unique_xor([4,1,2,1,2]))      # 4
a, b = find_two_unique_xor([1,2,1,3,2,5])
print(f"Unique pair: {a}, {b}")          # 3, 5
print(gray_code(3))
# [0, 1, 3, 2, 6, 7, 5, 4]
Interactive — Binary Bit Visualizer
§07

Computational Geometry

points · hulls · sweepline

Orientation Test — The Cross Product

The signed cross product of vectors OA and OB determines if three points make a left turn, right turn, or are collinear. It is the core primitive of all geometric algorithms.

Python · Geometry Primitives
O(1)
def cross(O, A, B):
    """2D cross product of vectors OA and OB.
    Returns:
      > 0: counterclockwise turn (A, B left of O→A)
      < 0: clockwise turn
      = 0: collinear
    Formula: (A-O) × (B-O) = (Ax-Ox)(By-Oy) - (Ay-Oy)(Bx-Ox)"""
    return (A[0]-O[0]) * (B[1]-O[1]) - (A[1]-O[1]) * (B[0]-O[0])


def dist_sq(p, q):
    """Squared Euclidean distance — avoids sqrt, exact for integers."""
    return (p[0]-q[0])**2 + (p[1]-q[1])**2


def segments_intersect(p1, p2, p3, p4):
    """Do segments p1p2 and p3p4 intersect?
    Uses orientation tests + collinear overlap check."""
    def on_segment(p, q, r):   # Is q on segment pr?
        return (min(p[0],r[0]) <= q[0] <= max(p[0],r[0]) and
                min(p[1],r[1]) <= q[1] <= max(p[1],r[1]))

    d1 = cross(p3, p4, p1)
    d2 = cross(p3, p4, p2)
    d3 = cross(p1, p2, p3)
    d4 = cross(p1, p2, p4)

    if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
       ((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0)):
        return True

    if d1 == 0 and on_segment(p3, p1, p4): return True
    if d2 == 0 and on_segment(p3, p2, p4): return True
    if d3 == 0 and on_segment(p1, p3, p2): return True
    if d4 == 0 and on_segment(p1, p4, p2): return True
    return False


def point_in_polygon(point, polygon):
    """Ray casting algorithm: cast ray from point to +∞ in x direction.
    Count crossings with polygon edges — odd count = inside.
    O(n) where n = number of polygon vertices."""
    x, y = point
    inside = False
    n = len(polygon)
    j = n - 1
    for i in range(n):
        xi, yi = polygon[i]
        xj, yj = polygon[j]
        if ((yi > y) != (yj > y)) and \
           x < (xj - xi) * (y - yi) / (yj - yi) + xi:
            inside = not inside
        j = i
    return inside


def polygon_area(polygon):
    """Area via Shoelace formula (Gauss's formula).
    Signed area — positive if CCW, negative if CW.
    |Area| = |Σ (x_i(y_{i+1} - y_{i-1})| / 2"""
    n = len(polygon)
    area = 0
    for i in range(n):
        j = (i + 1) % n
        area += polygon[i][0] * polygon[j][1]
        area -= polygon[j][0] * polygon[i][1]
    return abs(area) / 2


# Examples
O, A, B = (0,0), (1,0), (0,1)
print(cross(O, A, B))    # 1.0  → counterclockwise

poly = [(0,0),(4,0),(4,3),(0,3)]
print(point_in_polygon((2,1.5), poly))   # True
print(polygon_area(poly))                # 12.0

Graham Scan Convex Hull — O(n log n)

Python · Convex Hull (Graham Scan + Andrew's Monotone Chain)
O(n log n)
def convex_hull_monotone(points):
    """Andrew's Monotone Chain — simpler than Graham Scan, same complexity.
    Builds lower hull left-to-right, then upper hull right-to-left.
    Invariant: hull is always making left turns (CCW). If cross ≤ 0
    (right turn or collinear), pop the middle point — it's inside."""
    points = sorted(set(map(tuple, points)))   # remove duplicates, sort
    n = len(points)
    if n <= 1: return points

    def build_half(pts):
        hull = []
        for p in pts:
            while len(hull) >= 2 and cross(hull[-2], hull[-1], p) <= 0:
                hull.pop()
            hull.append(p)
        return hull

    lower = build_half(points)
    upper = build_half(reversed(points))
    return lower[:-1] + upper[:-1]   # omit last point (duplicate of first)


def convex_hull_graham(points):
    """Classic Graham Scan: sort by polar angle from lowest-leftmost point.
    Process in order, maintaining CCW invariant with a stack."""
    import math
    pts = [tuple(p) for p in points]
    n = len(pts)
    if n < 3: return pts

    # Find anchor: lowest then leftmost
    anchor = min(pts, key=lambda p: (p[1], p[0]))

    def polar_angle(p):
        dx, dy = p[0] - anchor[0], p[1] - anchor[1]
        return math.atan2(dy, dx)

    def dist(p):
        return dist_sq(p, anchor)

    pts.remove(anchor)
    pts.sort(key=lambda p: (polar_angle(p), dist(p)))

    # Remove collinear points at same angle (keep farthest)
    filtered = []
    i = 0
    while i < len(pts):
        j = i
        while j < len(pts) - 1 and polar_angle(pts[j]) == polar_angle(pts[j+1]):
            j += 1
        filtered.append(pts[j])
        i = j + 1

    if len(filtered) < 2: return [anchor] + filtered

    stack = [anchor, filtered[0], filtered[1]]
    for p in filtered[2:]:
        while len(stack) > 1 and cross(stack[-2], stack[-1], p) <= 0:
            stack.pop()
        stack.append(p)
    return stack


def closest_pair_of_points(points):
    """Divide-and-conquer closest pair. O(n log n).
    Key insight: in the strip of width 2δ around the dividing line,
    each point needs to check at most 7 other points."""
    import math
    pts = sorted(points)
    n   = len(pts)

    def brute(p):
        d = float('inf')
        pair = None
        for i in range(len(p)):
            for j in range(i+1, len(p)):
                dd = math.dist(p[i], p[j])
                if dd < d: d = dd; pair = (p[i], p[j])
        return d, pair

    def rec(p_sorted_x):
        n = len(p_sorted_x)
        if n <= 3: return brute(p_sorted_x)

        mid = n // 2
        mid_x = p_sorted_x[mid][0]

        dl, pl = rec(p_sorted_x[:mid])
        dr, pr = rec(p_sorted_x[mid:])
        d, best = (dl, pl) if dl <= dr else (dr, pr)

        strip = [p for p in p_sorted_x if abs(p[0] - mid_x) < d]
        strip.sort(key=lambda p: p[1])

        for i in range(len(strip)):
            for j in range(i+1, len(strip)):
                if strip[j][1] - strip[i][1] >= d: break
                dd = math.dist(strip[i], strip[j])
                if dd < d: d = dd; best = (strip[i], strip[j])

        return d, best

    return rec(pts)


# Examples
pts = [(0,0),(1,1),(2,2),(0,2),(2,0),(-1,1),(1,-1),(3,1)]
hull = convex_hull_monotone(pts)
print(f"Hull: {hull}")

d, pair = closest_pair_of_points(pts)
print(f"Closest pair: {pair}, distance: {d:.4f}")
Proof — Convex Hull Correctness (Monotone Chain)

Invariant: The lower hull array always forms a CCW chain. When adding a new point p, if the last turn is not strictly CCW (cross product ≤ 0), the middle point is not on the hull and is removed.

Correctness of popping: If cross(A, B, C) ≤ 0, then B is not a left-turn vertex — B is either inside the triangle AC-p or collinear. No convex hull vertex can be interior to another triangle formed by three points on the hull, so B is correctly excluded.

Completeness: No valid hull vertex is ever wrongly popped. If B were on the true hull, then cross(A,B,C) > 0, contradicting the pop condition. ∎

Interactive — Convex Hull (Click to Add Points)
Click canvas to add points
Click anywhere on the canvas to add a point.
§08

Advanced Data Structures

trie · segment tree · DSU · Fenwick tree

Trie (Prefix Tree)

Python · Trie — Full Implementation
O(L) insert/search O(L) space/word
class TrieNode:
    __slots__ = ('children', 'is_end', 'count')
    def __init__(self):
        self.children = {}
        self.is_end   = False
        self.count    = 0    # number of words passing through


class Trie:
    """Prefix tree for O(L) insert, search, prefix-count, and
    autocomplete where L = word length. Space: O(total_chars)."""

    def __init__(self):
        self.root = TrieNode()

    def insert(self, word):
        node = self.root
        for ch in word:
            if ch not in node.children:
                node.children[ch] = TrieNode()
            node = node.children[ch]
            node.count += 1
        node.is_end = True

    def search(self, word):
        """Return True if exact word exists."""
        node = self._walk(word)
        return node is not None and node.is_end

    def starts_with(self, prefix):
        """Return True if any word starts with prefix."""
        return self._walk(prefix) is not None

    def count_prefix(self, prefix):
        """Count words with this prefix."""
        node = self._walk(prefix)
        return node.count if node else 0

    def autocomplete(self, prefix, limit=10):
        """Return up to 'limit' words with given prefix."""
        node = self._walk(prefix)
        if not node: return []
        results = []
        def dfs(n, path):
            if len(results) >= limit: return
            if n.is_end: results.append(prefix[:-len(path)] + path if path else prefix)
            for ch, child in sorted(n.children.items()):
                dfs(child, path + ch)
        # Rebuild path from prefix node
        def dfs2(n, current):
            if len(results) >= limit: return
            if n.is_end: results.append(current)
            for ch, child in sorted(n.children.items()):
                dfs2(child, current + ch)
        dfs2(node, prefix)
        return results

    def delete(self, word):
        """Remove word from trie. Returns True if deleted."""
        def _del(node, word, idx):
            if idx == len(word):
                if not node.is_end: return False
                node.is_end = False
                return len(node.children) == 0
            ch = word[idx]
            if ch not in node.children: return False
            child = node.children[ch]
            child.count -= 1
            should_delete = _del(child, word, idx + 1)
            if should_delete:
                del node.children[ch]
                return not node.is_end and len(node.children) == 0
            return False
        _del(self.root, word, 0)

    def _walk(self, prefix):
        node = self.root
        for ch in prefix:
            if ch not in node.children: return None
            node = node.children[ch]
        return node


# Example
t = Trie()
words = ["apple","app","apply","application","apt","banana","band","bandana"]
for w in words: t.insert(w)

print(t.search("app"))           # True
print(t.search("ap"))            # False
print(t.count_prefix("app"))     # 4 (app, apple, apply, application)
print(t.autocomplete("app"))     # ['app','apple','application','apply']
t.delete("app")
print(t.search("app"))           # False
print(t.search("apple"))         # True (apple still exists)

Segment Tree — Range Queries & Point Updates

Python · Segment Tree with Lazy Propagation
O(log n) query/update O(n) space
class SegTree:
    """Segment tree supporting:
      - Range sum / min / max queries: O(log n)
      - Point update: O(log n)
      - Range update (with lazy propagation): O(log n)
    Internal representation: 1-indexed array of size 4n."""

    def __init__(self, arr):
        self.n    = len(arr)
        self.tree = [0] * (4 * self.n)
        self.lazy = [0] * (4 * self.n)
        self._build(arr, 1, 0, self.n - 1)

    def _build(self, arr, node, lo, hi):
        if lo == hi:
            self.tree[node] = arr[lo]; return
        mid = (lo + hi) // 2
        self._build(arr, 2*node,   lo,    mid)
        self._build(arr, 2*node+1, mid+1, hi)
        self.tree[node] = self.tree[2*node] + self.tree[2*node+1]

    def _push_down(self, node, lo, hi):
        """Propagate lazy updates to children."""
        if self.lazy[node]:
            mid = (lo + hi) // 2
            self._apply(2*node,   lo,    mid, self.lazy[node])
            self._apply(2*node+1, mid+1, hi,  self.lazy[node])
            self.lazy[node] = 0

    def _apply(self, node, lo, hi, delta):
        self.tree[node] += delta * (hi - lo + 1)
        self.lazy[node] += delta

    def range_update(self, l, r, delta, node=1, lo=0, hi=None):
        """Add delta to all elements in [l, r]. O(log n)."""
        if hi is None: hi = self.n - 1
        if r < lo or hi < l: return
        if l <= lo and hi <= r:
            self._apply(node, lo, hi, delta); return
        self._push_down(node, lo, hi)
        mid = (lo + hi) // 2
        self.range_update(l, r, delta, 2*node,   lo,    mid)
        self.range_update(l, r, delta, 2*node+1, mid+1, hi)
        self.tree[node] = self.tree[2*node] + self.tree[2*node+1]

    def range_query(self, l, r, node=1, lo=0, hi=None):
        """Sum of elements in [l, r]. O(log n)."""
        if hi is None: hi = self.n - 1
        if r < lo or hi < l: return 0
        if l <= lo and hi <= r: return self.tree[node]
        self._push_down(node, lo, hi)
        mid = (lo + hi) // 2
        return (self.range_query(l, r, 2*node,   lo,    mid) +
                self.range_query(l, r, 2*node+1, mid+1, hi))

    def point_update(self, idx, val, node=1, lo=0, hi=None):
        """Set element at idx to val. O(log n)."""
        if hi is None: hi = self.n - 1
        if lo == hi:
            self.tree[node] = val; return
        self._push_down(node, lo, hi)
        mid = (lo + hi) // 2
        if idx <= mid: self.point_update(idx, val, 2*node,   lo,    mid)
        else:          self.point_update(idx, val, 2*node+1, mid+1, hi)
        self.tree[node] = self.tree[2*node] + self.tree[2*node+1]


# Example
st = SegTree([1, 3, 5, 7, 9, 11])
print(st.range_query(1, 3))        # 3+5+7 = 15
st.range_update(1, 3, 2)           # add 2 to positions 1-3
print(st.range_query(1, 3))        # 5+7+9 = 21
print(st.range_query(0, 5))        # 1+5+7+9+9+11 = 42

Disjoint Set Union — Path Compression + Union by Rank

Python · DSU / Union-Find
O(α(n)) per op O(n) space
class DSU:
    """Disjoint Set Union (Union-Find) with:
    - Path compression in find(): amortizes to O(α(n)) per call
    - Union by rank: keeps trees flat
    Combined: nearly O(1) per operation. α(n) ≤ 4 for all practical n.

    Applications: Kruskal's MST, connected components, cycle detection,
    network connectivity, image segmentation, percolation."""

    def __init__(self, n):
        self.parent    = list(range(n))
        self.rank      = [0] * n
        self.size      = [1] * n     # size of each component
        self.num_comp  = n           # number of components

    def find(self, x):
        """Find root with full path compression (two-pass)."""
        root = x
        while self.parent[root] != root:
            root = self.parent[root]
        while self.parent[x] != root:       # path compression
            self.parent[x], x = root, self.parent[x]
        return root

    def union(self, x, y):
        """Union by rank. Returns True if merged, False if same component."""
        rx, ry = self.find(x), self.find(y)
        if rx == ry: return False
        # Attach smaller-rank tree under larger-rank
        if self.rank[rx] < self.rank[ry]: rx, ry = ry, rx
        self.parent[ry]  = rx
        self.size[rx]   += self.size[ry]
        if self.rank[rx] == self.rank[ry]: self.rank[rx] += 1
        self.num_comp   -= 1
        return True

    def connected(self, x, y):  return self.find(x) == self.find(y)
    def comp_size(self, x):      return self.size[self.find(x)]


def count_connected_components(n, edges):
    dsu = DSU(n)
    for u, v in edges: dsu.union(u, v)
    return dsu.num_comp


def detect_cycle_undirected(n, edges):
    """A cycle exists iff we ever try to union two vertices
    already in the same component."""
    dsu = DSU(n)
    for u, v in edges:
        if not dsu.union(u, v): return True   # already connected → cycle
    return False


def redundant_connection(edges):
    """Given n nodes and n edges forming exactly one cycle,
    find and remove the last edge that creates the cycle."""
    n = max(max(e) for e in edges) + 1
    dsu = DSU(n)
    for u, v in edges:
        if not dsu.union(u, v): return (u, v)
    return None


# Examples
edges = [(0,1),(1,2),(2,3),(3,4)]
print(count_connected_components(6, edges))   # 2 (node 5 isolated)
print(detect_cycle_undirected(4, [(0,1),(1,2),(2,0)]))  # True
print(redundant_connection([(1,2),(1,3),(2,3)]))  # (2,3)
Proof — Path Compression + Union by Rank Gives O(α(n))

Rank bound: A tree of rank k has at least 2^k nodes. Proof by induction: rank increases only when two trees of equal rank k-1 are merged, giving ≥ 2·2^(k-1) = 2^k nodes. So rank ≤ log₂(n).

Amortized analysis (Tarjan 1975): Define a potential function Φ based on rank gaps along paths. Path compression "flattens" the tree, releasing potential that pays for the compression work. The amortized cost per operation is O(log*(n)) without rank, and O(α(n)) with both optimizations, where α is the inverse Ackermann function.

Practical bound: α(2^65536) = 5. For all inputs that exist in the physical universe, α(n) ≤ 4, making each operation effectively O(1). ∎

Fenwick Tree (Binary Indexed Tree)

Python · Fenwick Tree — Prefix Sums & Inversions
O(log n) update/query O(n) space
class Fenwick:
    """Binary Indexed Tree (Fenwick 1994).
    Stores partial sums in a 1-indexed array using bit magic.
    Each node i is 'responsible' for the range (i - LSB(i), i]
    where LSB(i) = i & (-i) (lowest set bit of i).

    tree[i] = sum of arr[i - LSB(i) + 1 .. i]

    To update index i: add to tree[i], tree[i+LSB(i)], ...
    To prefix-sum to i: sum tree[i], tree[i-LSB(i)], ..."""

    def __init__(self, n):
        self.n    = n
        self.tree = [0] * (n + 1)

    def update(self, i, delta):
        """Add delta to position i (1-indexed). O(log n)."""
        while i <= self.n:
            self.tree[i] += delta
            i += i & (-i)          # move to next responsible node

    def prefix_sum(self, i):
        """Sum from position 1 to i (inclusive). O(log n)."""
        total = 0
        while i > 0:
            total += self.tree[i]
            i -= i & (-i)          # move to parent node
        return total

    def range_sum(self, l, r):
        """Sum from l to r (inclusive). O(log n)."""
        return self.prefix_sum(r) - self.prefix_sum(l - 1)

    def find_kth(self, k):
        """Find smallest index with prefix_sum ≥ k. O(log n).
        Uses binary lifting — walks down bit by bit."""
        pos = 0
        log = self.n.bit_length()
        for i in range(log, -1, -1):
            nxt = pos + (1 << i)
            if nxt <= self.n and self.tree[nxt] < k:
                pos  = nxt
                k   -= self.tree[nxt]
        return pos + 1


def count_inversions_fenwick(arr):
    """Count inversions in O(n log n) using coordinate compression
    and a Fenwick tree. An inversion at position i is the number
    of elements already processed that are > arr[i]."""
    # Coordinate compress to 1..n
    sorted_vals = sorted(set(arr))
    compress    = {v: i+1 for i, v in enumerate(sorted_vals)}
    n           = len(arr)
    fw          = Fenwick(n)

    inversions = 0
    for i, x in enumerate(arr):
        cx = compress[x]
        # Elements already inserted that are > x
        inversions += i - fw.prefix_sum(cx)
        fw.update(cx, 1)

    return inversions


# Range minimum query via sparse table — O(n log n) build, O(1) query
def build_sparse_table(arr):
    """Sparse table for static RMQ. Preprocessing O(n log n).
    Query O(1) using overlap of two ranges of equal power-of-2 length."""
    import math
    n  = len(arr)
    k  = max(1, int(math.log2(n)) + 1)
    st = [[float('inf')] * n for _ in range(k)]
    st[0] = arr[:]
    for j in range(1, k):
        for i in range(n - (1 << j) + 1):
            st[j][i] = min(st[j-1][i], st[j-1][i + (1 << (j-1))])
    return st

def rmq_query(st, l, r):
    """O(1) range minimum query using sparse table."""
    import math
    j = int(math.log2(r - l + 1))
    return min(st[j][l], st[j][r - (1 << j) + 1])


# Examples
fw = Fenwick(10)
for v in [3, 1, 6, 7, 2]: fw.update(v, 1)
print(fw.prefix_sum(5))    # 3 elements ≤ 5: {1,2,3}
print(fw.range_sum(3, 7))  # elements in [3,7]: {3,6,7} → 3

arr = [5, 2, 4, 1, 3]
print(count_inversions_fenwick(arr))  # 8

st = build_sparse_table([2,4,1,6,3,7,5])
print(rmq_query(st, 1, 4))   # min(4,1,6,3) = 1
🌲 Fenwick Tree — Why it Works

The key identity: every positive integer i has a unique binary representation, and i & (-i) isolates the lowest set bit. This means:

tree[i] stores the sum of arr[ i - LSB(i) + 1 .. i ] To get prefix_sum(i): repeatedly subtract LSB → visits O(log n) nodes To update position i: repeatedly add LSB → visits O(log n) nodes Because each position i appears in exactly O(log n) tree cells, and the cells partition the array with no overlaps in each traversal.