Insertion, Selection & Shellsort: Inversions & Diminishing Increments

The Combinatorics of Disorder: Inversions

Given an array $A[0 \dots n-1]$, an inversion is any ordered pair of indices $(i, j)$ such that $i < j$ and $A[i] > A[j]$. The total inversion count $I(A)$ quantifies exactly how far $A$ is from sorted order:

• A sorted array has $I(A) = 0$ inversions.
• A strictly reverse-sorted array has every pair inverted: $I(A) = \binom{n}{2} = \frac{n(n-1)}{2}$ inversions.
• Expected Inversions of a Random Permutation: By linearity of expectation over all $\binom{n}{2}$ pairs, each pair $(i, j)$ has probability exactly $1/2$ of satisfying $A[i] > A[j]$ in a uniformly random permutation of distinct keys. Thus: $$\mathbb{E}[I(A)] = \frac{1}{2}\binom{n}{2} = \frac{n(n-1)}{4} = \Theta(n^2)$$

Shellsort: Breaking the $\Omega(n^2)$ Barrier via $h$-Sorting

Donald Shell (1959) observed that moving an element across $h$ positions in a single swap can eliminate up to $2h - 1$ inversions simultaneously. An array is called $h$-sorted if every strided subsequence $A[i], A[i+h], A[i+2h], \dots$ is independently sorted via Insertion Sort with stride $h$. Shellsort executes $h$-Insertion Sort for a diminishing increment sequence $h_t > h_{t-1} > \dots > h_1 = 1$.

• The $h$-Sort Preservation Theorem: If an array is already $g$-sorted and we subsequently $h$-sort it, the array remains $g$-sorted! This nontrivial combinatorial lemma (proven via Gale-Ryser / Frobenius coin problem duality by Pratt and Knuth) guarantees that work done at coarse strides is never undone by finer strides.
• Why Powers of 2 Fail ($h = 2^k$): Shell's original sequence $N/2, N/4, \dots, 2, 1$ never compares elements at odd indices with elements at even indices until the very final $h = 1$ pass. Interleaving the largest $n/2$ elements at odd indices and the smallest $n/2$ elements at even indices forces the $h = 1$ pass to do $\Omega(n^2)$ work.
• Coprime Diminishing Increments:
  • Hibbard ($2^k - 1$): Adjacent increments are coprime; worst-case $\Theta(n^{3/2})$ (Papernov & Stasevich, 1965).
  • Knuth ($\frac{3^k - 1}{2} = 1, 4, 13, 40, 121, \dots$): Easy recurrence $h = 3h + 1$, worst-case $\Theta(n^{3/2})$.
  • Sedgewick ($4^k + 3 \cdot 2^{k-1} + 1$) / Tokuda / Ciura ($1, 4, 10, 23, 57, 132, 301, 701, 1750$): Achieves $O(n^{4/3})$ worst-case and matches $O(n \log n)$ algorithms in wall-clock time up to $n \approx 10^4$.

Production C++20 Implementation: Unguarded Insertion Sort & Ciura Shellsort

#include <algorithm>
#include <array>
#include <concepts>
#include <cstddef>
#include <utility>
#include <vector>

// Guarded + Unguarded Insertion Sort (used as base-case leaf in Introsort / pdqsort)
template <std::totally_ordered T>
void insertion_sort_unguarded(T* first, T* last) noexcept {
    if (last - first <= 1) return;

    // 1. Move the minimum of the first step, or guard only when touching A[0]
    for (T* i = first + 1; i != last; ++i) {
        T val = std::move(*i);
        // Fast path: if val < *first, it belongs at the very front (shift whole prefix)
        if (val < *first) {
            std::move_backward(first, i, i + 1);
            *first = std::move(val);
        } else {
            // Unguarded inner loop: *first <= val acts as a natural sentinel!
            // Eliminates the j > first bounds check branch on every iteration.
            T* j = i;
            while (val < *(j - 1)) {
                *j = std::move(*(j - 1));
                --j;
            }
            *j = std::move(val);
        }
    }
}

// Shellsort using Marcin Ciura's empirically optimal increment sequence (extended by *2.25)
template <std::totally_ordered T>
void shellsort_ciura(std::vector<T>& a) noexcept {
    const std::size_t n = a.size();
    if (n <= 1) return;

    // Ciura's sequence + geometric 2.25x extension for large N
    std::vector<std::size_t> gaps = {1, 4, 10, 23, 57, 132, 301, 701, 1750};
    while (gaps.back() < n / 2) {
        gaps.push_back(static_cast<std::size_t>(gaps.back() * 2.25));
    }

    for (auto it = gaps.rbegin(); it != gaps.rend(); ++it) {
        const std::size_t h = *it;
        for (std::size_t i = h; i < n; ++i) {
            T temp = std::move(a[i]);
            std::size_t j = i;
            while (j >= h && temp < a[j - h]) {
                a[j] = std::move(a[j - h]);
                j -= h;
            }
            a[j] = std::move(temp);
        }
    }
}