В этой статье рассмотрим "Алгоритм 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, и A_i,j = 1, если S_j находится в Y_i. Т.е. для примера выше:

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

• Входные данные: множества 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 (или до n², если решаются квадраты другого размера).
2. В каждой строке каждое число встречается по разу.
3. В каждом столбце каждое число встречается по разу.
4. В каждом квадранте каждое число встречается по разу.

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

Подмножества, которые рассматриваем, формируются подстановкой конкретного числа в конкретную клетку. Например, число 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

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

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