В конце второй статьи я попытался решить еще одну задачу, связанную с ходом коня и подсчитать количество замкнутых маршрутов в прямоугольнике m x n, но дальше квадрата 6x6 не продвинулся. После ряда оптимизаций удалось ускорить вычисления на шесть порядков, т.е. примерно в миллион раз и вплотную приблизиться к квадрату 8x8, вычислив количество циклов в прямоугольнике 7x8.

Пусть квадрат 8x8 по-прежнему кажется недоступным грубому перебору, но такое ускорение говорит о хорошем потенциале и языка и задачи в целом. И, собственно, опытом раскрытия этих потенциалов хотелось бы поделиться с читателями.

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

Одна эта проверка позволяет на два с лишним порядка сократить вычисления. И это при том, что в списковом варианте она выглядит достаточно мудрено и, как и проверка на связность, имеет квадратичную сложность. А сложность эта обусловлена тем, что для каждой клетки, раз за разом ищется почти всегда один и тот же список соседей.

Переходим ко второй идее – для всех клеток вычислить списки связности загодя один раз, и затем лишь изменять по мере надобности. Хранить пары ключ-значение можно по-разному, хоть в тех же списках, но у списков операции поиска и удаления имеют сложность O(n). Есть более оптимальные структуры данных, например сбалансированные бинарные деревья, и в языке Haskell издавна такая структура реализована в модуле Data.Map. Все основные операции с этой структурой имеют сложность не более O(log n), а, главное, функция поиска изолированных и висячих вершин теперь принимает элегантный вид

deadlocks = keys . filter ((<2).length)

И ее сложность улучшается до O(n).

Интерфейсная функция из-за объема подготовительных действий достаточно сильно преображается


kNCircles m n = 
    kNFromTo [(2,3)] (3,2) $ prepare $ 
    tail [(x,y) | x <- [1..m], y <- [1..n]]
    where
    prepare xs = M.fromList 
        [(x, filter (near x) xs) | x <- xs]
    near (x1,y1) (x2,y2) = 
        abs ((x2-x1)*(y2-y1)) == 2

Что же до основной рекурсии, то раз уж сами маршруты теперь неинтересны, нас интересует только их количество, от выстраивания цепочек можно перейти к простому подсчету числа успехов. И с учетом новой структуры данных функция принимает следующий вид


kNFromTo ks s xs
    | size xs == 1 = 1
    | otherwise = sum
        [kNFromTo (xs ! k) s (kDel k xs) | 
        k <- filter (/= s) ks,
        null $ deadlocks xs \\ [k,s]]

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


kDel x xs = delete x $ foldr 
    (adjust (delete x)) xs (xs ! x)

Выглядит наворочено, но суммарная сложность осталась O(log n), пусть и с достаточно большим коэффициентом

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

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

deadlocks xs = filter ((<2).length.(xs !))

И тем самым сложность проверки улучшить до O(log n).

Все, слабых мест не осталось, можно лишь слегка подкручивать то, что есть. Точнее, теоретически, можно было бы еще ускориться в пару-тройку раз, перейдя к массивам. Собственно массивы для прямоугольников вообще самая очевидная структура данных. И поиск элемента требует O(1) времени, а это самая востребованная операция. Но в Хаскеле для сохранения чистоты при изменении элементов создается новая копия массива, а это уже O(n). В результате вместо ускорения получается в полтора раза медленней.

Но! До сих пор мы копали проблему вглубь, а можно пойти и, так сказать, «вширь». Все вычисления в языке Haskell по умолчанию выполняются однопоточно, на одном ядре. А ядер, пусть и, порой, виртуальных, бывает несколько. И раз уж пошла борьба за производительнось, загрузить их хотелось бы все, это позволит еще на порядок (плюс-минус) ускорить вычисления.

В поисках материала по распараллеливанию вычислений, я наткнулся на интересную статью, в которой рассказывается, что чистые вычисления в языке Haskell распараллеливаются не просто, а очень просто. И у нас для этого уже почти все есть. Достаточно подгрузить модуль Control.Parallel.Strategies и в рекурсии после конструктора списков добавить магическую строчку `using` parList rdeepseq.

Ниже я привел листинг окончательного варианта. Программу необходимо скомпилировать с ключом -threaded и при запуске в параметрах кроме размеров прямоугольника указать ключи +RTS -N, например knc 7 6 +RTS -N.

knc.hs


import Data.List(delete, (\\))
import qualified Data.Map.Lazy as M
import Control.Parallel.Strategies
import System.Environment

type Cell = (Int, Int)
type Pool = M.Map Cell [Cell]

kDel :: Cell -> Pool -> Pool
kDel x xs = M.delete x $ foldr 
    (M.adjust (delete x)) xs (xs M.! x)
    
deadlocks :: Pool -> [Cell] -> [Cell]
deadlocks xs = filter ((<2).length.(xs M.!)) 

kNFromTo :: [Cell] -> Cell -> Pool -> Int
kNFromTo ks s xs
    | M.size xs == 1 = 1
    | otherwise = sum ( 
        [kNFromTo (xs M.! k) s (kDel k xs) | 
        k <- ks, k /= s,
        null $ deadlocks xs ks \\ [k,s]]
        `using` parList rdeepseq )

kNCircles :: Int -> Int -> Int
kNCircles m n = 
    kNFromTo [(2,3)] (3,2) $ prepare $ 
    tail [(x,y) | x <- [1..m], y <- [1..n]]
    where
    prepare xs = M.fromList 
        [(x, filter (near x) xs) | x <- xs]
    near (x1,y1) (x2,y2) = 
        abs ((x2-x1)*(y2-y1)) == 2

main = do
    [m,n] <- getArgs
    print $ kNCircles (read m) (read n)

Простота реализации параллелизма и его результат, конечно, впечатляют. Теоретически программа способна распараллелиться на число потоков, равное числу ветвлений, т.е. общему числу замкнутых маршрутов. Такого количества ядер пока, к сожалению, нет, но имеющиеся 8 (пробовал и 24) загружаются на 100% с кратным ускорением. И это финальное ускорение позволило за неделю вычислений расколоть прямоугольник 7x8, найдя в нем 34 524 432 316 циклов. Это оказалось даже больше ожидаемого, и теперь оценка из википедии для квадрата 8x8 видится вполне реальной.

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

3x2n:

0, 0, 0, 0, 16, 176, 1536, 15424, 147728, 1448416, 14060048, 136947616, 1332257856, 12965578752, …

5x2n:

0, 0, 8, 44202, 13311268, 4557702762, …

6xn:

0, 0, 0, 0, 8, 9862, 1067638, 55488142, 3374967940, 239187240144, …

7x2n:

0, 0, 1067638, 34524432316, …

И хоть я и написал в заголовке слово «заключение», в самой задаче точку ставить рано.

P.S. Задача таки сразу не отпустила...

… периодически приходили идеи по дальнейшей оптимизации вычислений. Срабатывало, разумеется, не все, но тонким тюнингом удалось еще в три раза ускориться. Не очень много, перебор квадрата 8х8 на одной машине теперь требует три года работы, а вот на ферме из двух десятков узлов результата удалось достичь за два с небольшим месяца. И этот результат полностью совпал с Википедией, квадрат 8х8 содержит 13 267 364 410 532 замкнутых маршрута. Сенсации не произошло, на пять лет опоздал, но теперь, надеюсь, отпустит )


import Data.List (delete)
import qualified Data.Map.Lazy as M
import Control.Parallel.Strategies
import System.Environment (getArgs)

type Cell = Int
type Pool = M.Map Cell [Cell]

kDel :: Cell -> Pool -> Pool
kDel x xs = M.delete x $ foldr 
    (M.adjust (delete x)) xs (xs M.! x)

kNC :: [Cell] -> Pool -> Integer
kNC ks xs
    | M.size xs == 4 = 1
    | otherwise = 
    let ds = filter (null.tail.(xs M.!)) ks
    in  if null ds
        then sum ( 
            [ kNC (xs M.! k) (kDel k xs) |
            k <- ks,  k /= 1 ]
            `using` parList rseq )
        else let k = head ds in
            if null (tail ds) && k /= 1
            then kNC (xs M.! k) (kDel k xs)
        else 0

kNCircles :: Int -> Int -> Integer
kNCircles m n = 
    kNC [m] $ prepare $ tail
        [(x,y) | x <- [1..m], y <- [1..n]]
    where
    prepare xs = 
        M.adjust (++[1]) 1 $ M.fromList 
        [(enc x, enc <$> filter (near x) xs) |
        x <- xs]
    near (x1,y1) (x2,y2) = 
        abs ((x2 - x1) * (y2 - y1)) == 2
    enc (x,y) = (y - 2) * m + x - 2

main = do
    [m,n] <- getArgs
    print $ kNCircles (read m) (read n)

→ Часть первая
→ Часть вторая

Хаскель — ход конем 3. Заключение