Задача коммивояжёра

Коммивояжёр (travelling salesman) хочет посетить каждый из $n$ городов по одному разу и вернуться в свой родной город. Он знает расстояния между всеми парами городов. Какой порядок посещения городов оптимален?

Более формально: дан полный (есть ребро между каждой парой вершин) взвешенный неориентированный граф на $n$ вершинах; необходимо найти гамильтонов цикл (цикл, проходящий по каждой вершине ровно один раз), имеющий минимально возможную длину. Расстояние между вершинами $u$ и $v$ будем обозначать как $d_{u, v}$.

Всего различных маршрутов $(n-1)$ !, поэтому перебор их всех возможен лишь при небольших $n$. Динамическим программированием эту задачу можно решить за $O\left(2^{n} \cdot n^{2}\right)$, что уже значительно быстрее (хотя оценка времени всё ещё экспоненциальна).

Естественная подзадача в данном случае - задача нахождения “префикса” маршрута (первых нескольких рёбер маршрута, начиная из родного города коммивояжёра; будем считать, что этот город имеет номер 0). Как можно закодировать начало маршрута? Для того, чтобы префикс маршрута можно было попытаться продлить (или, наоборот, рассмотреть префикс без последнего ребра), необходимо знать, какие вершины входят в префикс, а также какая из них является последней. Таким образом, состояние параметризуется парой $(A, v)$, где $A$ - подмножество вершин графа, содержащее вершину 0 , $v \in A$ - номер последней вершины в префиксе маршрута.

Начальные состояния: $l[\{0\}, 0]=0 ; l[A, 0]=\infty$ при $|A|>1$. Для перехода нужно перебрать последнее ребро префикса маршрута: пусть $v \neq 0, v \in A$, тогда

$$ l[A, v]=\min _{\substack{u \in A, u \neq v}}\left(l[A \backslash\{v\}, u]+d_{u, v}\right) $$

Для получения ответа на задачу нужно перебрать последнюю вершину маршрута: ответ равен

$$ \min _{v>0}\left(l[A, v]+d_{v, 0}\right) . $$

Прежде, чем перейти к реализации, нужно научиться использовать множества в качестве индексов массива, а также быстро производить над множествами различные операции (например, удаление элемента из множества).

Работа с множествами

Пронумеруем все подмножества $n$-элементного множества $U=\{0, \ldots, n-1\} n$-битными числами: $A \subset U$ сопоставим $\operatorname{mask}(A)=\sum_{i \in A} 2^{i}$. Номер множества удобно использовать в качестве индекса массива; также с помощью битовых операций над числами можно быстро

осуществлять различные операции над множествами с соответствующими номерами. Приведём несколько примеров:

$\operatorname{mask}(A \cup B)=\operatorname{mask}(A) \mid \operatorname{mask}(B)$, где $\mid$ - битовое “ИЛИ”;
$\operatorname{mask}(A \cap B)=\operatorname{mask}(A) \& \operatorname{mask}(B)$, где & - битовое “И”;
$\operatorname{mask}(A \oplus B)=\operatorname{mask}(A) \oplus \operatorname{mask}(B)$, где $\oplus$ слева - симметрическая разность множеств, а $\oplus$ справа - битовый “XOR”;
$\operatorname{mask}(A \backslash B)=\operatorname{mask}(A) \& \sim \operatorname{mask}(B)$, где $\sim$ - битовое отрицание;
$\operatorname{mask}(U)=(1 \ll n)-1$, где $\ll k$ - битовый сдвиг на $k$ бит влево;
$\operatorname{mask}(\{k\})=1 \ll k$;
$(\operatorname{mask}(A) \gg k) \& 1$ - выражение, равное единице тогда и только тогда, когда $k \in A$ (здесь $\gg k$ - битовый сдвиг на $k$ бит вправо).

В C++ все битовые операторы, за исключением “XOR”, обозначаются теми же символами, что и выше; ^ - битовый “XOR” в C++. Заметим, что битовые операции над стандартными числовыми типами осуществляются за $O(1)$. Мы для простоты будем считать, что размер множества позволяет использовать стандартный тип данных для хранения номеров подмножеств, поэтому все приведённые выше выражения можно вычислить за $O(1)$.

Размер подмножеств

В качестве простого примера вычислим размер каждого подмножества $A \subset U$. Размер множества совпадает с количеством ненулевых бит в его номере. Количество единичных бит в числе $x$ равняется сумме количества единичных бит в $\left\lfloor\frac{x}{2}\right\rfloor$ и значения младшего бита $x$.

1
2
3


size[0] = 0 # размер пустого множества равен 0
for mask = 1..((1 << n) - 1):
    size[mask] = size[mask >> 1] + (mask & 1)

✍️ В компиляторе C++ GCC есть встроенная функция $\texttt{__builtin_popcount}(x)$, возвращающая число бит в числе $x$ типа $\texttt{int}$ (есть подобные функции и для других стандартных типов); она работает за $O(1)$, но значительно медленнее, чем, например, элементарные битовые операции.

Сумма элементов в подмножестве

Пусть есть $n$ чисел $a_{0}, \ldots, a_{n-1}$; необходимо для каждого подмножества $A \subset\{0, \ldots, n-1\}$ вычислить $s(A)=\sum_{i \in A} a_{i}$. Можно вычислить каждую сумму по определению:

for mask = 0..((1 << n) - 1):
    s[mask] = 0
    for i = 0..(n - 1):
        if ((mask >> i) & 1): # i лежит в множестве с номером mask
            s[mask] += a[i]

Это потребует $O\left(2^{n} \cdot n\right)$ времени. Можно справиться быстрее с помощью простого динамического программирования: пусть $x \in A$ - какой-нибудь элемент $A$, тогда $s(A)=$ $s(A \backslash\{x\})+a_{x}$. При этом, поскольку мы вычисляем значения $s$ в порядке возрастания номеров множеств, а $\operatorname{mask}(A \backslash\{x\})<\operatorname{mask}(A), s(A \backslash\{x\})$ уже будет известно к моменту вычисления $s(A)$. В качестве $x$ удобно использовать максимальный элемент $A$ : такой элемент соответствует старшему биту $\operatorname{mask}(A)$, номер которого может только возрастать при возрастании номера. Получаем следующее решение за $O\left(2^{n}\right)$ :

1
2
3
4
5
6


x = 0
s[0] = 0
for mask = 1..((1<< n) - 1):
    if mask == (1 << (x + 1)): # номер старшего бита увеличился
        x += 1
    s[mask] = s[mask - (1 << x)] + a[x]

Сумма по подмножествам

Пусть теперь для каждого подмножества $A \subset\{0, \ldots, n-1\}$ дано число $f_{A}$, и требуется для каждого $A$ вычислить $s(A)=\sum_{B \subset A} f_{B}$.

Можно для каждого $A$ перебрать все возможные $B$ и проверить, какие из них являются подмножествами $A$; получится решение за $O\left(2^{n} \cdot 2^{n}\right)=O\left(4^{n}\right)$.

1
2
3
4
5


for A = 0..((1 << n) - 1):
    g[A] = 0
    for B = 0..((1<< n) - 1):
        if (A & B) == B: # B - подмножество A
            g[A] += f[B]

Можно действовать аккуратнее, и для фиксированного $A$ сразу перебирать только $B$, являющиеся подмножествами $A$. Будем перебирать подмножества $A$ в порядке уменьшения их номеров. Пусть $A=\left\{a_{0}, \ldots, a_{k-1}\right\}, a_{0}<\cdots < a_{k-1}$. Если мы мысленно сдвинем разряды $a_{0}, \ldots, a_{k-1}$ на позиции $0, \ldots, k-1$, то номера подмножеств $A$ перейдут в числа от 0 до $2^{k}-1$. Значит, мы хотим перебрать числа от 0 до $2^{k}-1$ в порядке убывания, только со сдвинутыми разрядами: нулевой разряд нужно сдвинуть в $a_{0}$-й,… $(k-1)$-й разряд — в $a_{k-1}$-й.

Для того, чтобы по числу $i$ получить число $i-1$, нужно занулить младший единичный бит $i$, а во все биты младше его записать единицу. Значит, для того, чтобы получить по подмножеству $B \subset A$ следующее подмножество $A$ в порядке убывания номеров, нужно удалить из $B$ минимальный элемент $x \in B$ и добавить все элементы из $A \backslash B$, меньшие $x$. Это соответствует формуле $(\operatorname{mask}(B)-1) \& \operatorname{mask}(A)$.

1
2
3
4


for A = 0..((1<< n) - 1):
    g[A] = f[0] # учли пустое подмножество
    for (B = A; B > 0; B = (( B - 1) & A)): # все непустые подмножества A
        g[A] += f[B]

Время работы полученного решения пропорционально количеству пар $(A, B)$ таких, что $B \subset A$. Количество таких пар равняется $3^{n}$. Это можно понять следующим способом: множество $A$ размера $k$ имеет $2^{k}$ различных подмножеств, количество множеств размера $k$ равно $\binom{n}{k}$, значит, всего пар $(A, B)$

$$ \sum_{k=0}^{n}\binom{n}{k} \cdot 2^{k}=\sum_{k=0}^{n}\binom{n}{k} \cdot 2^{k} \cdot 1^{n-k}=(2+1)^{n}=3^{n} $$

Альтернативный способ: если пара $(A, B)$ зафиксирована, то каждый элемент находится в одном из трёх состояний: лежит в $A$, лежит в $A \backslash B$ или же не лежит в $A$. Всего элементов $n$, значит, всего $3^{n}$ вариантов.

Итак, мы получили решение за $O\left(3^{n}\right)$. Можно воспользоваться динамическим программированием и решить задачу ещё быстрее. Пусть

$$ g[A, k]=\sum_{\substack{B \subset A, A \backslash B \subset\{0, \ldots, k-1\}}} f_{B} $$

Нас интересуют значения $g[A, n]=g(A)$. Начальные значения: $g[A, 0]=f_{A}$. Переход:

$$ g[A, k+1]= \begin{cases}g[A, k], & \text { если } k \notin A, \\ g[A, k]+g[A \backslash\{k\}, k], & \text { если } k \in A .\end{cases} $$

Почему это так? $g[A, k+1]$ - сумма $f_{B}$ по подмножествам $B \subset A$, которые могут отличаться от $A$ только принадлежностью элементов $0, \ldots, k$. Если $k \notin A$, все такие $B$ отличаются от $A$ только принадлежностью $0, \ldots, k-1$, значит, $g[A, k+1]=g[A, k]$. Если же $k \in A$, то все такие $B$ делятся на две группы: те, которые содержат $k$, и те, которые не содержат. Сумма по первой группе равна $g[A, k]$. Все множества во второй группе являются подмножествами $A \backslash\{k\}$; сумма по второй группе равна $g[A \backslash\{k\}, k]$.

Заметим, что при вычислении состояния $(A, k+1)$ нас интересуют только состояния $(A, k)$, и, возможно $(A \backslash\{k\}, k)$. Однако мы пользуемся значением $g[A \backslash\{k\}, k]$ только в случае, когда $g[A \backslash\{k\}, k]=g[A \backslash\{k\}, k+1]$. Значит, можно обойтись одномерным массивом вместо двумерного. Получаем решение за $O\left(2^{n} \cdot n\right)$, использующее $O\left(2^{n}\right)$ памяти.

1
2
3
4
5
6


for mask = 0..((1 << n) - 1):
    g[mask] = f[mask]
for k = 0..(n - 1):
    for mask = 0..((1 << n) - 1):
        if (mask >> k) & 1:
            g[mask] += g[mask - (1 << k)]

Реализация решения задачи коммивояжёра

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14


fill(l, inf)
l[1, 0] = 0 # 1 - номер множества 0
# номера множеств, содержащих 0 - все нечётные числа
for (mask = 1; mask < (1 << n); mask += 2):
    for v = 1..(n - 1): # перебираем последнюю вершину префикса марщрута
        if ((mask >> v) & 1) == 0: # v не лежит в префиксе
            continue
        for u = 0..(n - 1): # перебираем второй конец последнего ребра
            if u == v or ((mask >> u) & 1) == 0:
                continue
            l[mask, v] = min(l[mask, v], l[mask - (1 << v), u] + d[u, v])
res = inf
for v = 1..(n - 1):
    res = min(res, l[(1 << n) - 1, v] + d[v, 0])

Восстановление ответа осуществляется стандартными методами.