В этой статье рассмотрим "Алгоритм X" Кнута и его применение для решения судоку. Прелесть алгоритма в том, что судоку при этом решается быстро без программирования каких-то продвинутых техник решения.


Началось всё, собственно, с задачки из Project Euler, где, чтобы получить ответ, нужно решить 50 судоку. И вроде ответ на неё получил, написав программку для решения довольно тупым перебором, но как-то осталась неудовлетворённость скоростью решения. Посмотрев, как решают судоку нормальные люди, я обнаружил, что сейчас для этого используется Алгоритм X, придуманный тем самым Дональдом Кнутом.


Алгоритм X


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


  • Есть множество S
  • Есть набор подмножеств Y этого множества
  • Задача состоит в том, чтобы выбрать из Y такие элементы Y*, что каждый элемент из S содержится только в одном из множеств, входящих в Y*
  • То есть объединение всех множеств в Y* и составляет (покрывает) множество S, и при этом в Y* нет попарно пересекающихся множеств

Например, рассмотрим множества
S = {1, 2, 3, 4, 5} и
Y = { A = {1, 2},
         B = {2, 3},
         C = {1, 5},
         D = {1, 4},
         E = {5} }
Множества B, D и E формируют точное покрытие множества S.


Для алгоритма X Кнута множество Y представляется в виде двоичной матрицы, где строки соответствуют элементам Y, и Ai,j = 1, если Sj находится в Yi. Т.е. для примера выше:


Yi \ Sj 1 2 3 4 5
A 1 1 0 0 0
B 0 1 1 0 0
C 1 0 0 0 1
D 1 0 0 1 0
E 0 0 0 0 1

Алгоритм поиска точного покрытия следующий:


  • Входные данные: множества S и Y; стэк Stack множеств, потенциально входящих в покрытие (может изначально быть пустой или уже иметь какие-то элементы)
    1. Если множество S пустое — на стэке лежит искомое покрытие.
    2. Если множество Y пустое — покрытие не найдено.
    3. Ищем в множестве S элемент s, входящий в минимальное число множеств из Y.
    4. Выбираем из Y все строчки X, содержащие s.
    5. Для каждого множества X повторяем 6-9.
    6. Добавляем X в стэк Stack как потенциальный элемент покрытия.
    7. Формируем множества S' и Y': S' — это S, из которого удалены все элементы, содержащиеся в X, Y' — множества из Y, не пересекающиеся с X.
    8. Вызываем алгоритм X для S', Y' и Stack.
    9. Если на шаге 7 получено, что покрытие невозможно — снимаем с вершины Stackа элемент и переходим к следующему X. Если решение найдено — возвращаем его.
    10. Если ни для какого X решения нет — для этих входных данных задача не решается.

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


Тонкость в том, что на практике этот алгоритм применяется для задач, где множества в Y — "маленькие", т.е. матрица весьма разреженная, из-за чего, например, поиск пересечений между столбцами при стандартном хранении в виде матрицы занимает непозволительно много времени.
Поэтому Кнут дополняет этот алгоритм механизмом "пляшущих ссылок". Матрица представляется в виде двумерного двусвязного списка: для каждой строки в списке хранятся только номера столбцов, где в этой строке содержатся единицы. Также в списке хранятся ссылки на следующий и предыдущий элемент в строке и столбце. Такая организация позволяет удалять из разреженной матрицы столбцы и строки за время O(1) по сравнению с O(m * n) при хранении в двумерном массиве.


Интересно, что Ali Assaf предлагает альтернативу механизму пляшущих ссылок с использованием ассоциативных списков, что позволяет на высокоуровневых языках реализовывать алгоритм X буквально в несколько десятков строчек.


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


Рассмотрим реализацию алгоритма на Julia.
Матрица из примера будет выглядеть теперь так:


Y = Dict(
    'A' => [1, 2],
    'B' => [2, 3],
    'C' => [1, 5],
    'D' => [1, 4],
    'E' => [5]
)

S = Dict(
    1 => Set(['A', 'C', 'D']),
    2 => Set(['A', 'B']),
    3 => Set(['B']),
    4 => Set(['D']),
    5 => Set(['C', 'E'])
)

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


function extract_intersects!(rows, columns, base_row)
    buf = []
    for elt in rows[base_row]
        # вынимаем текущий столбец из таблицы в буфер
        push!(buf, pop!(columns, elt))
        # удаляем все пересекающиеся строки из всех оставшихся столбцов
        for intersecting_row in buf[end]
            for other_elt in rows[intersecting_row]
                if other_elt != elt
                    delete!(columns[other_elt], intersecting_row)
                end
            end
        end
    end
    return buf
end

function restore_intersects!(rows, columns, base_row, buf)
    # удаляли столбцы от первого пересечения с base_row к последнему, восстанавливать надо в обратном порядке
    for elt in Iterators.reverse(rows[base_row])
        columns[elt] = pop!(buf)
        for added_row in columns[elt]
            for col in rows[added_row]
                push!(columns[col], added_row)
            end
        end
    end
end

Чтобы эти две функции работали как надо, как раз и требовалось упорядоченное хранение элементов в строках матрицы. В функции extract_intersects!() на каждой итерации внешнего цикла из матрицы убираются те строки, которые пересекаются с base_row, но не содержат элементов, просмотренных на предыдущих итерациях. Это гарантирует, что, когда мы в restore_intersects!() вставляем столбцы в обратном порядке, в самом внутреннем цикле на момент вызова push!(columns[col], added_row) столбец columns[col] в матрицу уже будет возвращён, и все удалённые в extract_intersects!() элементы из столбцов будут возвращены на прежнее место.


Теперь сам алгоритм X:


function algorithm_x(rows, columns, cover = [])
    if isempty(columns)
        return cover
    else
        # ищем столбец с минимальным числом элементов
        m, c = findmin(Dict(k => length(v) for (k, v) in columns))
        for subset in collect(columns[c])
            push!(cover, subset)
            # удаляем пересекающиеся подмножества и
            # содержащиеся в subset элементы
            buf_cols = extract_intersects!(rows, columns, subset)
            s = algorithm_x(rows, columns, cover)
            # если нашлось непустое решение - готово, выходим
            s == nothing || return s
            restore_intersects!(rows, columns, subset, buf_cols)
            pop!(cover)
        end
        # сюда дойдём либо если в columns[c] пусто,
        # либо когда рекурсивный поиск не нашёл решения
        return nothing
    end
end

Судоку


Алгоритм есть, дело за малым — представить судоку как задачу поиска точного покрытия.


Сформулируем требования, которым должно удовлетворять решённое судоку:


  1. В каждой клетке стоит цифра от 1 до 9 (или до n2, если решаются квадраты другого размера).
  2. В каждой строке каждое число встречается по разу.
  3. В каждом столбце каждое число встречается по разу.
  4. В каждом квадранте каждое число встречается по разу.

Каждое из этих требований должно выполняться ровно по 1 разу, т.е. они и формируют множество, которое надо покрыть. В нём ровно 4n2 элементов (столбцов в матрице).


Подмножества, которые рассматриваем, формируются подстановкой конкретного числа в конкретную клетку. Например, число 9 на пересечении 1 строки и 4 столбца "накрывает" подмножество "в клетке (1,4) есть число, в 1 строке есть число 9, в 4 столбце есть число 9, во 2 квадранте есть число 9" (подразумевая обычное судоку 9?9).


После этого алгоритм решения пишется тривиально.


# судоку задаётся матрицей 9?9, на месте неизвестных чисел нули
# идентификаторы строк - кортежи вида (row, col, num)
# идентификаторы столбцов:
# (0, row, col) - на пересечении row и col стоит число
# (1, row, num) - в строке row есть число num
# (2, col, num) - в столбце col есть число num
# (3, q, num) - в квадранте q есть число num
function xsudoku(puzzle::AbstractMatrix{Int})
    rows = Dict()
    cols = Dict()
    # заполняем строки
    for row in 1:9, col in 1:9, num in 1:9
        r = []
        quad = ((row-1)?3)*3 + (col-1)?3 + 1
        push!(r, (0, row, col), (1, row, num), (2, col, num), (3, quad, num))
        rows[(row, col, num)] = r
    end
    # заполняем столбцы
    for type in 0:3, n1 in 1:9, n2 in 1:9
        cols[(type, n1, n2)] = Set()
    end
    for (rk, rv) in rows
        for v in rv
            push!(cols[v], rk)
        end
    end

    # s - заготовка для ответа
    # для начала, туда надо внести те цифры, которые уже заполнены
    s = []
    for i in 1:9, j in 1:9
        if puzzle[i, j] > 0
            elt = (i, j, puzzle[i,j])
            push!(s, elt)
            # добавив клетку в решение, удаляем из матрицы все несовместимые элементы
            extract_intersects!(rows, cols, elt)
        end
    end

    # всё, что осталось - найти покрытие
    success = algorithm_x(rows, cols, s)
    success != nothing || return nothing
    # ответ выдадим в виде матрицы
    ret = similar(puzzle)
    for (i, j, num) in s
        ret[i,j] = num
    end
    return ret
end

Проверим на каком-нибудь примере:


 julia> @time xsudoku([0 0 0 0 0 0 4 0 0;
       3 0 6 0 0 0 0 0 0;
       0 0 0 1 9 6 0 3 0;
       0 7 0 0 0 0 0 1 0;
       8 0 0 2 5 0 0 9 0;
       0 4 0 0 0 0 8 0 0;
       0 6 0 4 0 9 0 0 8;
       0 0 5 0 0 0 0 2 0;
       0 0 0 5 0 0 0 0 7])
  0.006326 seconds (54.91 k allocations: 3.321 MiB)
9?9 Array{Int64,2}:
 1  5  7  8  3  2  4  6  9
 3  9  6  7  4  5  2  8  1
 2  8  4  1  9  6  7  3  5
 6  7  2  9  8  4  5  1  3
 8  3  1  2  5  7  6  9  4
 5  4  9  6  1  3  8  7  2
 7  6  3  4  2  9  1  5  8
 4  1  5  3  7  8  9  2  6
 9  2  8  5  6  1  3  4  7

Вроде работает, и скорость приемлемая.


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

Комментарии (14)


  1. Pochemuk
    04.08.2019 23:05

    Не на правах рекламы:

    Есть старая (но мне очень нравится) программа:

    www.angusj.com/sudoku

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

    К сожалению, справка (а в ней как раз и описаны пояснения к этим изящным приемам) будет работать только под Win9x/2K/XP. Для более старших надо скачать патч:

    support.microsoft.com/ru-ru/help/917607/feature-not-included-help-not-supported-error-opening-help-windows

    А для Win10 и патча такого нет…


    1. Pand5461 Автор
      05.08.2019 14:56

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

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


      Меня тут поразило, что настолько лобовой подход настолько эффективен при использовании структуры данных, позволяющей быстрый откат из рекурсии.


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


      1. Pochemuk
        05.08.2019 18:13

        … найти клетку с минимумом возможных вариантов

        Кстати, а почему именно с минимумом?

        Как изменится эффективность алгоритма, если брать случайную клетку или даже с максимумом вариантов?


        1. Pand5461 Автор
          05.08.2019 18:50

          Ну это просто жадная стратегия поиска.
          Логика, как мне кажется, тут следующая. Поиск ограничения, которое можно удовлетворить минимальным числом способов, всё равно нужен. Потому что если этот минимум ноль — то ветка тупиковая, и надо возвращаться. Если он равен единице — то элемент, удовлетворяющий ограничению, точно входит в покрытие. Ну а дальше, раз минимум всё равно нашли, давайте его и будем использовать.
          В терминах О-большого асимптотика худшего случая не меняется, т.к. задача NP-полная. Для лучшего случая просто получается, что поиск одиночек не нужно программировать как какую-то особую стратегию, выбор клетки с минимальным числом вариантов именно их в первую очередь и будет заполнять.


  1. domix32
    05.08.2019 12:30

    Интересно было бы посмотреть на рост времени в зависимости от сложности.


    1. Pand5461 Автор
      05.08.2019 14:02

      От сложности в каком смысле? Уровень сложности классического судоку или при переходе к полям 16?16, 25?25 и т.д.?


      В первом смысле: из той самой задачки с Project Euler самое простое судоку решается за 3.8 мс, самое сложное — за 4.6. Пример из комментария про "самое сложное судоку" в статье mrHobbY решается за 69 мс. Если в функциях нормально прописать аннотации типов, то цифры быстрее раза в 3-4, но соотношение между самым простым и самым сложным вариантом на порядок так и остаётся.


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


      1. Pochemuk
        05.08.2019 19:12

        Эффективная очередь это как? Отсортированный по количеству возможных покрытий массив клеток? Но тогда придется его перестраивать при каждой попытке. Перестроение тоже не O(n) потребует, а поболее. Несмотря на то, что он уже частично отсортирован.

        И таких перестроений будет превеликое множество. А еще возвраты надо делать после откатов из тупиков.

        Так что, слухи об эффективности O(N*log(N)) сильно преувеличены.


        1. Pand5461 Автор
          05.08.2019 19:47

          Если говорить про лучший случай, когда всё решается последовательным заполнением одиночек, то тупиков не будет. Это я специально оговорил.
          Эффективная очередь — это когда гарантирован только минимум в голове, а дальше может быть что угодно. Тогда смену приоритета одного элемента можно сделать за O(logN). Каждая строка в матрице имеет не больше 4 единиц, а столбец — не больше N1/2. Т.е. дописали циферку в решение — удаляем O(N1/2) строк из матрицы и тратим ещё O(N1/2logN) на обновление приоритетов в очереди.
          Ну да — O(N3/2logN) (где N — число клеток, а цифр — N1/2). Но всё-таки лучше, чем квадрат. Хотя на классическом поле 9?9 все эти очереди с перестроениями могут тупо против линейного поиска ещё не сыграть.


          1. Pochemuk
            05.08.2019 20:37

            Не совсем понял по поводу вычеркивания…

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

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


            1. Pand5461 Автор
              05.08.2019 21:36

              Ну да, примерно так.


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


              Только на самом деле там кроме ситуации "на клетку Г, В осталось мало кандидатов" нужно просмотреть "на горизонтали Г для цифры А осталось мало возможных мест" и аналоги для вертикали и квадрата. Вычёркивается, соответственно, не только цифра из кандидатов в другие клетки, но и вертикаль из доступных для других цифр на этой горизонтали и т.д.


      1. domix32
        06.08.2019 09:39

        Более сложные судоку имеют меньшее количество стартовых цифр. Собственно все. Генерализация на разные основания систем счисления уже не настолько интересно.


  1. Kyushu
    05.08.2019 16:49

    Как соотносится алгоритм с существованием и единственностью решений?


    1. Pand5461 Автор
      05.08.2019 18:54

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


  1. spice_harj
    07.08.2019 15:55

    Просто удивительно! Иногда, кажется, что Хабр резонирует… Статья появилась, как раз в то время, когда я занялся решением задачки «exact cover problem» и был в потенциальном поиске алгоритмов (ещё не знал к какому классу принадлежит задача).
    С псевдо-кодом было бы великолепно! (не все мы Julia)