Сортировка слиянием

Сортировка слиянием один из тех редких алгоритмов, которые не устаревают, а напротив год от года развивается, получает новые публикации и исследования. Это действительно чудесный со многих сторон алгоритм, который находит применение в современном мире, а исследования, связанные с ним предлагают решения, применимые далеко за пределами сортировки. Это своего рода жемчужина среди всех алгоритмов сортировки.

Базовый алгоритм

Mergesort является двойственным к быстрой сортировке и относится к классу алгоритмов «Разделяй и властвуй». Основной функцией этого алгоритма является функция слияния, которая получает на вход два отсортированных массива и возвращает их отсортированное объединение. В простейшем варианте её можно записать так:

def merge(A, B):
    i, j, C = 0, 0, []
    while True:
        if A[i] < B[j]:
            C.append(A[i])
            i += 1
            if i == len(A):
                C.extend(B[j:])
                break

        else:
            C.append(B[j])
            j += 1
            if j == len(B):
                C.extend(A[i:])
                break

    return C

Имея такую функцию можно построить функцию сортировки двумя путями.

Top-Down и Bottom-Up

Например, можно поступить как в быстрой сортировке: делить массив пополам до тех пор, пока не получатся подмассивы длины \(1\). Такие подмассивы всегда отсортированы, поэтому к ним можно применить функцию слияния:

def top_down_merge_sort(A):
    if len(A) == 1:
        return A

    d = len(A) // 2
    left = top_down_merge_sort(A[:d])
    right = top_down_merge_sort(A[d:])

    return merge(left, right)

Конечно, в отличии от быстрой сортировки, глубина стека этой рекурсивной реализации гарантированно будет \(\log n\), но можно решить эту задачу вообще не используя рекурсию.

def bottom_up_merge_sort(A):
    k = 1
    while k < len(A):
        for i in range(0, len(A)-k, 2*k):
            A[i:i+2*k] = merge(A[i:i+k], A[i+k:i+2*k])
        k *= 2

    return A

То есть пойти снизу вверх, от одноэлементных массивов через массивы длинны \(k = 2^j\) до тех пор, пока сливаемая часть не будет содержать весь исходный массив.

Эти два подхода к решению задачи, как не трудно догадаться, носят названия Top-Down (сверху вниз) и Bottom-Up (снизу вверх). Подход сверху вниз обычно является рекурсивным. Но, если для какого-то алгоритма есть решение рекурсивное и нерекурсивное, однозначно никогда не нужно использовать рекурсивное. Собственно и здесь рекурсивное решение приведено исключительно как дань уважения быстрой сортировке.

В текущей реализации получаем следующие оценки сложности:

И по памяти: \(S(n) = 2n\).

Экономим место

Вернёмся к функции merge. В текущей реализации эта функция требует \(n\) дополнительной памяти, потому что результат пишется в совершенно новый массив. Впрочем, если обратить внимание на функцию bottom_up_merge_sort, результат хорошо бы получить на месте исходных массивов. Попробуем сократить аппетиты функции merge, при учёте того, что сливаемые массивы идут в памяти друг за другом. Пусть даны два массива \(A\) и \(B\), являющиеся по сути двумя подмассивами одного большего массива, идущие друг за другом. А также дан пустой участок памяти, размером \(|A|\). Тогда несложно видеть, что, скопировав \(A\) в свободный участок памяти, можно записать ответ поверх начального местоположения \(A\) и \(B\). Для этого воспользуемся тремя указателями: *A — текущая голова копии \(A\), *B — текущая голова массива \(B\) и *Res — указатель на конец результата слияния.

При этом очевидно, что указатель *Res «догонит» указатель *B только тогда, когда указатель *A исчерпает весь массив \(A\). Также несложно адоптировать этот алгоритм для переноса в дополнительную память массива \(B\). Что, кстати сказать, будет выгоднее для нашей реализации Bottom-Up (?).

Итак, пусть AB — массив, в котором располагаются элементы массива \(A\), а затем элементы массива \(B\). n — индекс первого элемента массива \(B\). И C — пустой массив, размера \(B\).

def merge(AB, n, C):
    C[:] = AB[n:]

    a, b, r = n - 1, len(C) - 1, len(AB) - 1
    while True:
        if AB[a] > C[b]:
            AB[r] = AB[a]
            a -= 1
            if a < 0:
                AB[a+1:r] = C[:b+1]
                return
        else:
            AB[r] = C[b]
            b -= 1
            if b < 0:
                return
        r -= 1

А можно вообще без дополнительной памяти?

Слияние без дополнительной памяти сделать не получится, а вот всю сортировку целиком можно, хотя и не нужно (сложность будет всё также \(O(n\log n)\), вот только откидываемый скалярный коэффициент будет слишком большой). Для начала стоит отметить, что для слияния массивов разной длинны описанным выше способом объём дополнительной памяти требуется соразмерно меньшему из сливаемых массивов. Кроме того, если заменить все операции на обмены, то нетрудно видеть, что можно сохранить данные в дополнительной памяти (потеряв при этом их порядок). То есть, при перемещении \(A\) в дополнительную память данные из дополнительной памяти оказываются на месте \(A\) (\(a_i \leftrightarrow c_i\)). А при слиянии элементы из дополнительной памяти будут обмениваться с элементами \(A\) и \(B\), но в итоге вернутся в дополнительную память. Таким образом мы можем использовать часть оригинального массива как дополнительную память при сортировке другой его части. Не стоит забывать, что порядок всё же теряется, поэтому поступим следующим образом:

1. отсортируем правую половину массива, воспользовавшись левой как дополнительной памятью;
2. отсортируем левую неотсортированной части, воспользовавшись правой половиной неотсортированной части;
3. сольём отсортированные части, записывая результат поверх неотсортированной части;
4. повторим шаги 2 и 3 до тех пор, пока в неотсортированной части не останется одноэлементный массив, который мы сольём с остатком при помощи алгоритма вставки.

Оптимизируем по времени

В текущей реализации функция merge для двух списков \(A\) и \(B\) длин \(n\) и \(m\) соответственно имеет линейную сложность, то есть \(O(n+m)\). Это можно оптимизировать несколькими очевидными способами, но прежде чем это делать, давайте определимся с функцией, к которой будем стремиться. Теоретически можно оценить функцию слияния двух отсортированных массивов способом похожим на доказательство минимальной сложности сортировки.

Немного теории

Итак, пусть у нас есть два отсортированных массива: \[ A = ( a_i \colon i < n ) \quad\text{и}\quad B = ( b_j \colon j < m ). \]

Заметим, что в результирующем массиве \(a_i\) будет стоять левее \(a_{i+1}\) и аналогично с элементами \(B\). Таким образом это сводится к задаче Шаров и перегородок. То есть всевозможных результатов слияния будет \[ k = \tbinom{n + m}{m} = \tbinom{n + m}{n}. \]

Если на каждом шаге алгоритма мы делаем бинарный выбор и дерево выбора идеально сбалансированно, то его глубина будет \(\log k\). Следовательно минимальное число сравнений необходимое для слияния двух списков: \[ T(n, m) = O\left(\log \tbinom{n + m}{m}\right). \]

Применив формулу Стирлинга и упростив выражение по \(O\), получим следующую оценку: \[ T(n, m) = O\left(n \log \frac{m}{n} + n\right). \]

При этом без ограничения общности можем считать, что \(n < m\).

Эта оценка совпадает с наивной в крайних случаях: \[ T(1, m) = O(\log m) \quad\text{и}\quad T(n, n) = O(n). \] То есть при \(n = 1\) это можно свести к бинарному поиску места вставки одного элемента в массив \(B\). А для списков равной длины получаем линейную сложность.

А на практике как?

Можно ли получить подобную оценку на практике? Можно. Для этого нам понадобится немного видоизменить бинарный поиск. Пусть нам даны отсортированный массив \(B = ( b_i \colon i < m )\) и элемент \(x\). Надо найти \(i < m\) такое, что \(b_i \leqslant x\) и \(x \leqslant b_{i+1}\), если \(i+1 < m\).

Вместо того, чтобы делать 2 указателя и итеративно их сближать, сделаем следующим образом:

1. \(k = 0\);
2. пока \(2^k - 1 < m\) и \(b_{2^k-1} < x\) увеличиваем \(k\) на \(1\);
3. если \(2^k - 1 \geqslant m\), то \(l = m - 1\), иначе \(l = 2^k - 1\);
4. пока \(l \geqslant 0\) и \(b_l > x\) уменьшаем \(l\) на \(1\);
5. \(l + 1\) — искомая позиция;

Асимптотически такой поиск будет иметь логарифмическую сложность, при этом нетрудно заметить, что \(T\) будет равно не \(O(\log m)\), а \(O(\log i)\). Ведь мы не рассматриваем элементы правее \(2^{\lceil \log i \rceil}\).

Также надо помнить про то, что мы будем последовательно искать места для вставки элементов монотонной последовательности, так что начинать поиск надо не с начала массива \(B\), а с места предыдущей вставки. Таким образом для вставки \(n\) элементов между \(m\) элементами потребуется: \[ T(n, m) = O(\log i_0 + \log (i_1 - i_0) + \ldots + \log (i_{n-1} - i_{n-2})), \] где \(i_j\) — места для вставки \(j\)-го элемента.

Теперь поделим всё на \(n\) и вспомним, что функция \(\log\) выпукла вверх, а значит можно применить неравенство Йенсена: \[ \frac{T(n, m)}{n} = O\left(\frac{\log i_0 + \log (i_1 - i_0) + \ldots + \log (i_{n-1} - i_{n-2})}{n}\right) \leqslant \] \[ \leqslant O\left(\log\left(\frac{i_0 + (i_1 - i_0) + \ldots + (i_{n-1} - i_{n-2})}{n}\right)\right) = \] \[ = O\left(\log\left(\frac{i_{n-1}}{n}\right)\right) \leqslant O\left(\log\left(\frac{m}{n}\right)\right). \]

В получившейся формуле есть один недостаток: при \(n = m\) мы получаем \(\log 1\) равный \(0\). Ноль это несколько излишне оптимистичная оценка, поэтому добавим \(1\) на этот случай (?). Домножим обратно на \(n\) и получим нашу теоретическую оценку: \[ T(n, m) = O(n \log \frac{m}{n} + n). \]

Такое слияние называется галопирование (galloping). Но ближе к коду. В данном примере копировать в дополнительную память будем \(A\), а проходить массивы будем слева направо.

def galloping(AB, n, C):
    C[:] = AB[:n]

    # r — указатель на конец результата
    # j — место последней вставки
    # m — длина остатка B
    r, j, m = 0, n, len(AB) - n
    for i in range(n):
        # k — степень двойки
        # l — указатель на 2^k-1 элемент
        k, l = 0, 0
        while l < m and AB[j+l] < C[i]:
            k += 1
            l = 2**k - 1

        if l >= m:
            l = m - 1

        while l >= 0 and AB[j+l] > C[i]:
            l -= 1

        l += 1
        AB[r:r+l], AB[r+l] = AB[j:j+l], C[i]
        r, j, m = r + l + 1, j + l, m - l

Оптимальный порядок слияния

С отдельными слияниями разобрались, но что с сортировкой в целом? Первое, что стоит сделать, это улучшить случай отсортированного массива. Для этого необходимо не просто разрезать массив на подмассивы одинаковой длинны \(2^j\), а выбирать максимально длинные отсортированные части. Кроме того можно выбирать не только подмассивы, отсортированные в нужном нам порядке, но и обратно отсортированные, меняя для них направление итератора. Таким образом можно получить список подмассивов оригинального массива, для каждого подмассива необходима пара чисел: начало и конец. При этом, если начало больше конца, значит подмассив обратно отсортирован.

def chunking(A):
    chunks = []
    a, d = 0, 0
    for b in range(1, len(A)):
        if d == 0:
            d = A[b] - A[a]
            continue

        if (A[b] - A[b-1])*d < 0:
            chunks.append((a, b-1) if d > 0 else (b-1, a))
            a, d = b, 0

    chunks.append((a, b) if d > 0 else (b, a))

    return chunks

И тут возникает следующий вопрос: в каком порядке сливать полученные подмассивы? Для того, чтобы более наглядно фиксировать порядок слияния, будем использовать дерево слияния.

Ниже в узлах дерева будем записывать только размеры подмассивов.

В данном случае это минимальное по глубине дерево. Такое дерево легко получить при обработке подмассивов в порядке расположения. Однако, это не всегда выигрышный подход, как и в этом примере. Надо помнить, что сложность функции слияния зависит от размеров сливаемых частей. И это дерево не является оптимальным, потому что самый большой подмассив, состоящий из \(16\) элементов мы сливали дважды. Более оптимальным в данном случае будет следующее дерево слияния:

Лучшим деревом слияния будет то, где каждый раз сливаются массивы минимальной длинны из оставшихся. Для его построения воспользуемся следующей леммой:

Лемма

Пусть \(X = ( x_i \colon i < k )\) конечная последовательность натуральных чисел. Будем строить последовательность \(Y\) следующим образом:

• выберем и вычеркнем из \(X \cup Y\) два минимальных числа \(a\) и \(b\);
• добавим в \(Y\) число \(a + b\);
• повторяем до тех пор, пока \(|X \cup Y| \geqslant 2\).

Тогда на каждом шаге \(a + b \geqslant \max Y\).

Доказательство достаточно тривиально и оставляется на откуп читателю (?).

Куда труднее эту лемму применить. Для оптимального слияния необходимо правильно расположить начальные подмассивы друг относительно друга так, чтобы всегда сливались соседние. Расположить подмассивы просто в порядке возрастания их длин неверное решение, например, если у нас есть подмассивы длин \((1, 2, 2, 2, 3)\), то последний подмассив надо перенести в середину:

Релизация подобной перестановки оставляется читателю.

TimSort

Все описанные оптимизации нашли применение в гибридном алгоритме TimSort. Помимо этого там применяется сортировка вставками для того, чтобы избежать слишком коротких начальных подмассивов, а также оригинальный способ строить дерево слияния близкое к идеальному без перестановки подмассивов, используя стек. Этот алгоритм является стандартным для языков Python и Java, а также имеет реализации для других языков.

Вопросы

1. ^ Почему для приведённой Bottom-Up реализации выгоднее копировать в дополнительную память именно массив \(B\)?
2. ^ Почему мы вольны добавить \(1\) при оценке сложности, не теряя строгости переходов?
3. ^ Докажите Лемму 1.

Задания

1. Напишите сортировку слиянием, использующую \(n / 2\) дополнительной памяти.
2. Напишите функцию-итератор, возвращающую подмассивы слияния в оптимальном порядке.
3. * Напишите функцию, переставляющую подмассивы так, что в оптимальном дереве слияний всегда будут сливаться соседние подмассивы.

Дополнительное чтение

1. Peters, Tim. “Original original explanation of TimSort.” Python bug tracker. 2002. https://bugs.python.org/file4451/timsort.txt
2. Auger, Nicolas, Nicaud, Cyril, and Pivoteau, Carine. “Merge Strategies: from Merge Sort to TimSort.” HAL UPEC - UPEM Portal. 2015. https://hal-upec-upem.archives-ouvertes.fr/hal-01212839v2/document

Ссылки

1. Максим Александрович Бабенко. Видеолекции курса «Алгоритмы и структуры данных». https://yandexdataschool.ru/edu-process/courses/algorithms#item-3