Всем доброго времени суток.

В этой статье я опишу метод получения нормальных магических квадратов порядка n^m, где n и m - положительные натуральные числа, при условии, что нам известен нормальный магический квадрат порядка n

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

Как это начиналось

Я никогда не изучал магические квадраты, составленные ранее, изображения которых были в открытом доступе, и не следовал известным "рецептам".

Алгоритм, которого я придерживался был прост, хотя и весьма сложен в исполнении, учитывая, что все это происходило на обычной тетрадной бумаге в клетку:

1. Рисуем квадрат и вписываем в него числа от 1 до n², где n - сторона квадрата, по порядку, построчно, слева направо, сверху вниз.

2. Подписываем исходные суммы в строках, столбцах и диагоналях.

3. Выписываем два столбца, расположенных симметрично относительно центра квадрата.

4. Меняем местами числа в столбцах до тех пор, пока суммы чисел не окажутся равны (следуя интуиции и расчетам). Здесь стоит отметить, что этот алгоритм подходит только для ассоциативных магических квадратов, но об этом я тогда не знал, из-за чего серьезно застрял на квадрате 6х6.

5. Рядом с исходным квадратом рисуем пустой квадрат, в котором стрелками соединяем клетки, числа в которых "переехали".

6. Применяем схему перемещения чисел к строкам, находящимся на таком же удалении от центра и проверяем суммы в строках.

7. Хотел бы я сформулировать и этот пункт, но, боюсь, здесь начинается подгонка: скрупулезно подсчитывая суммы, менять местами числа, находящиеся на одинаковом удалении от центра, считая по количеству строк/столбцов.

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

С маленькими квадратами все было донельзя просто (n=3 и n=4), и меня заинтриговали схемы, получающиеся в итоге. Это были симметричные рисунки, в которых мне виделось нечто загадочное. Загадкой для меня было скорее не то, отчего они так выглядят, а то, какие рисунки получатся от квадратов с большими сторонами, отчего я сразу перешел к квадратам 5х5, 6х6, 7х7 и 9х9.

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

Закончив с квадратом 9х9, я понял, что бумаги уходит реально много: на один тетрадный лист помещалось только три-четыре таких квадрата, а учитывая мою любовь к зарисовыванию стрелочками "путей чисел", места оставалось еще меньше.

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

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

Наши дни

Забросив рисование квадратов на первом курсе, отучившись, отслужив в армии и устроившись на работу, в один из долгих зимних вечеров, от нечего делать, я решил систематизировать свои знания о магических квадратах. К моему большому огорчению, все мои записи о магических квадратах пропали, и я едва не опустил руки. Нарисовав самый простой квадрат 3х3, я нарисовал рядом его схему перемещения, или "шаблон". Как истинно ленивый программист, я подумал: что, если применить этот шаблон к базовому квадрату 9х9, разделив его на 9 секторов, по три строки и столбца, а затем применить этот же шаблон к каждому сектору. Что, если получится магический квадрат?

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

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

NumPy + PyGame = Profit!

Итак, собственно, код. Для написания программы использовался Python версии 3.9 и opensource-библиотеки Numpy и PyGame.

import numpy as np
import pygame as pg
# SortedList uses binary search instead full-list comparison cycle.
import sortedcontainers as scs

Базовый квадрат для магического - это элементарная матрица, где m_x,y == x + y * n + 1, где n - длина стороны квадрата. Напишем функцию-генератор:

def create_matrix(n: int):
    return np.arange(1, n ** 2 + 1).reshape(n, n)

Далее я записал шаблоны магических квадратов 3x3 и 4x4, в виде списка последовательностей перестановок:

m_3 = [[[1, 0], [0, 0], [1, 2], [2, 2]], [[2, 0], [2, 1], [0, 2], [0, 1]]]
m_4 = [[[0, 0], [2, 2]], [[1, 0], [0, 1]],
       [[2, 0], [3, 1]], [[3, 0], [1, 2]],
       [[0, 2], [1, 3]], [[1, 1], [3, 3]],
       [[2, 1], [0, 3]], [[3, 2], [2, 3]]]

Затем, после вдумчивого разбора предлагаемого NumPy функционала, были созданы функции для разбивания матрицы на сектора, и обратной сборки:

# _m - input matrix, _x - base row length
def m_split(_m, _x):
    return [np.vsplit(x, _x) for x in np.hsplit(_m, _x)]


# _m - input "previous splitted" matrix, _x - base row length
def m_join(_m, _x):
    return np.vstack([np.hstack(np.vstack(_m)[x::_x]) for x in range(_x)])

После чего реализованы функции, выполняющие последовательную перестановку и цепь перестановок, в соответствии с шаблоном:

# _m - input "previous splitted" matrix, idxs - indices of matrix cells.
def m_rot(_m, idxs):
  	#      [-----X----][-----Y----]
    a0 = _m[idxs[0][0]][idxs[0][1]]
    for i in range(len(idxs) - 1):
        _m[idxs[i][0]][idxs[i][1]] = _m[idxs[i + 1][0]][idxs[i + 1][1]]
    _m[idxs[-1][0]][idxs[-1][1]] = a0
    return _m


#_m - input "previous splitted" matrix, _steps - pattern of digits swaps (m_3, m_r, ...)...
def magic_steps(_m, _steps):
    m0 = _m
    for step in _steps:
        m0 = m_rot(m0, step)
    return m0

И, наконец, функция генерации магического квадрата со стороной n^pow:

# _mss - magic steps (m_3, m_4, or ...)
def pow_square(_root, _pow, _mss):
    m_dmag = lambda x, d: m_join(magic_steps(m_split(x, d), _mss), d)
    _p = _root ** _pow
    m0 = m_dmag(create_matrix(_p), _root)
    for i in range(1, _pow):
        m0 = m_join([[m_dmag(x, _root) for x in y] for y in m_split(m0, _root ** i)], _root ** i)
    return m0

Чтобы потешить ностальгию, и нарисовать глобальный шаблон для полученного на выходе квадрата, я использовал библиотеку PyGame, как в разы более шуструю альтернативу модулю turtle, который я по недоразумению попытался использовать вначале:

# New to_star function computes and draw all lines from matrix on fly.
def to_star(_m):
    dim = len(_m)
    used = scs.SortedList()
    max_sz = 980
    tp = 10	# brush transparency.
    sz = float(max_sz / dim)
    j_to_xy = lambda _j, _d: [(_j - 1) % _d, int((_j - 1) / _d)]
    # [matrix position] -> (drawing field point)
    c_to_pg = lambda _xy: (float(_xy[1] * sz + 10), float(_xy[0] * sz + 10))
    # pygame prepare:
    pg.init()
    pg.fastevent.init()
    d_sz = (max_sz + 20, max_sz + 20)
    _sc = pg.display.set_mode(d_sz)
    _sc.fill((255, 255, 255))
    pg.display.flip()
    c = pg.time.Clock()
    # start compute & draw
    for _y in range(dim):
        for events in pg.fastevent.get():
            if events.type == pg.QUIT:
                pg.quit()
                return
        sc = _sc.convert_alpha()
        sc.fill([0, 0, 0, 0])
        for _x in range(1, dim + 1):
            j = _x + _y * dim
            x, y = j_to_xy(j, dim)
            if _m[y][x] == j or _m[y][x] in used:
                continue
            a = _m[y][x]
            x, y = j_to_xy(a, dim)
            b = _m[y][x]
            used += [a]
            pg.draw.line(sc, (0, 0, 0, tp), c_to_pg([x, y]), c_to_pg(j_to_xy(b, dim)))
        _sc.blit(sc, (0, 0))
        pg.display.update()
    sc = _sc.convert_alpha()
    sc.fill([0, 0, 0, 0])
    for j in range(1, dim * dim + 1):
        x, y = j_to_xy(j, dim)
        if _m[y][x] == j:
          	# static point transparency is 50% higher than lines transparency.
            pg.draw.circle(sc, (255, 255, 255, (tp + 50) % 100), c_to_pg([x, y]), 1.5)
    _sc.blit(sc, (0, 0))
    pg.display.update()
    while True:
        c.tick(60)
        for events in pg.event.get():
            if events.type == pg.QUIT:
                pg.quit()
                return
        pg.display.update()

И в конце - функция main, которая, собственно и соберет все это в кучу и запустит:

if __name__ == '__main__':
    pow = 3
    root = 4
    m1 = pow_square(root, pow, m_4)
    try:
        to_star(m1)
    except pygame.error:
        print('Drawing finished.')
    print(m1)
    # Compute sums of all digits in rows, columns and diagonals.
    m2 = np.rot90(m1)
    print([sum(x) for x in m1])
    print([sum(x) for x in m2])
    print(sum([m1[x, x] for x in range(root ** pow)]))
    print(sum([m2[x, x] for x in range(root ** pow)]))

Быстрая проверка показала, что метод действительно работает. И при этом - не только на 3х3, но и на 4х4, 5х5, и, скорее всего, на всех других натуральных магических квадратах. Мой метод получения ассоциативных магических квадратов порядка 3n, 4n и 5n безукоризненно работает, и без проблем сработает с любым другим натуральным магическим квадратом (причем, не обязательно ассоциативным - все зависит только от заложенного шаблона)!

Несколько примеров работы программы

Оригинал магического квадрата порядка 81 (.csv)

Исходный код программы на GitHub

Примечание: программа работает довольно медленно с квадратами порядка > 1000. Не рекомендуется запускать на слабых ПК при больших размерностях. Подвисает во время отрисовки.

P.S. Скорее всего, код представленный в статье, имеет огромное количество точек оптимизации, и, вероятно, следовало обрабатывать матрицу многопоточно, но так как этот код практически не имеет иной области применения, кроме демонстрации метода генерации магических квадратов, он оставлен, как есть. Заранее извиняюсь за некоторое отсутствие эстетики в коде...

UPD1. Полевые испытания показали, что метод применим только к ассоциативным магическим квадратам. Проверка производилась только на нормальных магических квадратах. Код в статье обновлен в целях оптимизации. Были успешно сформированы магические квадраты 2048х2048 и 2401х2401.