Sorting Algorithms

Before reaching for a clever algorithm, it pays to understand the naive ones deeply — they reveal the cost structure of sorting and give us the loop-invariant reasoning technique.

Insertion Sort works exactly like sorting a hand of playing cards. You maintain a sorted prefix of the array and, one by one, take the next unsorted element and insert it into its correct position by shifting larger elements rightward.

• Loop invariant: after iteration i, A[0..i] is a sorted permutation of the original A[0..i].
• Time: O(n²) worst/average (reverse-sorted input), O(n) best (already sorted). The constant is small, making it practical for small n (typically n ≤ 10–20).
• Space: O(1) in-place.
• Stable: yes — the inner loop stops when it finds an equal element, preserving relative order.

Selection Sort finds the minimum of the unsorted suffix and swaps it to the front, repeatedly.

• Time: always exactly O(n²) — it does n(n-1)/2 comparisons regardless of input. This makes it worse than insertion sort in practice.
• Stable: in the naive swap-based version, no — the swap can move a later element before an equal earlier one.

Bubble Sort repeatedly sweeps the array, swapping adjacent out-of-order pairs. Each sweep bubbles the largest unsorted element to its final position.

• Time: O(n²) worst/average, O(n) best (with early-exit optimization).
• Stable: yes.
• The key insight: bubble sort performs at most n(n-1)/2 swaps, but also at most n(n-1)/2 comparisons. It is rarely used in practice because it has higher constant factors than insertion sort.

The shared bottleneck of all three: each inversion (a pair (i,j) where i < j but A[i] > A[j]) requires at least one swap or shift to fix, and a reverse-sorted array of n elements has n(n-1)/2 inversions. Any algorithm that fixes one inversion per operation is O(n²).

def insertion_sort(A):
    for i in range(1, len(A)):
        key = A[i]
        j = i - 1
        # shift elements of A[0..i-1] that are > key one position ahead
        while j >= 0 and A[j] > key:
            A[j + 1] = A[j]
            j -= 1
        A[j + 1] = key
    return A

# Trace on [5, 2, 4, 6, 1, 3]:
# i=1: key=2, shift 5 → [2,5,4,6,1,3]
# i=2: key=4, shift 5 → [2,4,5,6,1,3]
# i=3: key=6, no shift  → [2,4,5,6,1,3]
# i=4: key=1, shift 2,4,5,6 → [1,2,4,5,6,3]
# i=5: key=3, shift 4,5 → [1,2,3,4,5,6]

Merge sort is the textbook divide-and-conquer algorithm. The key insight is that merging two sorted lists is easy and linear, so if we can recursively sort two halves, we can combine them in O(n) time — yielding a total cost of O(n log n).

Algorithm:

1. Divide: split the array into two halves of size n/2.
2. Conquer: recursively sort each half.
3. Combine: merge the two sorted halves into a single sorted array by repeatedly taking the smaller front element.

Recurrence: T(n) = 2T(n/2) + O(n). By the Master Theorem (case 2), T(n) = O(n log n).

You can also see this geometrically: the recursion has log₂ n levels, and at each level the total work across all merge operations is exactly O(n). So the total is n × log n.

Properties:

• Stable: yes, provided the merge step breaks ties by taking from the left half first.
• Not in-place: requires O(n) auxiliary space for the merge step. This is a real cost — it doubles memory usage, which matters for cache efficiency.
• External sorting: merge sort is the algorithm of choice when data does not fit in RAM, because its access pattern is sequential (cache-friendly in an external context).

Timsort (Python's built-in sort(), Java's Arrays.sort for objects) is a production-hardened hybrid that finds existing sorted runs in the input and merges them intelligently. It achieves O(n) on nearly-sorted inputs and O(n log n) worst case, while remaining stable.

def merge_sort(A):
    if len(A) <= 1:
        return A
    mid = len(A) // 2
    left  = merge_sort(A[:mid])
    right = merge_sort(A[mid:])
    return merge(left, right)

def merge(L, R):
    result = []
    i = j = 0
    while i < len(L) and j < len(R):
        # '<=' (not '<') preserves stability
        if L[i] <= R[j]:
            result.append(L[i]); i += 1
        else:
            result.append(R[j]); j += 1
    result.extend(L[i:])
    result.extend(R[j:])
    return result

Quicksort is, in practice, the fastest comparison sort for in-memory data on modern hardware. Invented by Tony Hoare in 1959, it is a beautiful example of probabilistic algorithm analysis.

Algorithm:

1. Choose a pivot element p.
2. Partition: rearrange the array so all elements ≤ p come before all elements > p (with p in its final position).
3. Recurse on the left and right sub-arrays.

Worst case: If the pivot is always the minimum or maximum (e.g., always picking the first element of an already-sorted array), one partition has n−1 elements and the other has 0. Recurrence: T(n) = T(n−1) + O(n) → O(n²).

Fix — randomized quicksort: choose the pivot uniformly at random. Now the analysis becomes a beautiful probabilistic argument:

• The probability that elements i and j (in sorted order) are compared is exactly 2/(j−i+1).
• The expected total comparisons = Σᵢ Σⱼ>ᵢ 2/(j−i+1) = 2n ln n + O(n) ≈ 1.386 n log₂ n.

So randomized quicksort runs in expected O(n log n) for any input. The probability of significantly exceeding this (say, > cn log n for large c) is extremely small.

Why faster than merge sort in practice?

• In-place: no auxiliary array required, much better cache behavior.
• Smaller constants: the partition step is a simple linear scan with swaps.
• Locality: once a sub-array fits in cache, it stays there.

Partition schemes:

• Lomuto partition: simpler to implement, slightly more swaps.
• Hoare partition: original, slightly fewer swaps, but trickier to implement correctly (especially the pivot-final-position logic).

Three-way quicksort (Bentley-McIlroy): when there are many equal elements, a standard partition puts them all in one side and degrades. Three-way partitioning (Dutch National Flag) splits into < p, = p, > p in O(n), achieving O(n) on inputs with only O(1) distinct keys.

Stability: quicksort as typically implemented is not stable.

import random

def quicksort(A, lo, hi):
    if lo < hi:
        p = partition(A, lo, hi)
        quicksort(A, lo, p - 1)
        quicksort(A, p + 1, hi)

def partition(A, lo, hi):
    # randomize pivot to avoid worst case
    r = random.randint(lo, hi)
    A[r], A[hi] = A[hi], A[r]
    pivot = A[hi]
    i = lo - 1
    for j in range(lo, hi):
        if A[j] <= pivot:
            i += 1
            A[i], A[j] = A[j], A[i]
    A[i + 1], A[hi] = A[hi], A[i + 1]
    return i + 1

# Usage: quicksort(A, 0, len(A)-1)

Heapsort achieves the best of both worlds that the other comparison sorts could not simultaneously deliver: O(n log n) worst-case (unlike quicksort) and O(1) extra space (unlike merge sort). The price is a larger constant factor and poor cache performance.

The heap data structure: a binary heap is a complete binary tree stored in an array where for a max-heap, every node is ≥ its children. In a zero-indexed array, the children of node i are at 2i+1 and 2i+2; the parent is at (i−1)//2.

Algorithm (two phases):

1. Build-max-heap: call max_heapify on every internal node from bottom to top. This takes O(n) time — not O(n log n) as you might expect! The key insight is that nodes near the bottom have small subtrees and sift down quickly.
2. Sort: repeatedly extract the maximum (swap root with last element, reduce heap size by 1, sift root down). Each extraction is O(log n), and we do n of them: O(n log n) total.

Build-heap O(n) proof: The work at height h is O(h), and there are at most ⌈n/2^(h+1)⌉ nodes at height h. Total work = Σ_{h=0}^{⌊log n⌋} ⌈n/2^(h+1)⌉ × O(h) = O(n Σ h/2^h) = O(n × 2) = O(n). (The sum Σ h/2^h converges to 2.)

Cache behavior: heapsort is notoriously cache-unfriendly. The sift-down operation jumps between widely separated array positions (root to leaf), causing frequent cache misses. This is why quicksort is faster in practice despite heapsort's better worst-case guarantee.

When to use: heapsort is the right choice when you need a guaranteed O(n log n) bound and cannot afford O(n) extra space (e.g., embedded systems). It also forms the basis of priority queues, which are used pervasively in graph algorithms like Dijkstra's.

def heapsort(A):
    n = len(A)

    # Phase 1: build max-heap in O(n)
    for i in range(n // 2 - 1, -1, -1):
        sift_down(A, i, n)

    # Phase 2: extract max repeatedly
    for end in range(n - 1, 0, -1):
        A[0], A[end] = A[end], A[0]   # move current max to sorted suffix
        sift_down(A, 0, end)

def sift_down(A, root, end):
    """Restore max-heap property for subtree rooted at 'root', heap size 'end'."""
    while True:
        left  = 2 * root + 1
        right = 2 * root + 2
        largest = root
        if left  < end and A[left]  > A[largest]: largest = left
        if right < end and A[right] > A[largest]: largest = right
        if largest == root:
            break
        A[root], A[largest] = A[largest], A[root]
        root = largest

Here is one of the most beautiful impossibility proofs in all of computer science: any algorithm that sorts n elements purely by comparing pairs of elements must make at least ⌊log₂(n!)⌋ = Ω(n log n) comparisons in the worst case. No matter how clever you are, you cannot beat this.

The decision tree model: model any deterministic comparison sort as a binary tree. Each internal node represents a comparison A[i] : A[j]; its two children represent the two outcomes (≤ or >). Each leaf represents a complete sorted order — one of the n! possible permutations of the input. Any correct sort must reach a distinct leaf for each permutation (otherwise two different inputs would be sorted identically, which is wrong).

Lower bound argument:

1. The decision tree has at least n! leaves.
2. A binary tree of height h has at most 2^h leaves.
3. Therefore: 2^h ≥ n! → h ≥ log₂(n!).
4. By Stirling's approximation: log₂(n!) = n log₂ n − n log₂ e + O(log n) = Ω(n log n).

The height of the decision tree equals the worst-case number of comparisons, so every comparison-based sort requires Ω(n log n) comparisons.

Information-theoretic intuition: there are n! possible inputs (permutations). Each comparison gives you 1 bit of information. To distinguish among n! possibilities, you need at least log₂(n!) bits, i.e., log₂(n!) comparisons. Merge sort and heapsort achieve this bound (up to constants), making them asymptotically optimal comparison sorts.

What this does NOT rule out: sorting algorithms that do not restrict themselves to pairwise comparisons — they can read individual bits of the keys, use arithmetic on keys, or exploit the fact that keys come from a bounded universe. This is the door through which radix sort and counting sort walk.

For randomized algorithms: the lower bound also applies in expectation — even a randomized comparison sort must use Ω(n log n) comparisons in expectation on some input.

# Concrete illustration of the decision tree bound:
import math

def lower_bound_comparisons(n):
    """Minimum comparisons any comparison sort must make."""
    log_n_factorial = sum(math.log2(i) for i in range(1, n + 1))
    return math.ceil(log_n_factorial)

for n in [4, 8, 16, 32, 64]:
    lb = lower_bound_comparisons(n)
    nlogn = n * math.log2(n)
    print(f"n={n:3d}: lower bound = {lb:5d}, n*log2(n) = {nlogn:7.1f}")

# n=  4: lower bound =     5, n*log2(n) =     8.0
# n=  8: lower bound =    16, n*log2(n) =    24.0
# n= 16: lower bound =    44, n*log2(n) =    64.0
# n= 32: lower bound =   118, n*log2(n) =   160.0
# n= 64: lower bound =   296, n*log2(n) =   384.0

Counting sort breaks the Ω(n log n) barrier by exploiting knowledge about the range of keys, not just their relative order. If all keys are integers in [0, k), we can sort in O(n + k) time.

Algorithm:

1. Allocate a count array C[0..k-1] initialized to zero.
2. For each element in the input, increment C[A[i]].
3. Compute prefix sums of C so that C[j] now stores the number of elements ≤ j (i.e., the correct final position of the last element equal to j).
4. Iterate the input in reverse (for stability), placing each element at position C[A[i]] - 1 in the output, then decrement C[A[i]].

Time: O(n + k). Space: O(n + k).

Why reverse iteration for stability? By processing the input right-to-left in step 4, elements with equal keys appear in the output in the same left-to-right order as they appeared in the input — i.e., the sort is stable. This stability is what makes counting sort a useful subroutine inside radix sort.

When is counting sort practical? When k = O(n), it is truly O(n). But if k = 10⁹ (e.g., sorting 32-bit integers), allocating an array of size 10⁹ is infeasible. The solution: radix sort.

Intuition: counting sort does not compare elements at all — it simply counts occurrences. It only works because integer keys can be used as array indices. This is not a comparison sort in any sense, so the Ω(n log n) lower bound does not apply.

def counting_sort(A, k):
    """Sort A where all elements are integers in [0, k). Stable."""
    C = [0] * k
    for x in A:
        C[x] += 1
    # prefix sums: C[j] = number of elements <= j
    for j in range(1, k):
        C[j] += C[j - 1]
    B = [0] * len(A)
    # iterate in reverse for stability
    for i in range(len(A) - 1, -1, -1):
        B[C[A[i]] - 1] = A[i]
        C[A[i]] -= 1
    return B

# Example: sort [3,1,4,1,5,9,2,6] with k=10
# Returns: [1,1,2,3,4,5,6,9]

Radix sort extends counting sort to handle large keys by processing them one digit at a time. The key insight, first used in punch-card sorting machines in the 1890s, is: sort by the least significant digit first, then the next digit, and so on. Stability ensures that each pass correctly refines the order established by previous passes.

LSD (Least Significant Digit) Radix Sort:

1. For d passes (d = number of digits in base b):
   • Apply a stable sort (counting sort) to the current digit position.
2. After all d passes, the array is sorted.

Why LSD and not MSD? If you sort by the most significant digit first, you'd need to recurse into each bucket independently. LSD avoids this by using stability: after sorting by digit i, elements with the same digit i are in the order produced by sorting digit i−1. This invariant propagates correctly.

Time analysis: each pass of counting sort on a single digit takes O(n + b) where b is the base (alphabet size for that digit). With d passes: O(d(n + b)).

For n 32-bit integers, choose b = 256 (8-bit digits), so d = 4 passes. Each pass is O(n + 256) = O(n). Total: O(4n) = O(n). In practice this is extremely fast.

Space: O(n + b) for the output and count arrays.

MSD Radix Sort: sort by most significant digit first, recursing into each bucket. More complex but useful for variable-length strings and when early termination is possible (most strings are distinguished by early characters).

The real world: radix sort is used in GPU sort implementations, network packet classification, and integer sorting libraries where keys are bounded integers.

def radix_sort(A):
    """LSD radix sort for non-negative integers."""
    if not A:
        return A
    max_val = max(A)
    exp = 1  # 1, 10, 100, ...
    while max_val // exp > 0:
        A = counting_sort_by_digit(A, exp)
        exp *= 10
    return A

def counting_sort_by_digit(A, exp):
    """Stable sort by the digit at position exp (base 10)."""
    n = len(A)
    output = [0] * n
    count  = [0] * 10
    for x in A:
        count[(x // exp) % 10] += 1
    for i in range(1, 10):
        count[i] += count[i - 1]
    for i in range(n - 1, -1, -1):  # reverse for stability
        d = (A[i] // exp) % 10
        output[count[d] - 1] = A[i]
        count[d] -= 1
    return output

# radix_sort([170, 45, 75, 90, 802, 24, 2, 66])
# → [2, 24, 45, 66, 75, 90, 170, 802]

Bucket sort takes a different approach to beating Ω(n log n): it assumes the input is drawn from a known, roughly uniform distribution over [0, 1). By distributing elements into n equally-spaced buckets and sorting each bucket independently, we achieve expected O(n) time.

Algorithm:

1. Create n empty buckets corresponding to [0, 1/n), [1/n, 2/n), ..., [(n-1)/n, 1).
2. For each element x ∈ [0,1), place it in bucket ⌊n × x⌋.
3. Sort each bucket (insertion sort is fine — each bucket has O(1) elements in expectation).
4. Concatenate the buckets.

Time analysis (expected O(n) for uniform input):

• The expected number of elements in each bucket is 1.
• The expected cost to sort a bucket of k elements with insertion sort is O(k²) in the worst case, but the expected cost over uniformly distributed elements is O(1) per bucket.
• Total expected time: O(n).

More carefully: let nᵢ be the number of elements in bucket i. Sort cost = Σᵢ O(nᵢ²). The expected value E[nᵢ²] = (n−1)/n × (1/n) + (1/n)² × (using second moment of binomial) = 2 − 1/n. So E[Σ nᵢ²] = O(n).

When does it degrade? If the input is not uniformly distributed, many elements may land in the same bucket, and the sort degrades to O(n²). Bucket sort is the right tool when you have a good model of the input distribution.

Generalizations: The idea extends to any distribution — if you can compute a rank function that approximately evenly distributes keys, you can bucket sort effectively. This is the basis of some parallel sorting algorithms.

def bucket_sort(A):
    """Sort floats in [0, 1). Expected O(n) for uniform input."""
    n = len(A)
    if n == 0:
        return A
    buckets = [[] for _ in range(n)]
    for x in A:
        idx = int(n * x)          # bucket index in [0, n-1]
        if idx == n: idx = n - 1  # handle x == 1.0 edge case
        buckets[idx].append(x)
    result = []
    for bucket in buckets:
        bucket.sort()             # insertion sort is fine for small buckets
        result.extend(bucket)
    return result

import random
data = [random.random() for _ in range(10000)]
sorted_data = bucket_sort(data)
assert sorted_data == sorted(data)

A sort is stable if elements with equal keys appear in the output in the same relative order as they appeared in the input. This property may seem like a minor nicety, but it is crucial in practice.

Why stability matters — the multi-key sort: Suppose you have a table of employees, and you want to sort first by department, then by salary within each department. The way to do this with a single-key sort is:

1. First, sort by salary (the secondary key).
2. Then, sort by department (the primary key) using a stable sort.

Because the second sort is stable, employees with the same department retain the salary order from step 1. If the second sort were unstable, the within-department salary order would be scrambled.

This is exactly how database ORDER BY with multiple columns works, and how spreadsheet sorts work.

Which algorithms are stable?

Making an unstable sort stable: any sort can be made stable by augmenting each key with its original index as a tiebreaker. This adds O(n) space but guarantees stability. Python's sorted() uses this trick internally (via the Timsort algorithm, which is natively stable).

Real-world example: Python's list.sort() and sorted() are guaranteed stable. Java's Arrays.sort for objects uses merge sort (stable), but Arrays.sort for primitives uses a dual-pivot quicksort (not stable — but stability is meaningless for primitives since equal primitives are indistinguishable).

# Stability demonstration: sort records by department, then by salary
employees = [
    {"name": "Alice",   "dept": "Eng",  "salary": 90},
    {"name": "Bob",     "dept": "Mkt",  "salary": 70},
    {"name": "Charlie", "dept": "Eng",  "salary": 80},
    {"name": "Diana",   "dept": "Mkt",  "salary": 85},
]

# Step 1: stable sort by salary
step1 = sorted(employees, key=lambda e: e["salary"])
# [Bob/70, Charlie/80, Diana/85, Alice/90]

# Step 2: stable sort by department
step2 = sorted(step1, key=lambda e: e["dept"])
# Eng: [Charlie/80, Alice/90]  <- salary order preserved
# Mkt: [Bob/70, Diana/85]      <- salary order preserved

for e in step2:
    print(e["dept"], e["name"], e["salary"])
# Eng Charlie 80
# Eng Alice   90
# Mkt Bob     70
# Mkt Diana   85

No single sorting algorithm is best in all circumstances. The right choice depends on the size of the input, the distribution of keys, memory constraints, whether stability is needed, and the hardware.

Summary table:

Introsort: the algorithm used in C++ std::sort. It starts with quicksort for average-case speed, switches to heapsort when recursion depth exceeds 2 log n (to avoid quicksort's O(n²) worst case), and switches to insertion sort for sub-arrays of size ≤ 16. It combines the strengths of all three comparison sorts.

Timsort: Python/Java's stable sort. Finds natural sorted runs in the input, extends short runs with insertion sort, then merges runs using a sophisticated merge strategy that exploits real-world structure (nearly-sorted data, repeated elements).

Key decision rules:

• n < ~20: insertion sort
• Keys are bounded integers in a small range: counting sort or radix sort
• Keys are floating-point from a known distribution: bucket sort
• Need stability + external: merge sort
• General purpose, performance critical: introsort (quicksort + heapsort + insertion sort hybrid)
• Need stable + general: Timsort