Как я решал задачу 2025 года. Часть 1

Рис. 1. Пример заполнения квадрата 45x45. Количество квадратиков каждого типа равно стороне квадратика (от 1 до 9).

1-го января из сообщества Незадача дня я узнал про интересные равенства относительно числа 2025 и про задачу, которую на их основе можно сформулировать.

Равенства следующие:

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

На основе этих равенств можно сформулировать задачу:

Сколько существует способов расположить 1 квадратик со стороной 1, 2 квадратика со стороной 2, 3 квадратика со стороной 3, … , 8 квадратиков со стороной 8, 9 квадратиков со стороной 9 в квадрате со стороной 45, чтобы они не пересекались?

Как следует из равенства, суммарные площади будут равны, и остаётся только расположить квадратики на плоскости без пересечений. Один из способов представить эту задачу: каждый куб раскладывается на слои толщины 1 и этими плоскими фигурами нужно заполнить квадрат со стороной 45.

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

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

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

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

Сначала я начал думать в сторону алгоритма, который полностью заполняет одну четверть большого квадрата (таких вариантов не очень много) и потом стыкует полученные раскладки, чтобы в итоге получился полностью заполненный квадрат. Если вкратце, можно заполнять четверть сверху вниз, слева направо, так чтобы центры квадратиков лежали не правее правой границы четверти и выше нижней границы четверти. Отдельно обрабатываются 6 вариантов центрального квадратика. Как выяснилось в процессе реализации, граница такой раскладки может содержать разрывы, поэтому закодировать её как одну линию не удастся. Не считая этой проблемы, идея, вообще, мне очень нравилась: предполагалось хранить границу в виде префиксного дерева, чтобы потом можно было относительно быстро состыковать соседние четверти.

Рис. 2. Заполнение одной четверти квадрата. Полупрозрачные квадратики не включаются в раскладку, так как их центры лежат за пределами четверти. Но их всё равно надо расставлять, чтобы понять, какие области заполнены, а какие – нет.

Рис. 3. Случай, когда граница разрывна. Видно, что квадратик размера 2 не соединён с остальной границей включённых квадратиков. Это происходит, когда квадратик из правой верхней четверти касается квадратика из левой нижней четверти и при этом после них ещё идёт квадратик из левой верхней четверти.

Поняв, что граница может быть разрывная, я попробовал ещё одну вариацию алгоритма: заполнить верхнюю половину квадрата, сохранить все варианты и потом, аналогично предыдущему алгоритму, подобрать раскладку для нижней половины. Как оказалось, варианты заполнения половины квадрата 45x45 занимают слишком много оперативной памяти, даже с учётом того, что на сервере было более 200 Гб памяти. Я имел возможность примерно в 2 раза сократить использование памяти за счёт учёта симметрии, но при этом я ещё не хранил полный список квадратиков, а лишь границу и мультимножество используемых квадратиков. Пришлось бы запускать программу 2 раза, чтобы получить искомые расположения квадратиков.

Рис. 4. Если мы переберём все варианты заполнения верхней половины квадрата и сохраним форму границы, количества каждого квадратика и центральный квадратик, то найти заполнения всего квадрата можно с помощью поиска дополняющей границы. При этом в сумме квадратики должны насчитываться в количестве, определённом в условии задачи. Чтобы половины дополняли друг друга, левая половина центральной линии включается строго, а правая половина - нестрого.

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

Финальная идея была следующая - использовать графический процессор (GPU) для перебора вариантов. Такая идея витала с самого начала, но меня останавливало то, что нужно будет преобразовать рекурсивный алгоритм в итеративный, и что на GPU нет ассоциативных массивов, которые у меня использовались для хранения свободных промежутков. Из плюсов вычислений на GPU: 1. пониженное энергопотребление GPU по сравнению с вычислением на центральном процессоре и 2. большое число параллельных потоков. В итоге, собравшись с духом, я развернул рекурсию в итеративный алгоритм и заменил std::map<> на наивную реализацию std::flat_map<> - хранение упорядоченного списка ключей и соответствующих им значений.

Верхнеуровнево алгоритм выглядит следующим образом:

1. Находим первый незаполненный промежуток сверху вниз, слева направо;

2. Пробуем заполнить его оставшимися квадратиками, и в каждом успешном случае применяем рекурсию и переходим к шагу 1.

Рис. 5. Пример перебора вариантов следующего квадратика.

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

Рис. 6. Пример множества незаполненных промежутков — они отмечены зелёным цветом.

Допустим следующий свободный промежуток находится на уровне Y и продолжается от L до R по оси X.
Мы будем делать поиск в глубину: переберём оставшиеся размеры квадратиков и попробуем разместить каждый из них в положении (Y, L) - см. Рис. 5.
Пусть мы проверяем квадратик размера S.
Для простоты можно рассмотреть 2 случая:
1. S < R-L, то есть квадратик помещается с запасом
В этом случае нужно заменить промежуток (Y, (L, R)) на промежуток (Y, (L+S, R)), а также добавить промежуток (Y+S, (L, R)). При этом есть 2 варианта - на уровне Y+S есть промежуток, который заканчивается на позиции L или нет. Если есть, то его нужно объединить с промежутком (L, R).
2. S = R-L, то есть квадратик помещается "впритык"
В этом случае нужно полностью убрать промежуток (Y, (L, R)) и добавить промежуток (Y+S, (L,R)), который, возможно, объединится с левым и/или правым промежутком на уровне Y+S.

Рис. 7. Пример объединения промежутков на уровне Y+S, когда три промежутка объединяются в один.

Также нужно проверить, что квадратики влезают по высоте, то есть Y+S <= 45.

При откате, разумеется, нужно вернуть структуру данных промежутков в прежний вид.

Для хранения промежутков можно использовать два массива размера 45x12, которые будут хранить начала и концы промежутков на уровне Y в упорядоченном по X виде. Здесь используется 12 в качестве второго измерения, так как в одной строке точно не может возникнуть более 11 промежутков. Для обозначения окончания списка начал и концов можно остаток массивов заполнить значением «-1».

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

Я бы не сказал, что данный алгоритм хорошо ложится на архитектуру вычислений на GPU. Всё упирается в то, что каждое ядро GPU в пределах некоторого количества потоков (обычно 32 потока) выполняет в один момент времени ровно одну строчку кода. Поэтому сложности возникают в связи с тремя аспектами:

1. В один момент времени выполняется либо шаг вперёд рекурсии, либо откат. Соответственно, производительность падает примерно в 1,5-2 раза, так как скорее всего хотя бы один поток захочет сделать откат, в то время как большинству остальных потоков его делать не нужно.

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

3. Некоторые ветки могут выполняться существенно дольше остальных. В этом случае выполнение 32 параллельных задач может выродиться в выполнение 1 задачи, что будет означать 32-кратное сокращение эффективности вычислений.

Последний пункт удалось ускорить с помощью нехитрого алгоритма.

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

Что же сделать, чтобы справиться с неоднородностью размера задач? А давайте будем выполнять определённое количество итераций (например, ), и если за это время задача выполнена полностью, то ОК. А если нет, то добавляем к задаче ещё 1 квадратик (всевозможные его варианты) и тем самым подразбиваем задачу примерно на 9 подзадач. Таким образом задачи будут со временем решаться, а какие не решены - подразбиваться на более мелкие. Заметим, что при разбиении нерешённой задачи какая-то часть вычислений будет повторяться, но зато заполненность ядер графического процессора будет достаточно большой. К слову, я не пробовал менять лимит на число итераций, может быть в этом есть какой-то дополнительный выигрыш по времени.

Рис. 8. Разбиение задачи на подзадачи с ограничением по количеству итераций.

Для такой алгоритм сократил время решения задачи с 1 часа до 15 минут. Интересно наблюдать за алгоритмом: сначала число подзадач растёт (причём не в 9 раз, а обычно не более чем в 2 раза, так как некоторые подзадачи успешно решаются), а потом, начиная с какого-то момента, их число начинает сокращаться, пока все подзадачи в итоге не будут решены.

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

Ну и напоследок расскажу про ещё одну оптимизацию, которая ускорила вычисления более чем в 2 раза.
Несложно показать, что для каждого заполнения квадрата существуют ещё 7 различных вариантов заполнения, которые получаются из него: симметрия относительно вертикали; поворот на 90, 180, 270 градусов; симметрия + поворот на 90, 180, 270 градусов.
Таким образом нам достаточно найти хотя бы одно из этих заполнений, а остальные мы легко сможем построить.

Рис. 9. Заполнение и его 7 вариаций, которые получаются с помощью поворота на 90 градусов против часовой стрелки и симметрии. Повторений не будет из-за ограничений на количество квадратиков каждого типа.

Соответственно, если мы сможем каким-то способом упорядочить эти заполнения, то мы можем искать только те заполнения, которые идут под номером 1 (назовём их "каноническими" заполнениями). Как же это сделать? Оказывается, несложно. Возьмём заполнение и запишем его как последовательность , где , и - координаты и размер квадратиков. Далее, отсортируем эти тройки в лексикографическом порядке, то есть сначала по Y, а потом по X. Несложно заметить, что это будет соответствовать заполнению сверху вниз, слева направо, и, в случае полного заполнения, всю последовательность можно восстановить только по последовательности размеров квадратиков . Соответственно, для сравнения заполнений, можно сравнивать в лексикографическом порядке только их кодировки (как выяснилось, это называется Bouwkampcode или tablecode).

Таким образом, когда у нас есть какой-то набор подзадач, можно из них оставить только те, которые при применении симметрии и поворота являются наименьшими относительно лексикографического сравнения их кодировок. Разумеется, когда у нас заполнен только верх квадрата, нельзя применять повороты на 180 и 270 градусов, так как тогда начало кодировки будет неизвестно. Но как показала практика, сравнение 4-х вариантов ускорило вычисления более чем в 2 раза. (Не в 4 раза, из-за того, что отбор подзадач производился на Питоне в однопоточном режиме.) В конце нужно не забыть отобрать только канонические заполнения, чтобы правильно посчитать общее количество вариантов («каноничность» оценивалась только на уровне подзадач, сам алгоритм перебора при этом остался прежним).

Итог

У меня получилось 216'285 вариантов заполнения квадрата 45x45 квадратиками соответствующих размеров, или в 8 раз больше, если включить в ответ все повороты и симметрии.

В заключение приведу ссылку на свой код на Python Numba CUDA, который вычислил число заполнений квадрата 45x45 менее чем за 4 суток на довольно мощном графическом процессоре Nvidia A100. Объём потребляемой видеопамяти был менее 8 Гб, хотя этот момент можно ещё оптимизировать, заменив int32 на int8 при передаче данных на GPU.

Для того, чтобы убедиться в правильности результата, я посчитал часть заполнений более простым кодом на C++ и проверил, что все они были найдены и ничего более не было найдено (в определённых границах подзадач). Также я использовал следующее свойство заполнений: если в заполнении есть подпрямоугольник, полностью заполненный квадратиками разного размера, то их можно переставить и получить новое заполнение. Соответственно, во всех заполнениях я брал полностью заполненные прямоугольники площади 225, переставлял в них квадратики и убеждался, что полученные заполнения (или их повороты-симметрии) также были найдены. Ограничение на 225 чисто техническое, иначе бы перебор вариантов занял большее время. Например, вообще без ограничения, для некоторых заполнений пришлось бы решать полную задачу для n=8. Проверка прошла успешно, для всех таких пертурбаций вновь полученные заполнения уже были в списке найденных.

Рис. 10. Пример локальных пертурбаций в некоторых заполнениях. Выделенный жёлтым цветом прямоугольник содержит квадратики разных размеров, поэтому их можно переставить в другом порядке.

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

Спасибо за внимание!
Анализ интересных закономерностей заполнений приведён в части 2.