Counting Sort — Mathematical Foundations and Complexity Theory¶

Definition (Counting Sort). Let A[0..n-1] be an array of n integer keys
drawn from the universe U = {0, 1, ..., k-1}. Counting Sort is the
algorithm CS(A, k) -> B that produces B[0..n-1] satisfying:

  (i)  B is a permutation of A (multiset equality).
  (ii) For all 0 <= i < j < n: B[i] <= B[j]  (sorted ascending).
  (iii) Stability: for all i < j with A[i] == A[j], the element at A[i]
        appears at an index strictly less than that of A[j] in B.

The algorithm uses an auxiliary count array C[0..k-1] and an output array
B[0..n-1], runs in three straight-line passes, and performs zero
comparisons between input keys.

CS(A, k):
  C <- zero array of length k
  B <- uninitialised array of length n
  for i in 0..n-1:       C[A[i]] += 1                # Phase 1 (tally)
  for v in 1..k-1:       C[v]   += C[v-1]            # Phase 2 (prefix sum)
  for i in n-1..0 step -1:                            # Phase 3 (placement)
    v <- A[i]
    C[v] -= 1
    B[C[v]] <- A[i]
  return B

Invariant P1(i): At the start of iteration i (0 <= i <= n) of Phase 1,
                 for every v in [0, k):
                 C[v] == |{j in [0, i) : A[j] == v}|.

Base case (i=0):  C is all zeros and the set on the right is empty. P1(0) holds.

Inductive step:   Assume P1(i) for some i < n.
                  Iteration i executes C[A[i]] += 1.
                  Let w = A[i]. After the increment:
                    - For v != w: C[v] unchanged; the right-hand-side set
                      gains no member because we required A[j] == v != w.
                    - For v == w: C[w] increases by 1; the set gains element
                      A[i] = w.
                  Hence P1(i+1) holds.

Termination:      Loop runs exactly n iterations and terminates.

Postcondition:    P1(n) holds: for every v, C[v] = |{j in [0, n) : A[j] == v}|
                  = number of occurrences of v in A.

Invariant P2(v): At the start of iteration v (1 <= v <= k) of Phase 2,
                 for every u < v:
                 C[u] == |{j in [0, n) : A[j] <= u}|.

Base case (v=1):  C[0] still equals (occurrences of 0 in A) = |{j: A[j] <= 0}|.

Inductive step:   Assume P2(v).
                  Iteration v executes C[v] += C[v-1].
                  Before: C[v] = (occurrences of v in A);
                          C[v-1] = |{j: A[j] <= v-1}| by P2(v).
                  After:  C[v] = (occurrences of v) + |{j: A[j] <= v-1}|
                              = |{j: A[j] == v}| + |{j: A[j] <= v-1}|
                              = |{j: A[j] <= v}|.
                  Other entries unchanged. P2(v+1) holds.

Termination:      Loop runs k-1 iterations and terminates.

Postcondition:    P2(k) holds: for every v in [0, k):
                  C[v] = |{j: A[j] <= v}|.

Invariant P3(i): Just before iteration i (n-1 >= i >= 0) of Phase 3,
                 for every v in [0, k):
                 C[v] == T[v] - |{j in (i, n) : A[j] == v}|.

Base case (i=n-1): The set on the right is empty (j > n-1 impossible).
                   So C[v] = T[v], which is what Phase 2 left.

Inductive step:    Assume P3(i+1). Iteration (i+1) executes:
                     v <- A[i+1]
                     C[v] -= 1
                     B[C[v]] <- A[i+1]
                   Let w = A[i+1].
                   For u != w: C[u] unchanged; the set
                     {j in (i+1, n) : A[j] == u} unchanged when moving to i.
                   For u == w: C[w] decreases by 1; the set
                     {j in (i, n) : A[j] == w} = {j in (i+1, n) : A[j] == w} ∪ {i+1}
                     gains exactly one element.
                   Both sides change by 1 in the same direction. P3(i) holds.

Termination:       Loop runs n iterations (from i=n-1 down to i=0) and terminates.

Postcondition:     P3(-1) holds: for every v,
                   C[v] = T[v] - |{j in (-1, n) : A[j] == v}|
                        = T[v] - (occurrences of v)
                        = T[v-1] (or 0 if v=0).

  At every iteration where the algorithm wrote B[C[v]] <- A[i+1] with w=v=A[i+1],
  C[v] was decreased BEFORE the write. So the write happened at index
  C[v]_old - 1 in B. By P3, C[v]_old = T[v] - |{j > i+1 : A[j] == v}|,
  so the write went to position T[v] - 1 - |{j > i+1 : A[j] == v}|.

  As i decreases from n-1 to 0, the rightmost occurrence of each v lands
  at position T[v] - 1; the next at T[v] - 2; etc.
  Hence B[T[v-1] .. T[v] - 1] is filled with exactly the occurrences of v.

Claim: Counting Sort with right-to-left iteration in Phase 3 is stable.

Stability statement: For every pair (i, j) with i < j and A[i] == A[j],
                     the output position of A[i] is strictly less than
                     that of A[j].

Proof:  Let v = A[i] = A[j].
        In Phase 3, iteration j runs BEFORE iteration i (because we go
        from index n-1 down to 0 and j > i).

        Let C_j denote the value of C[v] right BEFORE iteration j executes
        its decrement. By invariant P3, after iteration j:
          C[v] becomes C_j - 1, and A[j] is written to B[C_j - 1].

        Between iteration j and iteration i, we may execute more
        decrement+write steps for value v, each at strictly smaller indices.
        Let p = number of such steps. Then right before iteration i:
          C[v] = C_j - 1 - p.

        Iteration i decrements C[v] to C_j - 2 - p and writes A[i] to
        position C_j - 2 - p.

        Compare: A[i] -> position C_j - 2 - p,
                 A[j] -> position C_j - 1.
        Since p >= 0:  C_j - 2 - p <= C_j - 2 < C_j - 1.
        So A[i]'s output position is strictly less than A[j]'s.

        Q.E.D.

Theorem (Ω(n log n) for comparison sorts):
  Any deterministic comparison-based sorting algorithm has worst-case
  running time Ω(n log n) on inputs of n distinct keys.

Proof sketch (decision tree):
  A comparison sort can be modelled as a binary decision tree whose
  internal nodes are key comparisons and whose leaves are output
  permutations. For n distinct keys there are n! possible orderings,
  so the tree must have at least n! leaves.
  A binary tree with L leaves has depth >= log2(L).
  Therefore depth >= log2(n!) = Θ(n log n).
  The depth equals the worst-case number of comparisons.
  Q.E.D.

Formally, Counting Sort lives in the Word-RAM model with word size
w >= log2(max(n, k)) bits. Its running time is O(n + k) word-RAM
operations, each of which is O(1). The comparison-model lower bound
says nothing about word-RAM algorithms that use keys as addresses.

Model:  External-Memory model (Aggarwal-Vitter)
        N = problem size in elements, M = cache size, B = block size.

Counting Sort cache-oblivious I/Os:
  Phase 1 (tally):  read A linearly -> O(N/B) reads;
                    write C at scattered positions -> O(N/B + k/B) writes
                    when C fits in cache, else O(N) when C is larger.

  Phase 2 (scan):   scan C once -> O(k/B).

  Phase 3 (place):  read A in reverse -> O(N/B);
                    write B at positions dictated by C -> O(N/B + k/B)
                    when both arrays fit in M; else O(N) random I/Os.

Total when C and B fit in M:  O(N/B + k/B).
Total when k > M:             O(N) -- one I/O per element.

External-memory comparison-sort lower bound (Aggarwal-Vitter):
  Sorting N elements requires Omega((N/B) * log_{M/B}(N/B)) I/Os.

Counting Sort I/Os when C fits in M: O(N/B + k/B).
For k <= O(N), this is O(N/B), which beats the comparison-sort lower
bound by a factor of log_{M/B}(N/B). Linear-time I/O behaviour.

Parallel Counting Sort (PRAM CREW model, p processors):

Phase 1 (parallel tally):
  Partition A into p slices of size n/p each.
  Each processor builds a local count array of size k.
  Work:  W1 = O(n + p*k)
  Span:  S1 = O(n/p + k)        -- per-processor pass + local-array init

Phase 2 (combine + prefix sum):
  Sum p local count arrays column-wise.
  Then run a work-efficient parallel prefix sum (Blelloch scan).
  Work:  W2 = O(p*k + k) = O(p*k)
  Span:  S2 = O(log p + log k)

Phase 3 (parallel scatter):
  Each processor scatters its slice to disjoint output ranges using
  per-processor offsets derived from the global prefix sum.
  Work:  W3 = O(n + p*k)        -- compute per-processor offsets + scatter
  Span:  S3 = O(n/p + k)

Totals:
  Total work:  W = O(n + p*k)
  Total span:  S = O(n/p + k + log p)
  Parallelism (W/S) = O((n + p*k) / (n/p + k + log p))

When k is constant and n >> p:
  W ~ O(n), S ~ O(n/p), parallelism ~ O(p) -- linear speedup.

When k > n/p:
  Span is dominated by the k term in Phase 1 and Phase 3 (each thread
  must process its k-sized local count array). Diminishing returns.

GPU Counting Sort (segmented scan, p = O(n) threads, k = 256):
  W = O(n)
  S = O(log n)
  Parallelism = O(n / log n)

This is span-optimal -- you cannot sort n elements in less than
Omega(log n) span on a PRAM, because the final output position of any
element depends on a count summed over the whole input.