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

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

где— множество всех параллелепипедов размера, где в каждой ячейке стоит 0 или 1. Итогослагаемых. Эти параллелепипеды мы назвали шляпами.

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

Таким образом, если мы объединим все шляпы в так называемые классы эквивалентности, где в каждом классе будут находиться одинаковые шляпы с точностью до перестановки слоёв, мы сможем вычислить ответ быстрее. Забегая вперёд — эти классы эквивалентности на языке комбинаторной теории групп называются орбитами. Тогда формула примет вид

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

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

Азы теории групп

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

1. ассоциативность:для любых,

2. наличие нейтрального элемента, такого, чтодля любого,

3. наличие обратного элементадля любого , при этом
   .

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

Для произвольного подмножества элементови произвольного элементаопределим множествокак. Аналогично можно определить.

— мощность множества, то есть количество элементов в. Далее мы будем работать только с конечными группами, то естьвсегда будет равно какому-то неотрицательному целому числу. Если— группа, тотакже называется порядком группы.

Лемма 1. Для любогои любогоимеет место.

Доказательство. Каждому элементупоставим в соответствие элемент, а каждому элементу поставим в соответствие элемент. Таким образом, получаем взаимно однозначное соответствие (биекцию) между элементамии, откуда. Аналогично можно доказать, что.

Смысл Леммы 1 в том, что мы можем прокручивать произвольные подмножествапо группе, и в процессе у нас ничего не слипнется. Мы также всегда можем открутитьобратно.

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

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

Лемма 2. Смежные классыпо подгруппелибо совпадают, либо не пересекаются.

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

Для подгруппыгруппы каждый элементпринадлежит ровно одному смежному классу, а именно —(двум смежным классампринадлежать не может по Лемме 2). Стало быть, все элементыразбиваются на непересекающиеся смежные классы. По Лемме 1 размер каждого смежного классаравен. Отсюда прямиком следует

Теорема Лагранжа. В конечных группах порядок любой подгруппы делит порядок группы.

Действие группы на множестве

Группы и их элементы зачастую довольно абстрактны в том смысле, что это — какие-то абстрактные действия в "вакууме". Чтобы получить наглядные примеры и больше пользы от происходящего — эти действия нужно к чему-то применять, к каким-то объектам.

Обозначим множество объектов, на которые мы будем применять нашу группу (с бинарной операцией), как. Мы также должны ввести операцию действия группы (слева!), которая каждой паре элемента группыи элемента множествабудет ставить в соответствие какой-то элемент из. Кроме, того, чтобы всё работало, нужно выполнение дополнительных аксиом:

1. ассоциативность: для любыхи любоговерно,

2. нейтральный элемент для любогоничего не меняет:.

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

Орбитой элемента, на который действует группа, называется множество.— количество элементов в орбите называется длиной орбиты.

Лемма 3. Две орбиты либо совпадают, либо не пересекаются.

Доказательство. Аналогично тому, что мы делали со смежными классами: пусть элементпринадлежит двум орбитами. Тогда существуют такие, что. Тогда. Для любогоимеем. Стало быть,. Аналогично доказываем, что, откуда следует, что.

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

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

Заможно взять множество всех вещественных чисел.— это значение линейной функциив точке. Орбит бесконечно много, поскольку всчётное число элементов, а в— континуум.

Пример 2. Пусть— множество невырожденных вещественных матриц размера, а— множество всех вещественных векторов длины. Операциии— это умножение матриц и умножение матрицы на вектор соответственно. Орбиты две. В одной из них нулевой вектор, в другой — все остальные векторы. В похожей модели (для) крутятся, сдвигаются и масштабируются трёхмерные модели в компьютерной графике.

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

Пример 4.— всевозможные повороты плоскости относительно начала координат, а— все точки этой плоскости. Орбитами будут окружности всевозможных радиусов с центром в начале координат. Ну и еще начало координат в своей личной орбите.

Пример 5. Предыдущие примеры были не совсем корректные в рамках данной статьи в том плане, чтоибыли бесконечными. Мы все таки хотим работать с конечными объектами и симметриями в множестве. Для этого модифицируем Пример 2. Возьмём втолько такие матрицы, в которых в каждой строке и каждом столбце ровно одна единица, а в остальных ячейках — нули (это, кстати, подгруппа группы всех невырожденных матриц). Ввозьмём только 0/1-вектора длины. Тогда оба множества будут конечного порядка. Нетрудно заметить, чтобудет изоморфно множеству перестановокдлины, и, действуя элементами этой группы на элементы из, мы будем каким-то образом перемешивать числа внутри соответствующих векторов. Количество орбит равно, в одной орбите все вектора с одним и тем же количеством единиц. В-й орбите будет ровновекторов.

Пример 6 совсем игрушечный, с которым мы будем работать дальше. Пусть— все 16 таблиц, каждая из ячеек которой покрашена либо в белый, либо в чёрный цвет. Будем считать две таблицы одинаковыми, если одну из них можно перевести в другую перестановками строк и столбцов. Это последнее правило позволяет нам построить группу преобразований, в которой будет 4 элемента:

1. Нейтральный элемент

2. Обмен строк

3. Обмен столбцов

4. Композиция элементов 2 и 3, когда мы обмениваем ячейки таблицы крест-накрест

Результат действия каждого элемента группына каждый элемент множестваможно показать в виде следующей таблицы:

Аккуратно проанализировав все 16 раскрасок, можно разбить их все на следующие 7 орбит:

Лемма Бернсайда и теорема Пойа

Стабилизатором элемента, на который действует группа, называется подмножество элементов, которые не меняют:.

Множеством неподвижных точек для элементав множественазывается подмножество элементов, которые не меняются под действием:.

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

Лемма 4. Имеет место тождество

Доказательство. Левая часть равенства — это красные кружочки, просуммированные по строкам, а правая часть — они же, просуммированные по столбцам.

Лемма 5. В рамках действия группыдля любого элементаверно

Доказательство. Сначала заметим, что— это подгруппа. В этом легко удостовериться, проверив все аксиомы. Значит, по теореме Лагранжа группараспадается на смежные классы по подгруппе. Выпишем эти классы:

где,, а. Действие каждого-го из этих смежных классов на элементдаёт множество из единственного элемента. Действительно, для любогоимеем.

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

Осталось осознать, что— это, собственно, различные элементы орбиты, то есть, а, как указано выше.

Лемма Бернсайда. Пусть— множество орбит, тогда

Доказательство. В одну строчку:

Теперь пройдёмся по формуле. Сначала бы проходимся по всем орбитам и всем элементам в орбите назначаем вес, где. Это эквивалентно просто суммированию по всем элементам. Далее мы последовательно применяем сначала Лемму 5, а затем Лемму 4.

Зачем вообще нужна лемма Бернсайда? С её помощью мы можем быстрее посчитать количество орбит когда порядок группысильно меньше, чем размер множества. Нужно только уметь считать количество неподвижных точек для каждого действия.

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

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

Теорема Пойа. Пусть— количество циклов в перестановке. Тогда

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

Для примера применим Теорему Пойа к нашему игрушечному примеру. У нас вчетыре перестановки, в одной из них 4 цикла, в трёх остальных — 2. Табличкимы красим в 2 цвета. Применяем формулу:

Сошлось!

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

Количество орбит шляп

Теперь можно применить теорему Пойа к нашей задачи имени Винтика и Шпунтика. У нас объекты в— это параллелепипеды размера, каждая ячейка которых покрашена в один из двух цветов, и у настаких параллелепипедов. Ву нас всевозможные действия, которые перемешивают слои параллелепипедов., так как мы можем независимо перемешивать слои по каждой из трёх осей координат, ну и, в принципе, у насизоморфна прямому произведению полных групп перестановок. Следующий код вычисляет количество орбит по формуле:

import itertools
import math

def get_num_of_classes( N, M, K ):
	Gsz = math.factorial(N) * math.factorial(M) * math.factorial(K)
	sum = 0
	for pi in itertools.permutations(range(N)):
		for pj in itertools.permutations(range(M)):
			for pk in itertools.permutations(range(K)):
				F = [False] * (N*M*K)
				cycles = 0
				for i in range(N):
					for j in range(M):
						for k in range(K):
							if F[i*M*K+j*K+k]: continue
							cycles += 1
							ii,jj,kk = i,j,k
							while not F[ii*M*K+jj*K+kk]:
								F[ii*M*K+jj*K+kk] = True
								ii,jj,kk = pi[ii],pj[jj],pk[kk]
				sum += 2 ** cycles
	print( (N,M,K), 2 ** (N*M*K), sum // Gsz )
    return sum // Gsz

Запустим код и получим следующее для наших 14 задач:

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

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

И вот, например, для задачипрограмма из предыдущей части считала бы ответ около 24 миллиардов лет. Теперь мы можем посчитать его "всего" за 2 миллиона лет. Прогресс налицо. Мы можем еще добавить в группу действия, которые вращают шляпу(ну или, что то же самое, меняют оси координат местами). Ускорение в 6 раз, итого — 300 тысяч лет. Но это программа на Python, если её переписать на C, хардкорно оптимизировать, а затем запустить параллельные вычисления на десятке видеокарт, то можно уложиться лет за 5. То есть, в принципе, задача уже решаемая.

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

Генерация орбит

Как генерировать орбиты? Один из самых простых способов следующий. Можно взять все элементыи представить их как вершины графа. При помощи действий издля любогоможно находить соседей, и таким образом проводить между вершинами рёбра. Если группабольшая, то можно найти в ней генераторы или порождающие элементы. Это такое подмножество действий, через которые можно выразить все остальные действия. Когда граф построен — нужно выделить в нём компоненты связности, это и будут орбиты. Минус такого подхода в том, что нужно очень много памяти для хранения всех посещённых вершин.

Другая идея состоит в том, чтобы выделить из каждой орбиты ровно одного представителя. Обычно это делается так: если объекты изсостоят из"штучек", каждая из которых красится в один изцветов, то можно на этих"штучках" задать порядок, а цвета пронумеровать числами, скажем, отдо. Тогда каждому объекту изможно поставить в соответствие-значное число в-ичной системе счисления (или строку длиныс алфавитом изсимволов). Представителем орбиты будет объект с максимальным числом (или с лексикографически максимальной строкой) по всем объектам в орбите. Будем называть такой объект максимальным. Теперь можно просто перебрать все объекты из, проверить максимальны ли они и отсеять все объекты, которые не максимальны. Таким образом можно получить ровно по одному представителю из каждой орбиты.

Как проверить, что объект максимальный? Нужно попробовать применить к нему все действия из группы. Если после какого-то действия в результате получится объект больше — то объект не максимальный. Если увеличения объекта не произошло для всех действий группы — он максимальный.

Посчитать количество объектов в орбите без затрат памяти тоже на самом деле не сложно. Нужно воспользоваться Леммой 5 и посчитать, то есть количество действий группы, которые переводят объект сам в себя. Тогда длина орбитыбудет равна.

Таким же способом можно просто искать максимальный объект в орбите. Возьмём какой-то объект, применим к нему действия группы всеми способами и возьмём максимум из того, что получилось. Например, если группа — это множество всех перестановок, то заможно найти максимальный объект в орбите, в котором все "штучки" упорядочены в порядке невозрастания. Ой, это же BogoSort. Человечеству известны чуть более быстрые алгоритмы сортировки, так что на самом деле в некоторых случаях искать максимальный объект в орбите можно быстрее. Но об этом поговорим в следующем параграфе.

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

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

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

Лемма 6. Если частично покрашенный объект не является максимальным, то не являются максимальными также все объекты в его поддереве.

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

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

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

Согласно Лемме 6, если в результате рекурсивного обхода дерева мы наткнёмся на не максимальную вершину, то нам следует откатиться, ибо во всём поддереве нам ловить нечего. Напротив, если частично покрашенный объект максимален — в его поддереве существует максимальный лист. Например, объект, в котором все не прокрашенные "штучки" покрашены в минимальный цвет.

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

Алгоритм можно дополнительно ускорить, используя секретную технику шаолиньских монахов, которую я называю халтурная проверка максимальности. Будем тестировать объекты на максимальность не всеми действиями группы, а только некоторыми, выбранными наугад. Тогда мы сэкономим время на проверку объектов, которые действительно являются максимальными (напомню, что на них надо испробовать все действия группы, чтобы убедиться в максимальности). Не максимальные объекты будут отсеиваться лишь с некоторой вероятностью, и мы будем вынуждены проверять некоторые дополнительные "лишние" не максимальные объекты. Но принцип здесь такой: пусть мы будем проверять в 10 раз больше объектов, зато каждая проверка будет в 100 раз быстрее, тогда в итоге получится ускорение в 10 раз. Понятно, что листы дерева поиска надо будет проверять полностью без всякой халтуры.

Работу алгоритма можно проиллюстрировать следующей картинкой для нашего игрушечного примера:

Как можно видеть, в этом дереве 7 листов, что совпадает с числом орбит.

Эксплуатируем структуру группы

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

Чуть выше уже упоминалось, чтоу нас изоморфна. То есть каждое действие можно представить как последовательность из трёх действий: сначала перемешаем слои шляпы, перпендикулярные оси(способов), потом перпендикулярные оси(тутспособов), и, наконец, перпендикулярные оси(способов —).

Вот именно на последней перестановке мы и сэкономим! А сэкономим как и описано выше — просто отсортируем уже перемешанные двумя первыми перестановками слои в порядке убывания. Порядок битов для лексикографического сравнения должен быть соответствующий: из более раннего по порядку слоя биты должны быть более значимы, чем из более позднего. Также биты должны быть одинакового упорядочены по значимости в каждом слое.

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

def is_maximal( N, M, K, hat ):
	for pi in itertools.permutations(range(N)):
		for pj in itertools.permutations(range(M)):
			hat2 = [0] * K
			for k in range(K):
				for i in pi:
					for j in pj:
						hat2[k] = hat2[k] * 2 + ((hat[k]>>(i*M+j))&1)
			hat2.sort( reverse=True )
			if hat2 > hat: return False
	return True

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

Подсчёт порядка стабилизатора тоже можно ускорить до. Действительно, пусть мы как-то перемешали слои по первым двум осям координат. Тогда по третьей оси координат нам опять нужно просто отсортировать все слои (мы предполагаем, что считаем порядок стабилизатора для максимального объекта). Если после сортировки получился другой список — тут элементов стабилизатора нет. Если списки совпали — то хотя бы один элемент стабилизатора есть, но их может быть больше одного. Чтобы найти их количество, нужно разбить все слои шляпы на группы одинаковых, взять факториалы размеров полученных групп, а потом всё перемножить. Код (сразу для вычисления длины орбиты):

def get_size_of_class( N, M, K, hat ):
	G_sz = math.factorial(N) * math.factorial(M) * math.factorial(K)
	Stab_sz = 0
	for pi in itertools.permutations(range(N)):
		for pj in itertools.permutations(range(M)):
			hat2 = [0] * K
			for k in range(K):
				for i in pi:
					for j in pj:
						hat2[k] = hat2[k] * 2 + ((hat[k]>>(i*M+j))&1)
			hat2.sort( reverse=True )
			if hat == hat2:
				ways = 1
				cnt = 1
				for i in range(1,K+1):
					if i==K or hat[i]!=hat[i-1]:
						ways *= math.factorial(cnt)
						cnt = 1
					else: cnt += 1
				Stab_sz += ways
	return G_sz // Stab_sz

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

Остаток кода и результаты запуска

Собственно, код ниже. Недостающие функции можно выковырять из текста статьи выше.

import sys
import itertools
import math
from timeit import default_timer as timer

def sign( x, L ):
	if (L-x.bit_count())&1:
		return -1
	return 1

def calc_dp2t( M, n, m ):
	trits = []
	for i in range(n*m):
		trits.append( M%3 )
		M //= 3
	res = 0
	for x in range(1<<n):
		for y in range(1<<m):
			prod = 1
			for i in range(n):
				for j in range(m):
					prod *= trits[i*m+j] + ((x>>i)&1) + ((y>>j)&1)
			res += sign(x, n) * sign(y, m) * prod
	return res

to_ternary = []
dp = []
max_mask = []
found_classes = 0
num_of_classes = 0
ans = 0

def init( N, M, K ):
	global to_ternary, dp, max_mask, found_classes, num_of_classes, ans
	# init to_ternary
	for i in range(1<<(N*M)):
		tmp = 0
		for j in range(N*M):
			if ((i>>j)&1):
				tmp += 3 ** j
		to_ternary.append( tmp )
	# init dp
	dp = [ calc_dp2t( i, N, M ) for i in range(3**(N*M)) ]
	# init max_mask
	for mask in range(1<<(N*M)):
		mx = mask
		for pi in itertools.permutations(range(N)):
			for pj in itertools.permutations(range(M)):
				mask2 = 0
				for i in pi:
					for j in pj:
						mask2 = mask2 * 2 + ((mask>>(i*M+j))&1)
				mx = max( mx, mask2 )
		max_mask.append( mx )
	# init num_of_classes
	found_classes = 0
	num_of_classes = get_num_of_classes( N, M, K )
	ans = 0

def calcF( N, M, K, T, hat ): # hat is list of K masks
	global to_ternary, dp
	sum = 0
	for apron in range(1<<(N*M)):
		prod = 1
		for part in range(K):
			prod *= dp[to_ternary[hat[part]]+to_ternary[apron]]
		sum += sign(apron, N*M) * prod
	res = sum ** T
	for i in range(K): res *= sign( hat[i], N*M )
	return res

def dfs( N, M, K, T, i, hat ):
	global max_mask, found_classes, num_of_classes, ans
	if i==K:
		found_classes += 1
		if found_classes % 10000 == 0:
			print( found_classes, "/", num_of_classes, ans, file=sys.stderr )
		ans += calcF( N, M, K, T, hat ) * get_size_of_class( N, M, K, hat )
		return
	for part in range(1<<(N*M)):
		if i>0:
  			if part > hat[i-1]: break				# line 82
			if max_mask[part] > hat[0]: continue	# line 83
		hat[i] = part
		if is_maximal( N, M, K, hat ):
			dfs( N, M, K, T, i+1, hat )
		hat[i] = 0

def solve4d( N, M, K, T ):
	global ans
	start = timer()
	print( (N, M, K, T) )
	init( N, M, K )
	dfs( N, M, K, T, 0, [0] * K )
	print( ans )
	end = timer()
	print( "time", end - start )

n, m, k, t = map(int, input().split())
solve4d( n, m, k, t )

Комментарии по коду:

В процедуре init инициализируется всё что можно — таблица dp, как в предыдущей части, считается число орбит через теорему Пойа, а самое главное — вычисляется список max_mask. В нём для каждой возможной маски-слоя хранится максимальная маска, которая может быть получена путём перестановки строк и столбцов. Эту предпросчитанную таблицу вы используем далее для ускорения проверки на максимальность.

Рекурсивная процедура dfs строит дерево перебора шляп. Для листа дерева (когда шляпа полностью сгенерирована) запускается вычисление функции(в коде это calcF), а также вычисляется размер класса эквивалентности (длина орбиты) для данной шляпы. Шляпы же в процессе рекурсивного обходя дерева строятся не побитово, а послойно. Также включены две дополнительные проверки, чтобы уменьшить количество запусков функции is_maximal. Первая из них (строка 82) заключается в том, что маски в слоях шляпы должны идти в порядке невозрастания — иначе мы можем поменять местами слои и получить шляпу больше. Вторая проверка (строка 83) использует посчитанную ранее таблицу max_mask: если маску в текущем слое можно преобразовать в другую, которая больше маски в самом первом слое шляпы, то можно и всю шляпу "перетасовать" этим же преобразованием, а затем перенести маску текущего слоя в самый первый слой. И тогда получится шляпа больше.

В результате этим кодом удалось посчитать ответ для задачиза 3 часа (мы уже решили эту задачу в прошлой части, но ценой 1 месяца вычислений; ответы совпали), задачиза 3.3 часа,за 9 часов,за 29.5 часов, а такжеза 256 часов (почти 11 дней). Какой ответ получился — в сводной таблице ниже.

Примечание. Уже в процессе написания статьи, я нашел в коде баг, где вычислялся список max_mask (тут уже исправленная версия). Этот баг не влияет на ответ, но может повлиять на время вычислений (мне было лень перезапускать расчеты). Из-за него по сути отсечение в строке 83 почти не работало, что приводило к "лишним" запускам полной проверки максимальности. Так что, возможно, на самом деле код выше считает ответ несколько быстрее.

Замечательная симметрия задач (2, *, *, *)

В ходе исследования получаемых орбит (в целях уменьшения их числа) был замечен и доказан следующий факт:

Лемма 7. Если в шляпе размераесть столбик вида, то при изменении его назначение функциине меняется.

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

Я только заменил вот этиина аналоги без штрихов. Этасоответствует некоторому срезуодного из слоёв четырёхмерного гиперпараллелепипеда. То есть— это маска избитов, а— маска из 2 битов. Преобразуем формулу, разрезав этот срез на столбикивсё тем же трюком, которым во второй части сворачивали всю формулу. Получим:

гдеи— это кусочки масокидля-го столбика, а— один из двух вариантов одного бита. Теперь вспомним, что такое. Мы для всех кубиков внутри столбика вычисляем сумму битов, которые на этот кубик смотрят, а затем перемножаем по всем кубикам. То есть вся функцияв своём фундаменте зависит от выражений внутри скобок, которые расписываются как

гдe— это биты, которые "пришли" в-й кубик. Используя заменуи, получим

то есть значениене зависит от порядкаи(как и от порядка других пар битов).

Таким образом, мы даже доказали чуть больше: если поменять в шляпе местами два бита в любом столбике высоты 2, то значениене изменится. И... размер группы у нас возрастает сдо. С соответствующим уменьшением числа орбит. Конечно же, с такой большой группой мы работать не будем, вместо этого рассмотрим троичные шляпы размера. Каждому 0 в такой троичной шляпе мы поставим в соответствие столбикв обычной шляпе, каждой 1 — либо, либо, а каждой 2 — столбик. Каждой троичной шляпе будет соответствовать ровнообычных шляп, где— количество единиц в троичной шляпе.

Я на скорую руку переделал код под задачи вида:

import sys
import itertools
import math
from timeit import default_timer as timer

def sign( x, L ):
	if (L-x.bit_count())&1:
		return -1
	return 1

def calc_dp2t( M, n, m ):
	trits = []
	for i in range(n*m):
		trits.append( M%3 )
		M //= 3
	res = 0
	for x in range(1<<n):
		for y in range(1<<m):
			prod = 1
			for i in range(n):
				for j in range(m):
					prod *= trits[i*m+j] + ((x>>i)&1) + ((y>>j)&1)
			res += sign(x, n) * sign(y, m) * prod
	return res

def get_num_of_classes( N, M, K ):
	Gsz = math.factorial(M) * math.factorial(K)
	sum = 0
	for pj in itertools.permutations(range(M)):
		for pk in itertools.permutations(range(K)):
			F = [False] * (M*K)
			cycles = 0
			for j in range(M):
				for k in range(K):
					if F[j*K+k]: continue
					cycles += 1
					jj,kk = j,k
					while not F[jj*K+kk]:
						F[jj*K+kk] = True
						jj,kk = pj[jj],pk[kk]
			sum += 3 ** cycles
	return sum // Gsz

to_ternary = []
dp = []
max_mask = []
found_classes = 0
num_of_classes = 0
ans = 0
T2 = []

def init( N, M, K ):
	global to_ternary, dp, max_mask, found_classes, num_of_classes, ans, T2
	# init to_ternary
	for i in range(1<<(N*M)):
		tmp = 0
		for j in range(N*M):
			if ((i>>j)&1):
				tmp += 3 ** j
		to_ternary.append( tmp )
	# init dp
	dp = [ calc_dp2t( i, N, M ) for i in range(3**(N*M)) ]
	# init max_mask
	for mask in range(1<<(N*M)):
		mx = mask
		for pi in itertools.permutations(range(N)):
			for pj in itertools.permutations(range(M)):
				mask2 = 0
				for i in pi:
					for j in pj:
						mask2 = mask2 * 2 + ((mask>>(i*M+j))&1)
				mx = max( mx, mask2 )
		max_mask.append( mx )
	# init num_of_classes
	found_classes = 0
	num_of_classes = get_num_of_classes( N, M, K )
	ans = 0
	# init T2 - layers with right ordered columns
	for mask in range(1<<(N*M)):
		flag = True
		for i in range(M):
			if ((mask>>i)&1)==1 and ((mask>>(i+M))&1)==0:
				flag = False
		if flag: T2.append( mask )

def calcF( N, M, K, T, hat ): # hat is list of C masks
	global to_ternary, dp
	sum = 0
	for apron in range(1<<(N*M)):
		prod = 1
		for part in range(K):
			prod *= dp[to_ternary[hat[part]]+to_ternary[apron]]
		sum += sign(apron, N*M) * prod
	res = sum ** T
	for i in range(K): res *= sign( hat[i], N*M )
	return res

def is_maximal( N, M, K, hat ):
	for pj in itertools.permutations(range(M)):
		hat2 = [0] * K
		for k in range(K):
			for i in range(N-1,-1,-1):
				for j in pj:
					hat2[k] = hat2[k] * 2 + ((hat[k]>>(i*M+j))&1)
		hat2.sort( reverse=True )
		if hat2 > hat: return False
	return True

def get_size_of_class( N, M, K, hat ):
	G_sz = math.factorial(M) * math.factorial(K)
	Stab_sz = 0
	for pj in itertools.permutations(range(M)):
		hat2 = [0] * K
		for k in range(K):
			for i in range(N-1,-1,-1):
				for j in pj:
					hat2[k] = hat2[k] * 2 + ((hat[k]>>(i*M+j))&1)
		hat2.sort( reverse=True )
		if hat == hat2:
			ways = 1
			cnt = 1
			for i in range(1,K+1):
				if i==K or hat[i]!=hat[i-1]:
					ways *= math.factorial(cnt)
					cnt = 1
				else: cnt += 1
			Stab_sz += ways
	pw = 0
	for part in hat:
		for i in range(M):
			if ((part>>i)&1)==0 and ((part>>(i+M))&1)==1:
				pw += 1
	return G_sz // Stab_sz * (2 ** pw)

def dfs( N, M, K, T, i, hat ):
	global max_mask, found_classes, num_of_classes, ans, T2
	if i==K:
		found_classes += 1
		if found_classes % 10000 == 0:
			print( found_classes, "/", num_of_classes, ans, file=sys.stderr )
		ans += calcF( N, M, K, T, hat ) * get_size_of_class( N, M, K, hat )
		return
	for part in T2:
		if i>0:
			if part > hat[i-1]: break
			if max_mask[part] > hat[0]: continue
		hat[i] = part
		if is_maximal( N, M, K, hat ):
			dfs( N, M, K, T, i+1, hat )
		hat[i] = 0

def solve4d( N, M, K, T ):
	global ans
	start = timer()
	print( (N, M, K, T) )
	init( N, M, K )
	dfs( N, M, K, T, 0, [0] * K )
	print( ans )
	end = timer()
	print( "time", end - start )

n, m, k, t = map(int, input().split())
solve4d( n, m, k, t )

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

К сожалению, для столбиков высоты 3 и выше эта идея не работает. Там получается

и вот эти вот произведения портят всю малину, и неясно что с ними делать.

Сводные таблички

Сначала таблица по количеству орбит, которая частично повторяет таблицу выше.

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

Теперь таблица с примерными сложностями задач и временем их решения:

Алгоритм с пометкой 1 — это алгоритм из предыдущей части, 2 — из текущей "универсальный", 3 — ускоренный для задач вида. Там, где звёздочка — оценочное время работы, эти задачи до конца не досчитались.

Как можно видеть, удалось решить ещё 5 задач, итого 10 из 14. Решение задачиускорилось с гипотетических 8 месяцев до одной минуты! Для других задач время можно оценить следующим способом: +1 к сложности увеличивает время счёта примерно в 2 раза, +10 к сложности — примерно в 1000 раз.

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

(2, 3, 7, 42)
45547211178376523606515359626169685427015437406338182347055767915998052221756516499184137717124855161194240736076929171622734146809635451918700496379367600348035629119141189588086855391729487239841351995490304000000000000000000000000000000000000000000000000000
(2, 3, 8, 24)
121586269643910956014028765107640267608939361023752056612894734379680384927084682817564023578571488626792890695921312297244842858226123380138967040000000000000000000000000000000000
(2, 3, 9, 18)
146596000834627443441076879480981326388518544901758348791578924651404670854667785842641219589019541525551079075385660136017100800000000000000000000000000000
(2, 3, 10, 15)
35566187691918270224227521123489566281757071652231902417839804745222054702486347418851606130239237409230079000785715200000000000000000000000000000
(2, 3, 12, 12)
538056340148878809909449179046766494467427176926530224806464998699344223285873662921531653464431895203684510288838656000000000000000000000000
(2, 4, 5, 20)
56090549803543632753658611845075457575564534394506257537192124254535376215201630785562247456494950039175717220420431894762857577840640000
(2, 4, 6, 12)
893969652710747906227960820785684673171795118709743136542188903012983716613196946220438854482723340288000
(2, 4, 8, 8)
544127608019523573372789992884706394616053807958652922862285917121018444715169215342633287680000
(2, 5, 5, 10)
25751380617252272796699701337640873547573107139149852088285440091020314414320091036205763788800
(3, 3, 4, 12)
1625717055687087185033963109305255897277923694223134108704431016700065145660728200396800

Порядок чисел, кстати, совпадает с оценками, полученными во второй части. Тем не менее, было бы здорово, если бы кто-нибудь независимо с нуля написал бы программу, аналогичную моей, и верифицировал бы ответы. А то мало ли где я баг мог посадить.

Заключение

В общем, мы решили 10 задачек из 14, но до решения задачипока далековато. На текущий момент у меня закончились эффективные математические и алгоритмические идеи по ускорению вычислений. И пока что дальнейшее продвижение видится только в аккуратном переписывании программы на C/C++, распараллеливание на видеокарты и последующие нудные технические оптимизации. Так что если следующая часть и будет, то весьма не скоро. Ну или вдруг в голову нагрянет какая-нибудь свежая идея.

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

А ещё можно поугадывать формулу. Но тут всё совсем неопределённо.