|
| 1 | +""" |
| 2 | +Smoothsort algorithm implementation. |
| 3 | +
|
| 4 | +Smoothsort is an adaptive, in-place comparison sort invented by Edsger W. Dijkstra. |
| 5 | +It runs in O(n log n) worst-case and degrades gracefully to O(n) for nearly sorted data. |
| 6 | +It uses a forest of Leonardo heaps to achieve this adaptive behaviour. |
| 7 | +
|
| 8 | +Reference: |
| 9 | + https://en.wikipedia.org/wiki/Smoothsort |
| 10 | + https://www.cs.utexas.edu/~EWD/ewd07xx/EWD796a.PDF |
| 11 | +""" |
| 12 | + |
| 13 | +# Precomputed Leonardo numbers: L(0)=1, L(1)=1, L(k)=L(k-1)+L(k-2)+1. |
| 14 | +# 46 values comfortably cover all practical list sizes. |
| 15 | +_LEONARDO: list[int] = [1, 1] |
| 16 | +while _LEONARDO[-1] < 2**31: |
| 17 | + _LEONARDO.append(_LEONARDO[-1] + _LEONARDO[-2] + 1) |
| 18 | + |
| 19 | + |
| 20 | +def _sift(seq: list[int], root: int, order: int) -> None: |
| 21 | + """ |
| 22 | + Restore the max-heap property within a Leonardo tree of the given ``order``. |
| 23 | +
|
| 24 | + Sifts ``seq[root]`` downward until the subtree satisfies the Leonardo |
| 25 | + max-heap invariant: every node is >= both of its children. |
| 26 | + Trees of order 0 or 1 are single nodes and already satisfy the invariant. |
| 27 | +
|
| 28 | + In a Leonardo tree of order k rooted at index ``root``: |
| 29 | + - the right child root is at ``root - 1`` |
| 30 | + - the left child root is at ``root - 1 - L(k-2)`` |
| 31 | +
|
| 32 | + Args: |
| 33 | + seq: The list being sorted (mutated in-place). |
| 34 | + root: Index of the root of the Leonardo tree to fix. |
| 35 | + order: Leonardo order of the tree rooted at ``root``. |
| 36 | +
|
| 37 | + Examples: |
| 38 | + >>> data = [3, 5, 4] |
| 39 | + >>> _sift(data, 2, 2) |
| 40 | + >>> data |
| 41 | + [3, 4, 5] |
| 42 | +
|
| 43 | + >>> data = [1, 2, 3] |
| 44 | + >>> _sift(data, 2, 2) |
| 45 | + >>> data |
| 46 | + [1, 2, 3] |
| 47 | +
|
| 48 | + >>> data = [7] |
| 49 | + >>> _sift(data, 0, 1) |
| 50 | + >>> data |
| 51 | + [7] |
| 52 | +
|
| 53 | + >>> data = [9, 1, 8, 5, 3] |
| 54 | + >>> _sift(data, 4, 3) |
| 55 | + >>> data |
| 56 | + [3, 1, 9, 5, 8] |
| 57 | + """ |
| 58 | + while order > 1: |
| 59 | + right = root - 1 # right child root |
| 60 | + left = root - 1 - _LEONARDO[order - 2] # left child root |
| 61 | + |
| 62 | + if seq[left] >= seq[right] and seq[left] > seq[root]: |
| 63 | + seq[root], seq[left] = seq[left], seq[root] |
| 64 | + root = left |
| 65 | + order -= 1 |
| 66 | + elif seq[right] > seq[left] and seq[right] > seq[root]: |
| 67 | + seq[root], seq[right] = seq[right], seq[root] |
| 68 | + root = right |
| 69 | + order -= 2 |
| 70 | + else: |
| 71 | + break |
| 72 | + |
| 73 | + |
| 74 | +def _trinkle( |
| 75 | + seq: list[int], |
| 76 | + pos: int, |
| 77 | + heap_sizes: list[int], |
| 78 | + idx: int, |
| 79 | +) -> None: |
| 80 | + """ |
| 81 | + Restore both the inter-heap root ordering and the intra-heap ordering. |
| 82 | +
|
| 83 | + Walks the value at ``pos`` leftwards through the forest-root chain as |
| 84 | + long as the left-neighbour root is larger, then calls ``_sift`` to fix |
| 85 | + the heap at the final resting position. |
| 86 | +
|
| 87 | + Args: |
| 88 | + seq: The list being sorted (mutated in-place). |
| 89 | + pos: Index of the root being inserted or newly exposed. |
| 90 | + heap_sizes: List of Leonardo orders for the current forest (left to |
| 91 | + right); ``heap_sizes[idx]`` is the order of the tree |
| 92 | + whose root is at ``pos``. |
| 93 | + idx: Position in ``heap_sizes`` for the tree rooted at ``pos``. |
| 94 | +
|
| 95 | + Examples: |
| 96 | + >>> data = [1, 5, 3] |
| 97 | + >>> _trinkle(data, 2, [1, 1], 1) |
| 98 | + >>> data |
| 99 | + [1, 3, 5] |
| 100 | +
|
| 101 | + >>> data = [3, 5, 4] |
| 102 | + >>> _trinkle(data, 2, [2], 0) |
| 103 | + >>> data |
| 104 | + [3, 4, 5] |
| 105 | + """ |
| 106 | + while idx > 0: |
| 107 | + prev_root = pos - _LEONARDO[heap_sizes[idx]] |
| 108 | + if seq[pos] >= seq[prev_root]: |
| 109 | + break |
| 110 | + # Only swap if prev_root is also >= its own children; otherwise |
| 111 | + # moving it would break the heap on the left side. |
| 112 | + if heap_sizes[idx] > 1: |
| 113 | + right = pos - 1 |
| 114 | + left = pos - 1 - _LEONARDO[heap_sizes[idx] - 2] |
| 115 | + if seq[prev_root] <= seq[right] or seq[prev_root] <= seq[left]: |
| 116 | + break |
| 117 | + seq[pos], seq[prev_root] = seq[prev_root], seq[pos] |
| 118 | + pos = prev_root |
| 119 | + idx -= 1 |
| 120 | + |
| 121 | + _sift(seq, pos, heap_sizes[idx]) |
| 122 | + |
| 123 | + |
| 124 | +def smoothsort(seq: list[int]) -> list[int]: |
| 125 | + """ |
| 126 | + Sort a list in-place using the Smoothsort algorithm and return it. |
| 127 | +
|
| 128 | + Smoothsort (Edsger W. Dijkstra, 1981) is an adaptive, in-place sort |
| 129 | + with O(n log n) worst-case time and O(n) best-case time on already-sorted |
| 130 | + input. It improves on Heapsort by maintaining a forest of Leonardo heaps |
| 131 | + whose structure mirrors the sorted prefix of the sequence. |
| 132 | +
|
| 133 | + Args: |
| 134 | + seq: A list of integers to sort. |
| 135 | +
|
| 136 | + Returns: |
| 137 | + The same list object, sorted in ascending order. |
| 138 | +
|
| 139 | + Examples: |
| 140 | + >>> smoothsort([4, 1, 3, 9, 7]) |
| 141 | + [1, 3, 4, 7, 9] |
| 142 | + >>> smoothsort([]) |
| 143 | + [] |
| 144 | + >>> smoothsort([1]) |
| 145 | + [1] |
| 146 | + >>> smoothsort([5, 4, 3, 2, 1]) |
| 147 | + [1, 2, 3, 4, 5] |
| 148 | + >>> smoothsort([3, 3, 2, 1, 2]) |
| 149 | + [1, 2, 2, 3, 3] |
| 150 | + >>> smoothsort([1, 2, 3, 4, 5]) |
| 151 | + [1, 2, 3, 4, 5] |
| 152 | + >>> smoothsort([-3, 0, -1, 5, 2]) |
| 153 | + [-3, -1, 0, 2, 5] |
| 154 | + """ |
| 155 | + n = len(seq) |
| 156 | + if n < 2: |
| 157 | + return seq |
| 158 | + |
| 159 | + # ``heap_sizes[i]`` is the Leonardo order of the i-th tree (left to right). |
| 160 | + heap_sizes: list[int] = [] |
| 161 | + |
| 162 | + # ------------------------------------------------------------------ |
| 163 | + # Phase 1 : Build the Leonardo heap forest over seq[0..n-1]. |
| 164 | + # ------------------------------------------------------------------ |
| 165 | + for i in range(n): |
| 166 | + # If the two rightmost trees have consecutive orders, merge them. |
| 167 | + if len(heap_sizes) >= 2 and heap_sizes[-2] == heap_sizes[-1] + 1: |
| 168 | + heap_sizes.pop() |
| 169 | + heap_sizes[-1] += 1 |
| 170 | + elif heap_sizes and heap_sizes[-1] == 1: |
| 171 | + heap_sizes.append(0) |
| 172 | + else: |
| 173 | + heap_sizes.append(1) |
| 174 | + |
| 175 | + _trinkle(seq, i, heap_sizes, len(heap_sizes) - 1) |
| 176 | + |
| 177 | + # ------------------------------------------------------------------ |
| 178 | + # Phase 2 : Extract maximum elements right-to-left. |
| 179 | + # ------------------------------------------------------------------ |
| 180 | + for i in range(n - 1, -1, -1): |
| 181 | + order = heap_sizes.pop() |
| 182 | + if order > 1: |
| 183 | + # Expose the two child roots and re-trinkle each. |
| 184 | + right_order = order - 2 |
| 185 | + left_order = order - 1 |
| 186 | + right_pos = i - 1 |
| 187 | + left_pos = i - 1 - _LEONARDO[right_order] |
| 188 | + |
| 189 | + heap_sizes.append(left_order) |
| 190 | + _trinkle(seq, left_pos, heap_sizes, len(heap_sizes) - 1) |
| 191 | + |
| 192 | + heap_sizes.append(right_order) |
| 193 | + _trinkle(seq, right_pos, heap_sizes, len(heap_sizes) - 1) |
| 194 | + |
| 195 | + return seq |
| 196 | + |
| 197 | + |
| 198 | +if __name__ == "__main__": |
| 199 | + import doctest |
| 200 | + import random |
| 201 | + |
| 202 | + results = doctest.testmod(verbose=False) |
| 203 | + assert results.failed == 0, f"{results.failed} doctest(s) failed" |
| 204 | + |
| 205 | + for trial in range(5000): |
| 206 | + sample = random.choices(range(-50, 50), k=random.randint(0, 30)) |
| 207 | + got = smoothsort(sample[:]) |
| 208 | + assert got == sorted(sample), ( |
| 209 | + f"Trial {trial}: smoothsort({sample!r}) -> {got!r}, " |
| 210 | + f"expected {sorted(sample)!r}" |
| 211 | + ) |
| 212 | + |
| 213 | + print("All doctests and 5 000 random trials passed.") |
0 commit comments