Pull to refresh

Метод динамического программирования для подсчёта числа циклов на прямоугольной решетке

Reading time 11 min
Views 13K
Эта статья адресована тем читателям, кто занимается программированием алгоритмов, и особенно интересуется труднорешаемыми задачами. Тем хабралюдям, которые против размещения алгоритмов на Хабре следует немедленно прекратить читать данную работу.

В статье я покажу как использовать метод динамического программирования по профилю для решения задачи о подсчёте количества гамильтоновых циклов на прямоугольной решётке размером m на n. На Хабре есть несколько статей, посвященных теме динамического программирования (например, эта), но нигде не идёт речь о более сложном применении метода. Данный подход также можно называть методом матрицы переноса, кому как нравится.

Предупреждаю, что статья содержит около 2000 слов (8 страниц А4), но дорогу осилит идущий.



Краткое содержание



  1. Немного о труднорешаемых задачах
  2. Определения и постановка задачи
  3. Разрезы и конечный автомат
  4. Метод матрицы переноса
  5. Решение без использования матрицы
  6. Сложность алгоритма
  7. Модулярная арифметика
  8. Обсуждения и результаты
  9. Список использованных источников


1. Немного о труднорешаемых задачах



Что такое труднорешаемая задача? Это такая задача, которая не имеет эффективного алгоритма решения. Причём это не обязательно задача с экспоненциальной сложностью решения, сложность может быть и полиномиальной, но с достаточно большой степенью полинома (например, задача о k ферзях, о которой я уже писал здесь). Больше о таких задачах и их классификации можно узнать из классической книги М. Гэри, Д. Джонсон «Вычислительные машины и труднорешаемые задачи».

Такие задачи обычно решаются методами перебора, динамического программирования и комбинаторного анализа. Часто приемлемое решение получается при использовании этих трёх составляющих одновременно, но к этому нужен ещё и опыт написания достаточно оптимального кода. Совсем хорошо, если под рукой есть кластер хотя бы на несколько тысяч ядер и хотя бы пара сотен гигабайт оперативной памяти. К сожалению, я не шучу, многие задачи могут быть решены только на такой машине…

Мне нравится решать труднорешаемые задачи и я даже время от времени устраиваю конкурсы для программистов-математиков на данную тему. Эта статья частично приурочена к конкурсу «Предгамильтоновы циклы», который начался 1 октября и закончится 31 октября 2010 г. Задача, предлагаемая в условиях конкурса, решается мной таким же методом (с незначительными изменениями), который будет описан здесь.

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


2. Определения и постановка задачи



В работе встречается граф, называемый прямоугольной решёткой Pm×Pn. Прямоугольная решётка – это целочисленные точки координатной плоскости, связанные по принципу ближайшего соседства. Обозначение Pm×Pn связано с тем, что целочисленный прямоугольник суть прямое произведение двух цепей, первая из которых имеет m звеньев (обозначается Pm), а вторая – n (обозначается Pn). Число m – ширина решётки (протяжённость слева направо), а n – её длина (протяжённость снизу вверх).

Гамильтоновым циклом в графе называется цикл, проходящий через каждую его вершину в точности один раз. На рис. 1 представлен пример гамильтонового цикла на решётке P6×P8.


Рисунок 1 – Пример гамильтонового цикла на решётке размером 6×8

Задача состоит в том, чтобы по заданным m, n ≥ 2 отыскать количество гамильтоновых циклов на решётке Pm×Pn. Например, на следующем рисунке изображены все 6 гамильтоновых циклов на решётке размером 4×4.


Рисунок 2 – Все 6 гамильтоновых циклов на решётке размером 4×4


3. Разрезы и конечный автомат



Любой цикл в графе Pm×Pn можно строить последовательно по «слоям» Pm×Pk (k=0,1,2,…,n). Для наглядности, на рис. 3 изображена линия, разрезающая граф на две части. Снизу от неё граф P6×P5, на котором еще недостроенный «цикл» представляет собой множество непересекающихся простых цепей. Эту построенную часть цикла будем называть «прошлым». Сверху от линии изображён один из возможных вариантов завершения недостроенного «цикла» (это «будущее»). Цикл замкнут и не касается линии разреза в точках своего перегиба. А поскольку, обходя цикл в любом направлении, мы вернёмся в исходную точку, линия разреза будет пересечена циклом чётное число раз.


Рисунок 3 – Разрез цикла

Само пересечение цикла линией разреза можно представить в виде правильной скобочной последовательности, разряжённой нулями. Например, на рис. 3 такая последовательность имеет вид (0000). В этой последовательности пара соответствующих друг другу открывающей и закрывающей скобок соотносятся с двумя концам одной и той же цепи, походящей через «прошлое», а нули означают, что в данных точках нет пересечения линии разреза с циклом. Таким образом, любой разрез можно зашифровать правильной скобочной последовательностью, разряжённой нулями. Поскольку начальный и конечный разрезы должны быть нулевыми (состоящими целиком из нулей), мы должны различать их, обозначая один с чертой снизу, а второй – с чертой сверху.

В процессе построения цикла нам не нужно знать ни «прошлое», ни «будущее», а достаточно хранить только информацию о том, какой разрез в данный момент времени получился, и как из данного разреза можно получить остальные. Действительно, как переходить от текущего разреза к следующим? Например, как на рис. 3 получается переход от разреза (0000) к следующему разрезу 00()00? Это будет показано на примерах.


Рисунок 4 – Примеры перехода от одного разреза к другому

На рис. 4 изображены 3 примера переходов от одного разреза к другому. Дугообразные линии показывают, что концы соответствующих им вертикальных дуг соединяются через «прошлое». Первый случай как раз показывает переход на рис. 3 Остальные два случая чуть более сложные, но наглядно демонстрируют, какие новые разрезы могут получаться путём добавления разных комбинаций горизонтальных дуг.

Однако не все комбинации горизонтальных дуг могут привести к новому разрезу. Некоторые комбинации могут оказаться недопустимыми по трём причинам, обозначенным на рис 5.


Рисунок 5 – Примеры недопустимых переходов. Три возможных случая

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

Итак, предположим, что мы располагаем информацией обо всех возможных разрезах нашей решётки и о том, какие разрезы из каких могут быть получены путём включения горизонтальных дуг. Всю эту информацию можно представить в виде конечного автомата, узлы которого будут разрезами, а дуги будут показывать на возможность перехода. На рис. 6 представлен пример такого конечного автомата, построенного для решётки размером 3×n.


Рисунок 6 – Конечный автомат для решётки размером 3×n

Ответом на вопрос о том, сколько существует циклов на решётке P3×Pn будет число способов перейти от начального состояния автомата к финальному за n шагов. Как это сделать? Я знаю два способа: плохой и получше. Начну с плохого.


4. Метод матрицы переноса



Один из способов решения поставленной задачи называется методом матрицы переноса, который состоит в следующем. Пронумеруем все возможные разрезы от 1 до N любым способом, но так, чтобы начальное состояние имело номер 1, а конечное – номер N. Построим матрицу T, в которой числа Tij будут равны количеству способов получить из разреза с номером i разрез с номером j. Решением задачи для решётки Pm×Pn будет число (Tn)1,N.
Например, для решётки P3×Pn мы уже построили автомат (рис. 6), поэтому можем выписать матрицу переноса.


Рисунок 7 – Матрица переноса, или матрица смежности конечного автомата рис. 6 и её шестая степень

На рис. 7 изображена матрица переноса T для гамильтоновых циклов на решётке P3×Pn. Решением задачи будет число (Tn)15. Например, (T6)15=4, таково и количество гамильтоновых циклов на решётке P3×P6.

Однако на практике такой способ решения задачи применим, самое большее, для m=12, когда матрица переноса для гамильтоновых циклов на решётке имеет размер 3774×3774 (см. табл. 1 в разделе с обсуждениями) и умещается в оперативной памяти персональной ЭВМ. В самой матрице хранятся длинные числа, поскольку, например, ответ для решётки P12×P20 содержит 33 цифры, а для решётки P12×P100 количество цифр достигает 174.

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


5. Решение без использования матрицы



Данный способ решения, как это не странно, тоже называется методом матрицы переноса, но матрица в нём явно не фигурирует, поэтому в русскоязычной литературе такой способ принято называть методом динамического программирования [ по конечному автомату ].

Смысл метода заключается в хранении количества способов добраться до каждой вершины конечного автомата за k шагов. На шаге k=0 это количество равно 0 во всех вершинах, кроме стартовой. Число способов оказаться в стартовой вершине на шаге k=0 по определению полагается равным 1. Пусть на некотором на шаге k в вершине i можно находиться P способами, тогда, переходя к шагу k+1 мы добавляем число P к числу способов находится во всех вершинах j, в которые можно попасть из i. Но для этого нужно знать, в какие вершины можно попасть из i, то есть опять хранить матрицу! На самом деле не нужно: попадая в вершину i, и зная, какой разрез она обозначает, мы заново строим все такие разрезы, которые могут из него получиться.

Как это реализуется в программе – дело индивидуальное. Поясню свой подход. Пусть имеется очередь Q, в которой хранятся разрезы (в виде скобочной последовательности) и число способов оказаться в них. В начале алгоритма в очереди хранится разрез, состоящий из нулей и число способов, равное 1. Вместе с очередью имеется хэш-таблица, которая тоже хранит разрезы и число способов оказаться в них. В начале алгоритма таблица пустая.

На шаге k достаём из очереди очередной разрез i, и число способов P оказаться в нём. Строим множество J всех возможных разрезов, получаемых из i. Для каждого j из J проверяем, находится ли оно в хэш-таблице. Если нет – добавляем j в таблицу вместе с числом способов P. Если да, то из таблицы берём число способов Q, которым можно придти в j, и записываем на его место число P+Q. Исключение составляет тот случай, когда j равно нулевому разрезу, это означает, что мы попали в финальное состояние и должны прибавить к ответу число P, но не добавлять нулевой разрез в хэш.

В конце шага k очередь пуста, поэтому перемещаем данные из хэш-таблицы в очередь, хэш очищаем и переходим на шаг k+1.

Я использую хэш с открытой адресацией и двойным хэшированием. Число коллизий равно 0.6 на каждое обращение. То есть в 60% случаях, делается один шаг по таблице, а в остальных случаях попадает сразу. Я считаю это хорошим хэшем для такой задачи.




6. Сложность алгоритма



Сначала я укажу на верхнюю границу сложности, а потом покажу что реально получается на практике.

Множество всех правильных скобочных последовательностей, разряжённых нулями и имеющих длину m являются словами Моцкина длиной m над алфавитом {0, (, )}, порождаемыми контекстно-свободной грамматикой

Z → Y | Y(Z) Z,
Y → ε | 0Y,

где ε обозначает пустое слово. Количество таких слов можно посчитать из комбинаторных соображений:



В этой формуле k – число пар скобок. Первое выражение под знаком суммы (дробь и самый первый биномиальный коэффициент) – это число Каталана, которое показывает количество правильных скобочных последовательностей, а второй биномиальный коэффициент – число способов поставить 2k скобок в m позиций (заполнив остальные позиции нулями).

Таким образом, количество всех возможных разрезов для гамильтоновых циклов на решётке, по крайней мере, не превышает количество слов Моцкина, а слова Моцкина растут со скоростью почти 3m. На самом деле, таких разрезов гораздо меньше, поскольку некоторые могут оказаться недопустимыми. Например, недопустимыми для решётки шириной 6 будут такие разрезы как (0)000 или (()0)0, поскольку их невозможно получить при построении гамильтонового цикла. В разделе с обсуждениями будет показана табл. 1, в которой указано общее число разрезов для некоторых m. В реальности получается, что только треть слов Моцкина допустимы, а с учетом симметрии, достаточно хранить только пятую их часть.

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

Для каждого разреза мы должны перебрать 2m-1 способов поставить горизонтальные дуги. И все это вместе мы должны проделать n раз.

Итого, алгоритм имеет сложность O(6m·m-1/2·n). А на самом деле, из-за того, что многие разрезы недопустимы, на практике сложность алгоритма оказывается равной O(5m·n). Так получилось в ходе моего наблюдения за алгоритмом.

Что касается используемой памяти, то основную часть ресурсов занимает очередь и хэш таблица. Каждый разрез для m≤32 можно хранить в 64-битовом целом числе, так как необходимо 2 бита на хранение скобок и числа 0. Плюс длинная арифметика, которая съедает все ресурсы окончательно. Но существует один приём, позволяющий замедлить вычисления, отказавшись от длинной арифметики.


7. Модулярная арифметика



Давайте вычислять ответ по модулю различных простых чисел: 231-1 = 2147483647, 2147483629, 2147483587, … Для этого программу необходимо запускать столько раз, сколько различных модулей мы выберем (хотя, конечно, можно считать сразу по 2-3 модулям, покуда влезает в память). По Китайской теореме об остатках можно восстановить ответ с точностью до произведения всех используемых простых чисел. Если произведение всех простых чисел больше, чем ответ к задаче, то ответ восстанавливается однозначно. Как выбрать нужное количество модулей?

Например, давайте считать число циклов на решётке размером 16 на n по двум модулям. Пусть n=1,2,…,16. Тогда получим ответы

[0, 1, 128, 405688, 24980352, 776813457, 729683652, 1087605227, 2000673777, 456710131, 1550214608, 568568229, 2047094091, 1175631455, 380271385, 1536681549] (mod 2147483647).

[0, 1, 128, 405688, 24980352, 776813727, 729709644, 1107434405, 301217473, 1373982040, 103268356, 218837622, 1185113726, 2085126539, 1315887233, 2008046410] (mod 2147483629).

По китайской теореме об остатках, имеем:

[0, 1, 128, 405688, 24980352, 32989068162, 3101696069920, 2365714170297014, 309656520296472068, 2415277789552788286, 3926649012293853406, 726889843182193849, 153366515247378747, 1645735649663585962, 3698490188721496226, 1337259901989820598] (mod 2147483647 · 2147483629).

Как проверить, хватило ли нам 2 модулей? В этой задаче понятно, что числа должны расти экспоненциально, значит их логарифмы должны расти линейно. Возьмём логарифмы (например, натуральные):

[-INF, 0, 4.852, 12.91, 17.03, 24.22, 28.76, 35.40, 40.27, 42.33, 42.81, 41.13, 39.57, 41.94, 42.75, 41.74].

Логарифм произведения наших двух модулей равен 42.98. Мы видим, что разница между числами монотонно растет примерно на одинаковую величину, а, начиная с n=9, они почему-то перестают расти. Более того, в этом месте был достигнут предел, больше которого числа уже быть не могут. Это означает, что двух модулей не хватило. Давайте возьмем 4 модуля и восстановим ответ:

[0, 1, 128, 405688, 24980352, 32989068162, 3101696069920, 2365714170297014, 309656520296472068, 168435972906750526954, 27738606105535271640888, 12142048779807437697982030, 2344813362310160031818110686, 888511465584607682074513271223, 191678405883294971709423926242394, 65882516522625836326159786165530572] (mod 2147483647 · 2147483629 · 2147483587 · 2147483579).

Логарифмы равны:

[-IFN., 0, 4.852, 12.91, 17.03, 24.22, 28.76, 35.40, 40.27, 46.57, 51.68, 57.76, 63.02, 68.96, 74.33, 80.17].

Числа строго возрастают, увеличиваясь примерно на одинаковую величину. Логарифм произведения 4 выбранных модулей равен 85.95, что заметно больше последнего числа в последовательности. Это означает, что ответ правильный. Можете в этом убедиться, набрав в гугле запрос «65882516522625836326159786165530572».


8. Обсуждения и результаты



Этим алгоритмом мне удалось посчитать число циклов на решётке размером 22×100. Вычисления по одному модулю длились около 30 часов на 32 ядрах кластера. Всего потребовалось около 50 модулей. Памяти требовалось меньше 1 Гб на ядро. В процессе вычислений я собрал кое-какие сведения о количественных характеристиках решения. В табл. 1 ниже представлена данная информация.

m Реальное число разрезов Число слов Моцкина
3 3 4
4 6 9
5 12 21
6 23 51
7 62 127
8 109 323
9 365 835
10 607 2188
11 2355 5798
12 3774 15511
13 16020 41835
14 25188 113634
15 113198 310572
16 176061 853467
17 821923 2356779
18 1270562 6536382
19 6097041 18199284
10 9387784 50852019
21 46013564 142547559
22 70652188 400763223

Таблица 1 – Сравнение чисел Моцкина и реального числа разрезов, необходимых для хранения в памяти

В первой колонке ширина решётки m. Во третьей – количество слов Моцкина длиной m. В средней колонке указано реальное число разрезов, которые оказываются допустимыми (симметрия учтена и не посчитано финальное состояние).

Если угодно посмотреть число гамильтоновых циклов на квадратной решетке P2n×P2n, то ссылка на эту последовательность есть в списке источников.


9. Список использованных источников


  1. М. Гэри, Д. Джонсон. Вычислительные машины и труднорешаемые задачи.
  2. A003763 – Number of Hamiltonian cycles on 2n X 2n square grid of points.
  3. Китайская теорема об остатках.
  4. Числа Каталана, вывод формулы и асимптотики.


Спасибо за внимание.
Tags:
Hubs:
+88
Comments 16
Comments Comments 16

Articles