работа над ошибками

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

import flint

def count_tilings(n,m,ext):
    ans = flint.fmpz_mat(2**n,1)
    ans[0,0] = 1
    for _ in range(m):
        ans = ext * ans
    return ans[0,0]

def newext(n,ext):
    ans = flint.fmpz_mat(2**n,2**n)
    for mask in range(2**n):
        for perp in range(2**n):
            if mask&perp == 0:
                if (mask+perp)%2 == 1:
                    ans[mask,perp] = ext[n-1][mask>>1,perp>>1]
                elif (mask+perp)%4 == 0:
                    ans[mask,perp] = ext[n-2][mask>>2,perp>>2]
    return ans

N = 5
ext = [flint.fmpz_mat(1,1,[1]),flint.fmpz_mat(2,2,[0,1,1,0])]
for n in range(2,2*N+1):
    ext.append(newext(n,ext))
    if n%2 == 0:
        ans = count_tilings(n,n,ext[n])
        print(n,ans)

12988816 разбиений квадрата 8×8 теперь подсчитываются за пару сотых секунды



теперь легко проверить, что разные делимости есть и для разбиений прямоугольников на доминошки — например, для 6×2 разбиений 13 => для (7k+6)×(3l+2) количество разбиений делится на 13 (и действительно, для 13×20 получается 5422065916132126528319352874496 разбиений, это число делится на 13)

количество разбиений прямоугольника 2×N — это просто число Фибоначчи, и что F_N | F_{2N+1} легко доказать комбинаторно (посмотреть как там покрыт доминошками средний столбец)

а как доказать комбинаторно аналогичное утверждения для разбиений прямоугольника, say, 6×N — я не знаю (само утверждение тоже прочитал у Проппа)

***

еще переписал генерацию случайного разбиения на доминошки ацтекского брильянта — через domino shuffling (aka urban renewal)

полный код не влезает в лимит телеграма, так что положу только основную функцию, которая делает из случайного брильянта порядка size случайный брильянт порядка size+1… point, наверное, в том, что получилось не сложнее, чем наивный код с флипами (но, конечно, флипы работают для разных областей, а не только для брильянтов)

def extend(field,size):
    newfield = [ [' ']*(2*size+2) for _ in range(2*size+2) ]
    for i in range(2*size+2):
        for j in range(2*size+2):
            if inside(i,j,size+1) and newfield[i][j]==" ":
                if inside(i,j-1,size) and field[i][j-1] in "Uu" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Dd"):
                    newfield[i][j] = field[i][j-1]
                elif inside(i-2,j-1,size) and field[i-2][j-1] in "Dd" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Uu"):
                    newfield[i][j] = field[i-2][j-1]
                elif inside(i-1,j,size) and field[i-1][j] in "Ll" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Rr"):
                    newfield[i][j] = field[i-1][j] 
                elif inside(i-1,j-2,size) and field[i-1][j-2] in "Rr" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Ll"):
                    newfield[i][j] = field[i-1][j-2]
                else:
                    bit = getrandbits(1)
                    newfield[i][j], newfield[i][j+1] = ("U", "u") if bit else ("L", "R")
                    newfield[i+1][j], newfield[i+1][j+1] = ("D", "d") if bit else ("l", "r")
    return newfield