Python/sorts/smoothsort.py

"""
Smoothsort algorithm implementation.

Smoothsort is an adaptive, in-place comparison sort invented by Edsger W. Dijkstra.
It runs in O(n log n) worst-case and degrades gracefully to O(n) for nearly sorted data.
It uses a forest of Leonardo heaps to achieve this adaptive behaviour.

Reference:
    https://en.wikipedia.org/wiki/Smoothsort
    https://www.cs.utexas.edu/~EWD/ewd07xx/EWD796a.PDF
"""

from typing import Any, Protocol


class Comparable(Protocol):
    def __lt__(self, other: Any, /) -> bool: ...


# Precomputed Leonardo numbers: L(0)=1, L(1)=1, L(k)=L(k-1)+L(k-2)+1.
# 46 values comfortably cover all practical list sizes.
_LEONARDO: list[int] = [1, 1]
while _LEONARDO[-1] < 2**31:
    _LEONARDO.append(_LEONARDO[-1] + _LEONARDO[-2] + 1)


def _sift[T: Comparable](seq: list[T], root: int, order: int) -> None:
    """
    Restore the max-heap property within a Leonardo tree of the given ``order``.

    Sifts ``seq[root]`` downward until the subtree satisfies the Leonardo
    max-heap invariant: every node is >= both of its children.
    Trees of order 0 or 1 are single nodes and already satisfy the invariant.

    In a Leonardo tree of order k rooted at index ``root``:
      - the right child root is at ``root - 1``
      - the left  child root is at ``root - 1 - L(k-2)``

    Args:
        seq:   The list being sorted (mutated in-place).
        root:  Index of the root of the Leonardo tree to fix.
        order: Leonardo order of the tree rooted at ``root``.

    Examples:
        >>> data = [3, 5, 4]
        >>> _sift(data, 2, 2)
        >>> data
        [3, 4, 5]

        >>> data = [1, 2, 3]
        >>> _sift(data, 2, 2)
        >>> data
        [1, 2, 3]

        >>> data = [7]
        >>> _sift(data, 0, 1)
        >>> data
        [7]

        >>> data = [9, 1, 8, 5, 3]
        >>> _sift(data, 4, 3)
        >>> data
        [3, 1, 9, 5, 8]
    """
    while order > 1:
        right = root - 1  # right child root
        left = root - 1 - _LEONARDO[order - 2]  # left child root

        if not (seq[left] < seq[right]) and seq[root] < seq[left]:
            seq[root], seq[left] = seq[left], seq[root]
            root = left
            order -= 1
        elif seq[left] < seq[right] and seq[root] < seq[right]:
            seq[root], seq[right] = seq[right], seq[root]
            root = right
            order -= 2
        else:
            break


def _trinkle[T: Comparable](
    seq: list[T],
    pos: int,
    heap_sizes: list[int],
    idx: int,
) -> None:
    """
    Restore both the inter-heap root ordering and the intra-heap ordering.

    Walks the value at ``pos`` leftwards through the forest-root chain as
    long as the left-neighbour root is larger, then calls ``_sift`` to fix
    the heap at the final resting position.

    Args:
        seq:        The list being sorted (mutated in-place).
        pos:        Index of the root being inserted or newly exposed.
        heap_sizes: List of Leonardo orders for the current forest (left to
                    right); ``heap_sizes[idx]`` is the order of the tree
                    whose root is at ``pos``.
        idx:        Position in ``heap_sizes`` for the tree rooted at ``pos``.

    Examples:
        >>> data = [1, 5, 3]
        >>> _trinkle(data, 2, [1, 1], 1)
        >>> data
        [1, 3, 5]

        >>> data = [3, 5, 4]
        >>> _trinkle(data, 2, [2], 0)
        >>> data
        [3, 4, 5]
    """
    while idx > 0:
        prev_root = pos - _LEONARDO[heap_sizes[idx]]
        if not (seq[pos] < seq[prev_root]):
            break
        # Only swap if prev_root is also > its own children; otherwise
        # moving it would break the heap on the left side.
        if heap_sizes[idx] > 1:
            right = pos - 1
            left = pos - 1 - _LEONARDO[heap_sizes[idx] - 2]
            if not (seq[right] < seq[prev_root]) or not (seq[left] < seq[prev_root]):
                break
        seq[pos], seq[prev_root] = seq[prev_root], seq[pos]
        pos = prev_root
        idx -= 1

    _sift(seq, pos, heap_sizes[idx])


def smoothsort[T: Comparable](seq: list[T]) -> list[T]:
    """
    Sort a list in-place using the Smoothsort algorithm and return it.

    Smoothsort (Edsger W. Dijkstra, 1981) is an adaptive, in-place sort
    with O(n log n) worst-case time and O(n) best-case time on already-sorted
    input.  It improves on Heapsort by maintaining a forest of Leonardo heaps
    whose structure mirrors the sorted prefix of the sequence.

    Args:
        seq: A list of mutually comparable items to sort.

    Returns:
        The same list object, sorted in ascending order.

    Examples:
        >>> smoothsort([4, 1, 3, 9, 7])
        [1, 3, 4, 7, 9]
        >>> smoothsort([])
        []
        >>> smoothsort([1])
        [1]
        >>> smoothsort([5, 4, 3, 2, 1])
        [1, 2, 3, 4, 5]
        >>> smoothsort([3, 3, 2, 1, 2])
        [1, 2, 2, 3, 3]
        >>> smoothsort(["d", "a", "c", "b"])
        ['a', 'b', 'c', 'd']
        >>> smoothsort([2.5, -1, 0.0])
        [-1, 0.0, 2.5]
        >>> smoothsort([1, 2, 3, 4, 5])
        [1, 2, 3, 4, 5]
        >>> smoothsort([-3, 0, -1, 5, 2])
        [-3, -1, 0, 2, 5]
        >>> smoothsort([1, "a"])
        Traceback (most recent call last):
            ...
        TypeError: '<' not supported between instances of 'str' and 'int'
    """
    n = len(seq)
    if n < 2:
        return seq

    # ``heap_sizes[i]`` is the Leonardo order of the i-th tree (left to right).
    heap_sizes: list[int] = []

    # ------------------------------------------------------------------
    # Phase 1 : Build the Leonardo heap forest over seq[0..n-1].
    # ------------------------------------------------------------------
    for i in range(n):
        # If the two rightmost trees have consecutive orders, merge them.
        if len(heap_sizes) >= 2 and heap_sizes[-2] == heap_sizes[-1] + 1:
            heap_sizes.pop()
            heap_sizes[-1] += 1
        elif heap_sizes and heap_sizes[-1] == 1:
            heap_sizes.append(0)
        else:
            heap_sizes.append(1)

        _trinkle(seq, i, heap_sizes, len(heap_sizes) - 1)

    # ------------------------------------------------------------------
    # Phase 2 : Extract maximum elements right-to-left.
    # ------------------------------------------------------------------
    for i in range(n - 1, -1, -1):
        order = heap_sizes.pop()
        if order > 1:
            # Expose the two child roots and re-trinkle each.
            right_order = order - 2
            left_order = order - 1
            right_pos = i - 1
            left_pos = i - 1 - _LEONARDO[right_order]

            heap_sizes.append(left_order)
            _trinkle(seq, left_pos, heap_sizes, len(heap_sizes) - 1)

            heap_sizes.append(right_order)
            _trinkle(seq, right_pos, heap_sizes, len(heap_sizes) - 1)

    return seq


if __name__ == "__main__":
    import doctest
    import random

    results = doctest.testmod(verbose=False)
    assert results.failed == 0, f"{results.failed} doctest(s) failed"

    for trial in range(5000):
        sample = random.choices(range(-50, 50), k=random.randint(0, 30))
        got = smoothsort(sample[:])
        assert got == sorted(sample), (
            f"Trial {trial}: smoothsort({sample!r}) -> {got!r}, "
            f"expected {sorted(sample)!r}"
        )

    print("All doctests and 5 000 random trials passed.")