Классификация комбинаторных объектов на примере латинских квадратов

Хочу поделиться некоторым опытом по написанию программ для перечисления комбинаторных объектов из заданного класса (например, латинских квадратов, систем Штейнера, кодов исправляющих ошибки, ортогональных массивов). Обычно нужно перечислить все объекты с заданными параметрами, например, таблицы заданного размера, заполненные числами согласно некоторому правилу. Под словом «все» можно подразумевать как «все различные», так и «принципиально различные» в смысле, специфическом для конкретной задачи, например, таблицы могут считаться принципиально одинаковыми (эквивалентными), если одна получается из другой перестановкой строк. Число объектов, полученных в результате подобного перечисления, может варьироваться от 0 (это также имеет смысл, т.к. в результате мы устанавливаем, что объекты с данными параметрами не существуют) до десятков миллиардов. Обычно вычисления подобного рода являются частью теоретического исследования, но аналогичные подходы могут использоваться и в практических задачах.

Язык программирования в рассматриваемых примерах – Sage, по сути это расширение Python, текст программ в примерах из данной заметки не отличается от пайтоновского за исключением использования некоторых библиотечных функций. Примеры могут быть запущены на онлайн-калькуляторе https://sagecell.sagemath.org/.

Описание задачи и подхода

Начнем с описания частной задачи, на примере которой мы и будем рассматривать общий алгоритм перечисления. Мы попробуем посчитать все латинские квадраты порядка 7. Латинский квадрат порядка n – это таблица размера n-на-n, заполненная n символами так, что в каждой строке и в каждом столбце все символы встречаются по одному разу. Мы будем использовать символы от 0 до n−1, хотя Эйлер пользовался для этой цели буквами латинского алфавита (отсюда и название). Пример латинского квадрата порядка 5:

$\begin{array}{|c|c|c|c|c|}\hline0&1&2&3&4\\\hline1&0&3&4&2\\\hline2&4&1&0&3\\\hline3&2&4&1&0\\\hline4&3&0&2&1\\\hline\end{array}$

Количество всех латинских квадратов порядка 7 – четырнадцатизначное число, и хотя в принципе возможно такое количество сгенерировать и даже где-то сохранить, мы несколько упростим задачу. Мы будем перечислять квадраты «с точностью до эквивалентности», то есть, если несколько различных квадратов получаются друг из друга некоторыми естественными преобразованиями, мы объединим их в одну группу, класс эквивалентности, и из каждого класса эквивалентности будем сохранять по одному представителю. В качестве таких естественных преобразований эквивалентности мы возьмем всевозможные комбинации перестановок строк квадрата, перестановок столбцов, и перестановок символов алфавита, который используется для заполнения ячеек. Очевидно, что применив любую из этих операций к латинскому квадрату, мы получим опять латинский квадрат. Операции рассматриваемого типа называются изотопиями, а латинские квадраты, получаемые друг из друга изотопиями, называются изотопными; соответствующие классы эквивалентности будем называть изотопическими классами. В принципе, можно расширить рассматриваемые преобразования эквивалентности, добавив например возможность отразить квадрат по диагонали (смена ролей строк и столбцов, аналогично можно менять роли столбцов, строк и символов), тогда число классов эквивалентности будет еще меньше, но мы начнем с изотопий.

Для примера, легко проверить, что все латинские квадраты порядка 3 изотопны. Латинские квадраты порядка 4 делятся уже на два изотопических класса, с представителями

$\begin{array}{|c|c|c|c|}\hline0&1&2&3\\\hline1&0&3&2\\\hline2&3&0&1\\\hline3&2&1&0\\\hline\end{array}$ и $\begin{array}{|c|c|c|c|}\hline0&1&2&3\\\hline1&2&3&0\\\hline2&3&0&1\\\hline3&0&1&2\\\hline\end{array}$ .

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

Общая схема перечисления состоит из следующих этапов:

I. Порождение объектов.
II. Проверка на эквивалентность (в нашем случае – изотопность) и сохранение ровно одного представителя из каждого класса эквивалентности.
III. Численная проверка корректности вычислений.

Можно написать программу, которая сразу породит все латинские квадраты маленького порядка, возможно даже 6, их число «всего» 812851200. Но 14-значное число всех латинских квадратов порядка 7 …, даже просто пересчитать такое число объектов в разумное время проблематично, не говоря уже о том, чтобы сгенерировать и сравнить на изотопность. Поэтому мы разобьем задачу на несколько шагов, на каждом шаге перечисляя некоторые промежуточные объекты, как бы «недостроенные» латинские квадраты. В кажестве таких промежуточных объектов идеально подходят латинские прямоугольники. Латинский прямоугольник – это недостроенный латинский квадрат, у которого не заполнено несколько нижних строк. Формально, таблица размера m на n, заполненная числами от 0 до n−1, называется латинским прямоугольником, если ни одна строка и ни один столбец не содержит повторяющихся чисел. Теорема Холла гласит, что любой латинский прямоугольник можно достроить до латинского квадрата, но нам этот замечательный факт не понадобится.

Итак, мы по порядку будем перечислять все неизотопные латинские прямоугольники, сначала 1 на n, потом 2 на n, 3 на n, и т.д., до прямоугольников n на n, которые являются латинскими квадратами. На каждом шаге мы будем повторять этапы I, II и III, но на этапе I шага m будем строить не все латинские прямоугольники, а те, которые являются продолжениями попарно неизотопных представителей прямоугольников m−1 на n, полученных на предыдущем шаге вычислений (по сути, нам нужно к каждому из имеющихся прямоугольников m−1 на n приписать одну строку всеми возможными способами, не противоречащими определению).

Например, предположим для n=4 и m=2 мы нашли два неизотопных латинских прямоугольника

$\begin{array}{|c|c|c|c|}\hline0&1&2&3\\\hline1&0&3&2\\\hline\end{array}$ и $\begin{array}{|c|c|c|c|}\hline0&1&2&3\\\hline1&2&3&0\\\hline\end{array}$ .

На следующем шаге m=3, добавляя на этапе I третью строку всеми возможными способами к каждому из них, мы получим

$\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&0&3&2\\ \hline2&3&0&1\\ \hline\end{array}$ , $\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&0&3&2\\ \hline2&3&1&0\\ \hline\end{array}$ , $\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&0&3&2\\ \hline3&2&0&1\\ \hline\end{array}$ , $\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&0&3&2\\ \hline3&2&1&0\\ \hline\end{array}$

и $\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&2&3&0\\ \hline2&3&0&1\\ \hline\end{array}$ , $\begin{array}{|c|c|c|c|} \hline0&1&2&3\\ \hline1&2&3&0\\ \hline3&0&1&2\\ \hline\end{array}$ .

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

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

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

Пример: рассмотрим два графа на одном множестве вершин {0,1,2,3}. У первого графа множество ребер {{0,1},{0,2},{1,2},{1,3}}, у второго – {{0,2},{1,2},{1,3},{2,3}}; отображение 0→1, 1→2, 2→3, 3→0 является изоморфизмом первого графа на второй, а отображение 0→2, 1→1, 2→0, 3→3 является автоморфизмом первого графа (но не второго).

Шаг I: Порождение латинских прямоугольников m на n, продолжающих данный латинский прямоугольник m−1 на n

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

Начнем со сведения к задаче о поиске клик в графе. Рассмотрим для примера латинский прямоугольник 2 на 5

$\begin{array}{|c|c|c|c|c|}\hline 0&1&2&3&4\\\hline 1&0&3&4&2\\\hline \phantom{1}&&&&\\\hline \end{array}$ .

Мы хотим всеми возможными способами продолжить его до латинского прямоугольника 3 на 5. Рассмотрим столбец с номером 0. В нем уже стоят числа 0 и 1, поэтому в третьей строке может стоять только 2, 3, или 4. Мы запишем эти возможности как (0,2), (0,3) и (0,4), где 0 на первом месте – это номер столбца, а на втором месте стоит символ, который можно вписать в этот столбец. Для столбца номер 1 имеем (1,2), (1,3), (1,4), далее (2,0), (2,1), (2,4), потом (3,0), (3,1), (3,2) и наконец (4,0), (4,1), (4,3). Итого имеем 15 пар (0,2), (0,3), (0,4), (1,2), (1,3), (1,4), (2,0), (2,1), (2,4), (3,0), (3,1), (3,2), (4,0), (4,1), (4,3). Чтобы дописать к данному латинскому прямоугольнику еще одну строку, нужно выбрать из этих 15 пар пять пар, у которых различны все первые компоненты (это означает, что мы впишем по одному новому символу в 5 различных столбцов) и все вторые компоненты (это означает, что в новой строке все значения различны). Для этого на этих 15 парах мы построим «граф совместимости», в котором две пары будут образовывать ребро если у них различны как первые компоненты, так и вторые. Клики размера 5 в этом графе как раз и будут соответствовать возможностям заполнить третью строку. Напишем теперь генератор, который будет делать то же самое в общем случае.

def addRow(LR,algorithm=0):     """Генератор всех латинских прямоугольников m-на-n,        продолжающих данный латинский прямоугольник (m-1)-на-n"""     n = len(LR[0]) # ширина прямоугольника     N = set(range(n))     # для каждого столбца, записываем не использованные символы:     available = ( N-set(col) for col in zip(*LR) )    # список доступных пар (столбец,символ):     candidates = [ (i,j) for i,avail in enumerate(available) \                          for j in avail ]     if algorithm==0: # при помощи поиска клик          # построение графа совместимости         compatible_pairs = [ [(i,j),(ii,jj)] \                              for i,j in candidates \                              for ii,jj in candidates \                              if i<ii and j!=jj ]         G = Graph(compatible_pairs) # Sage граф         if G.clique_number()==n: # если есть клики размера n             for C in G.cliques_maximum():                 # каждую такую клику преобразуем в строку                 yield LR + [[j for i,j in sorted(C)]]

Этот код использует встроенный в Sage тип «граф» и два его метода clique_number (возвращает максимальный размер клики в графе) и cliques_maximum (возвращает все клики максимального размера).

Теперь второй способ. Аналогично предыдущему, для рассматриваемого примера строим доступные 15 пар (столбец,символ). Заметим, чтокаждая пара «использует» один столбец и один символ. Нам нужно вы-брать 5 пар таких, чтобы каждый столбец и каждый символ был использован ровно один раз. Это классическая задача о точном покрытии (exact cover). Есть множество (в нашем примере – множество мощности 10, пять столбцов и пять символов) и набор его подмножеств (в нашем случае, подмножеств мощности 2, состоящих из одного столбца иодного символа); нужно выбрать поднабор покрывающий данное множество однократно, то есть чтобы каждый элемент данного множествавстречался ровно в одном подмножестве из поднабора. Такая задача решается алгоритмом Dancing Links Дональда Кнута, реализованным для разных платформ. Кроме того, она является частным случаем задачи целочисленного линейного программирования, ILP, для которого есть множество решателей (однако, многие находят только одно решение, а нам нужны все). Добавим альтернативный алгоритм для генератора addRow, используя решатель задачи о точном покрытии AllExactCovers, встроенный в Sage.

    if algorithm!=0: # при помощи точного покрытия         # строим exact-cover матрицу M         M = []         for i,j in candidates:             M.append([0]*2*n)             # первые n элементов нумеруют столбцы, последние n - символы             M[-1][i] = M[-1][n+j] = 1         for cover in AllExactCovers(matrix(M)):             yield LR + [[list(line).index(1,n)-n \                          for line in reversed(sorted(cover))]]

Шаг II: проверка на изотопность

Как было уже упомянуто выше, проверка на эквивалентность комбинаторных объектов часто сводится к проверке на изоморфность соответствующих графов (в общем случае как определить «соответствующие графы» для данного класса объектов может оказаться нетривиальной задачей, иногда решений несколько, и они сильно отличаются по скорости сравнения графов на изоморфность). Для латинских прямоугольников m на n мы будем строить граф следующим образом. Вершин mn+n+n+m, в качестве вершин берем mn ячеек латинского прямоугольника, m строк, n столбцов, и n символов (каким-нибудь образом закодированные, например числами от 0 до mn+n+n+m−1); каждую вершину-ячейку соединяем одним ребром с вершиной-строкой, в которой эта ячейка находится, одним ребром с вершиной-столбцом, и одним ребром с вершиной-символом, который записан в данной ячейке.

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

Доказывать теорему не будем, факт достаточно простой. Она сводит изотопность латинских прямоугольников к изоморфности соответствующих графов, но с учетом разбиения вершин по ролям (изоморфизм между графами должен сопоставлять вершины-ячейки вершинам-ячейкам, вершины-строки вершинам-строкам и т.д.). К счастью, программы по распознаванию изоморфности графов умеют работать с разбиениями. На практике проверка изоморфности графов обычно происходит следующим образом: для каждого графа строится изоморфный ему канонический граф. Для изоморфных графов канонические графы одинаковы, поэтому их можно сравнивать уже просто на равенство. Выгода от этого подхода следующая. Предположим, что у нас есть 449 графов и мы хотим проверить, что все они неизоморфны. Если мы просто будем проверять каждую пару на изоморфность, нам понадобится 449 · 448/2 ∼ 100500 таких операций, каждая из которых довольно трудоемка. Вместо этого мы можем всего 449 раз построить канонический граф для каждого из данных графов, и провести 100500 уже очень легких сравнений на равенство. Итак, запишем подпрограмму nonIsotopic, которая из данного набора латинских прямоугольников выбирает неизотопные представители и возвращает список таких представителей. В виде отдельной функции rectangleToGraph мы выделим построение по прямоугольнику графа, при помощи встроенной в Sage функции Graph.

def rectangleToGraph(L,m,n):     """     Возвращает характеристический граф      латинского m-на-n прямоугольника L     и разбиение множества вершин по ролям     """     V = [ (i,j) for i in range(m) for j in range(n) ] # вершины-ячейки     Vr = [ ("r",i) for i in range(m) ] # вершины-строки     Vc = [ ("c",j) for j in range(n) ] # вершины-столбцы     Vs = [ ("s",k) for k in range(n) ] # вершины-символы     E  = [ [(i,j),Vr[i]] for i,j in V ] # ребра ячейка-строка     E += [ [(i,j),Vc[j]] for i,j in V ] # ребра ячейка-столбец     E += [ [(i,j),Vs[L[i][j]]] for i,j in V ] # ...ячейка-символ     P = [V,Vr,Vc,Vs] # разбиение множества вершин     # E += [ [U[i],U[j]] for U in [Vc,Vs] \     #          for j in range(n) for i in range(j) ] # понадобится позже     # P = [V,Vr,Vc+Vs] # понадобится позже     return Graph(E), P # graph and vertex partition

def nonIsotopic(LRs):     """     На входе итератор, который выдает латинские прямоугольники LR.     На выходе -- неизотопные представители      """     GG = {} # словарь порядковых номеров найденных канонических графов     Repres = [] # неизотопные представители латинских прямоугольников     for L in LRs:         m, n = len(L), len(L[0])         G, P = rectangleToGraph(L,m,n) # граф и разбиение         Gc = G.canonical_label(partition=P).copy(immutable=True) # канонический         if Gc not in GG:             GG.update({Gc:len(Repr)})             Repres += [L]     return Repres

Шаг III: проверка корректности вычислений

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

При описанном перечислительном подходе контроль правильности рассчетов можно провести, двумя способами вычислив мощность каждого найденного класса эквивалентности K (в нашем случае, изотопического класса) латинских прямоугольников m на n. С одной стороны, мы знаем сколько раз латинские прямоугольники из этого класса нам встретились при попытке продолжить каждый латинский прямоугольник m−1 на n. Это нам дает формулу $M(K) = \sum_{L} M_L \cdot S_L(K)$ , где суммирование происходит по представителям L изотопических классов латинских прямоугольников m−1 на n, найденных на предыдущем этапе, M_L – мощность соответствующих классов, и – число, сколько раз прямоугольник из класса K встретился среди продолжений прямоугольника L.

С другой стороны, мощность M(K) изотопического класса латинских прямоугольников m на n можно вычислить через число автотопий Aut(T) любого из его представителей T, $T\in K$ . А именно, $M(K) = m! \cdot n! \cdot n! / Aut(T)$ . Действительно, число всевозможных изотопий $m! \cdot n! \cdot n!$ (число перестановок строк , число перестановок столбцов , число перестановок символов ), из них Aut(T) изотопий дают один и тот же прямоугольник.

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

Добавим нужную проверку прямо в конец функции nonIsotopic, вычислив предварительно необходимые величины (новые строки – с переменными Auts, Count, parentCount).

def nonIsotopic(LRs):     """     На входе итератор, который выдает латинские прямоугольники LR     вместе с величиной parentCount -- мощностью изотопического     класса прямоугольника меньшей высоты, продолжением которого     получен LR.     На выходе -- неизотопные представители      и мощности соответствующих изотопических классов     """     GG = {} # словарь порядковых номеров найденных канонических графов     Repr = [] # неизотопные представители латинских прямоугольников     Auts = [] # число автотопий для каждого представителя     Count = [] # мощность каждого изотопического класса     for LR,parentCount in LRs:         m, n = len(LR), len(LR[0])         G, P = rectangleToGraph(LR,m,n)         Gc = G.canonical_label(partition=P).copy(immutable=True) # канонический         if Gc not in GG:             GG.update({Gc:len(Repr)})             Repr += [LR]             Auts += [G.automorphism_group(partition=P).order()]             Count += [0]         Count[GG[Gc]]+=parentCount # еще раз встретился     for cnt,aut in zip(Count,Auts):         cnt0 = factorial(n)**2 * factorial(m) / aut         # cnt0 = factorial(n)**2 * factorial(m) * 2 / aut # понадобится позже         assert cnt==cnt0, "error m="+str(m)+": "+str(cnt)+"!="+str(cnt0)     return Repr, Count

Заключительные шаги

Осталось все собрать. Представителя латинских прямоугольников 1 на n подставляем вручную, далее поочереди добавляем строку при помощи addRow() и находим неизотопных представителей при помощи nonIsotopic(), куда также передаем мощности родительских изотопических классов, для проверки корректности вычислений.

def findLRs(n,alg=0):     # представители латинских m*n прямоугольников, m=0,1,2,...:     LRs = [ [ [] ], [ [list(range(n))] ]  ]     # мощность соответствующих изотопических классов:     CNTs = [ [1], [factorial(n)] ]     for m in range(1,n):         LRs_nxt,CNTs_nxt = nonIsotopic( (LRnxt,cnt) \                                         for LR,cnt in zip(LRs[-1],CNTs[-1]) \                                         for LRnxt in addRow(LR,alg) )         LRs += [LRs_nxt]         CNTs += [CNTs_nxt]         print("n:",n,"m:",m+1," classes:",len(LRs_nxt)," total:",sum(CNTs_nxt))     return LRs  LRs = findLRs(7)

Копируем все (определения функций addRow, rectangleToGraph, nonIsotopic, findLRs, и строку LRs=findLRs(7)) в форму онлайн-калькулятора https://sagecell.sagemath.org/, запускаем, нажав кнопочку «Evaluate», и …

m: 2 classes: 4 total: 9344160
m: 3 classes: 56 total: 5411750400
m: 4 classes: 1398 total: 782137036800
m: 5 classes: 6941 total: 20449013760000

… о, нет! За выделенное на ресурсе время программа не успевает завершить выполнение, посчитав только до прямоугольников высоты 5, иногда, даже 6. Не хватает буквально нескольких секунд. Но не будем расстраиваться. Во-первых, подставив найденные значения в онлайн энциклопецию целочисленных последовательностей
https://oeis.org/search?q=56%2C1398%2C6941&go=Search
https://oeis.org/search?q=782137036800+20449013760000&go=Search
мы убеждаемся что наша программа выдает верные величины как для числа изотопических классов, так и для числа всех латинских прямоугольников. Во-вторых, программа успешно завершается для меньших порядков, n=6, n=5, … ( for n in (3,4,5,6): findLRs(n) ). А для n=7 мы попытаемся несколько ускорить программу, расширив эквивалентность. Кроме изотопий, мы разрешим менять ролями столбцы и символы. А именно, для латинского прямоугольника A мы всегда можем построить латинский прямоугольник B такой, что A(i,j)=k (эта запись означает, что на пересечении i-й строки и j-го столбца прямоугольника A стоит символ k) тогда и только тогда, когда B(i,k)=j. Такие прямоугольники не обязательно изотопны, и если раньше мы их не объединяли в один класс эквивалентности, то теперь будем. Это уменьшит число классов эквивалентности и, хочется надеяться, общее время перебора. Нам нужно внести следующие два изменения (можно просто раскомментировать в двух местах строки с пометкой «понадобится позже»). Во-первых, предпоследнюю строку определения функции rectangleToGraph заменим на

    E += [ [U[i],U[j]] for U in [Vc,Vs] \                       for j in range(len(U)) for i in range(j) ]     P = [V,Vr,Vc+Vs] # разбиение множества вершин

Наше новое преобразование эквивалентности меняет местами множества вершин-столбцов и вершин-символов, поэтому мы их объединяем в разбиении P. С другой стороны, нам нужно избежать преобразований, которые перемешивают эти два множества. Именно для этого мы добавляем дополнительные ребра в список ребер E, делая каждое их этих множеств кликой. Второе изменение – при подсчете мощности класса через число автотопий в конце определения функции nonIsotopic. Теперь у нас преобразований эквивалентности вдвое больше, чем автотопий, поэтому мы вставим эту двойку в формулу:

        cnt0 = factorial(n)**2 * factorial(m) * 2 / aut

Запускаем.

n: 7 m: 2 classes: 4 total: 9344160
n: 7 m: 3 classes: 45 total: 5411750400
n: 7 m: 4 classes: 808 total: 782137036800
n: 7 m: 5 classes: 3712 total: 20449013760000
n: 7 m: 6 classes: 1895 total: 61479419904000
n: 7 m: 7 classes: 324 total: 61479419904000

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

Выводы и заключительные замечания

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

Когда дело реально упирается в вычислительные ресурсы, приходится переходить на низкоуровневые языки программирования, а потом при необходимости пытаться подобрать ключик к конкретной задаче. Пара универсальных рецептов: 1. Чтобы сэкономить место на хранении канонических графов, можно их обратно преобразовать в латинские прямоугольники (ну или что мы там классифицируем) и хранить канонические латинские прямоугольники. 2. Часто бывает полезно классифицировать отдельно объекты, содержащие некоторую подконфигурацию, и объекты избегающие ее. Например, для латинских прямоугольников такой подконфигурацией мог бы быть подквадрат 2 на 2, содержащий всего два символа. 3. Иногда большое количество промежуточных решений не имеет продолжения. Бывает полезно попытаться отсечь их все или большую часть каким-то более быстрым способом, чем решателем, которым мы находим все продолжения в общем случае (например, можно поэкспериментировать с решателями задач линейного программирования, MILP).

При написании программ на C/C++ я использую следующие библиотеки:

проверка графов на изоморфность – Nauty&Traces https://pallini.di.uniroma1.it/
поиск клик в графе – Cliquer https://users.aalto.fi/~pat/cliquer.html
нахождение точного покрытия (в том числе о кратного) – Libexact https://pottonen.kapsi.fi/libexact.html

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

P. Kaski and P. R. J. Östergård. Classification Algorithms for Codes and Designs, vol. 15 of Algorithms Comput. Math. Springer, Berlin, 2006. https://doi.org/10.1007/3-540-28991-7

В заключение, пара замечаний про Sage для «начинающих». Обычный код Python может не работать в Sage из-за некоторых различий (однако, python код можно корректно импортировать в Sage). Самые заметные различия: Во-первых, символ ^ обозначает возведение в степень, а не исключающее или; исключающее или – ^^. Во-вторых, целочисленные константы трактуются как объекты класса sage.rings.integer.Integer, для которого операции определены по-своему и вымолняются значительно дольше, чем для int. Чтобы обозначить константу типа int, можно использовать суффикс r, например 100500r.

Если работать на своем компьютере, то к Sage лучше установить пакет bliss – с ним будет гораздо быстрее проходить построение канонического графа (без этого пакета пример с латинскими квадратами может работать часы, а не минуты). Если Sage установлен в отдельную папку (рекомендуется), а не из репозитория, то в этой папке набрать ./sage -i bliss. Дистрибутивы и инструкции можно найти на https://www.sagemath.org/

ссылка на оригинал статьи https://habr.com/ru/articles/824886/