👶🏽 👩‍👩‍👧 👊🏿 Dominos paving 🍌 😒 👼🏽

One of the first really interesting problems in mathematics that I encountered was formulated as follows: "two opposite corner cells were cut out from a chessboard, can the rest be cut into" dominoes "- figures of two cells with one side in common?" . It has a very simple formulation, unlike Fermat’s great theorem, it has a simple, elegant, but not obvious solution (if you know the solution to the problem, then try to apply it to the figure on the right).

In this article, I will talk about several algorithms that can cover an arbitrary cellular figure on a plane with dominoes, if possible, find situations where this is impossible and consider the number of possible dominoes tilings.

Nota bene! The material in this article is a truncated version of this jupyter-notebook , all the pictures and animations that you see in the article are generated by code from this laptop (though there will be no animations in the github preview). By the way, images from the header are also generated using python / matplotlib

Helper code for rendering, it will come in handy

import matplotlib.pyplot as plt
from matplotlib import colors as mcolors
colors = dict(mcolors.BASE_COLORS, **mcolors.CSS4_COLORS)

# Sort colors by hue, saturation, value and name.
by_hsv = sorted((tuple(mcolors.rgb_to_hsv(mcolors.to_rgba(color)[:3])), name)
                for name, color in colors.items())
names = [name for hsv, name in by_hsv if name not in {'black', 'k', 'w', 'white', 'crimson', 'royalblue', 'limegreen', 'yellow', 'orange'}]
import random
random.shuffle(names)
names = ['crimson', 'royalblue', 'limegreen', 'yellow', 'orange', *names]
names.append('red')
names.append('white')
names.append('black')

def fill_cell(i, j, color, ax):
    ax.fill([i, i, i + 1, i + 1, i], [j, j + 1, j + 1, j, j], color=color)

def draw_filling(filling):
    if filling is not None:
        n = len(filling)
        m = len(filling[0])
        fig = plt.figure(figsize=(m * 0.75, n * 0.75))

        ax = fig.add_axes([0, 0, 1, 1])
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)      

        for name, spine in ax.spines.items():
            spine.set_visible(False)
            spine.set_visible(False)

        for i, row in enumerate(filling):
            i = n - i - 1
            for j, cell in enumerate(row):
                fill_cell(j, i, names[cell], ax)

        for i in range(n + 1):
            ax.plot([0, m], [i, i], color='black')
        for i in range(m + 1):
            ax.plot([i, i], [0, n], color='black')
        plt.close(fig)
        return fig
    else:
        return None

Before moving on, solving the original problem

It is impossible to do this, and there is a beautiful and simple explanation for this:

30 32
,

Dynamic Programming by Profile

, , .

, - , . , ,

F_{n + 1} = F_{n} + F_{n - 1} .

$F_{n+1}=F_n+F_{n-1}.$

"" $C_n^k$ ,

C_{n}^{k} = \frac{n}{k} C_{n - 1}^{k - 1},

$C_n^k=\frac{n}{k}C_{n-1}^{k-1},$

C_{n}^{k} = C_{n - 1}^{k} + C_{n - 1}^{k - 1} .

$C_n^k=C_{n-1}^k+C_{n-1}^{k-1}.$

" $n\times m$ ?" . , , , .

: $(k-1)\times m$ , $k$ $j$ , , $k+1$ $j$ , . — $m$ : $i$ - , , $2^m$ ( — , — ).

$[0, 2^m)$ ( $m$ ).

$n\times m$ $k=n$ , $j=0$ .

$g(k, j, m, p)$ , $g(k, m, m, p)=g(k+1, 0, m, p)$ $g(k, j, m, p)$ $g(k, j-1, m, \cdot)$ , $j$ $j+1$ — : ,

,

. , , , — . , , .

( , ). ( — , — ), , -

tiling = [
    '........',
    '........',
    '........',
    '........',
    '........',
    '........',
    '........',
    '........',
]

def count_tilings(tiling):
    n = len(tiling)
    m = len(tiling[0])
    if ((n + 1) * m * 2 ** m) <= 10000000:
        dp = [[[(0 if k != 0 or j != 0 or mask != 0 else 1) for mask in range(2 ** m)] for j in range(m)] for k in range(n + 1)]
    for k in range(n):
        for j in range(m):
            for mask in range(2 ** m):
                #  
                if k < n - 1 and tiling[k][j] == '.' and tiling[k + 1][j] == '.' and (mask &  (1 << j)) == 0:
                    dp[k + ((j + 1) // m)][(j + 1) % m][mask + (1 << j)] += dp[k][j][mask]
                #  
                if j < m - 1 and tiling[k][j] == '.' and tiling[k][j + 1] == '.' and (mask & (3 << j)) == 0:
                    dp[k + ((j + 1) // m)][(j + 1) % m][mask + (2 << j)] += dp[k][j][mask]
                #  
                if ((1 << j) &  mask) != 0 or tiling[k][j] != '.':
                    dp[k + ((j + 1) // m)][(j + 1) % m][(mask | (1 << j)) - (1 << j)] += dp[k][j][mask]
    return dp

dp = count_tilings(tiling)
print(dp[8][0][0])

12988816

, . , $n\times 2$ — .

tiling_fib = [
    '..',
    '..',
    '..',
    '..',
    '..',
    '..',
    '..',
    '..'
]
dp = count_tilings(tiling_fib)
for i in range(8):
    print(dp[i][0][0], end=' ')

1 1 2 3 5 8 13 21

, ,

tiling_no_corners_opposite = [
    '.......#',
    '........',
    '........',
    '........',
    '........',
    '........',
    '........',
    '#.......',
]
dp = count_tilings(tiling_no_corners_opposite)
print(dp[8][0][0])

. , ?

def cover_if_possible(tiling, dp=None):
    if dp is None:
        dp = count_tilings(tiling)
    n = len(dp) - 1
    m = len(dp[0])
    if dp[n][0][0] == 0:
        return None

    result = [[-1 if tiling[i][j] == '#' else 0 for j in range(m)] for i in range(n)]
    num = 0

    k = n
    j = 0
    mask = 0

    while k > 0 or j > 0:
        #print(k, j, mask)
        prev_j = j - 1
        prev_k = k
        if prev_j == -1:
            prev_j += m
            prev_k -= 1

        #   ,       i, j, mask
        #     ,         
        #   
        for prev_mask in range(2 ** m):
            if prev_k < n - 1 and tiling[prev_k][prev_j] == '.' and tiling[prev_k + 1][prev_j] == '.' and \
                (prev_mask & (1 << prev_j)) == 0 and (prev_mask + (1 << prev_j)) == mask and dp[prev_k][prev_j][prev_mask] != 0:
                mask = prev_mask
                result[prev_k][prev_j] = num
                result[prev_k + 1][prev_j] = num
                #print(f'Vertical at ({prev_k}, {prev_j}) ({prev_k + 1}, {prev_j})')
                num += 1
                break
            elif prev_j < m - 1 and tiling[prev_k][prev_j] == '.' and tiling[prev_k][prev_j + 1] == '.' and (prev_mask & (3 << prev_j)) == 0 and \
                prev_mask + (2 << prev_j) == mask and dp[prev_k][prev_j][prev_mask] != 0:
                mask = prev_mask
                result[prev_k][prev_j] = num
                result[prev_k][prev_j + 1] = num
                #print(f'Horisontal at ({prev_k}, {prev_j}) ({prev_k}, {prev_j + 1})')
                num += 1
                break
            elif (((1 << prev_j) & prev_mask) != 0 or tiling[prev_k][prev_j] != '.') and \
                (prev_mask | (1 << prev_j)) - (1 << prev_j) == mask and dp[prev_k][prev_j][prev_mask] != 0:
                mask = prev_mask
                break

        j = prev_j
        k = prev_k

    return result

filling = cover_if_possible(tiling)
draw_filling(filling)

tiling_random = [
    '........',
    '#.#.....',
    '..#.....',
    '........',
    '........',
    '........',
    '........',
    '...#....'
]
filling_random = cover_if_possible(tiling_random)
draw_filling(filling_random)

def maxmimum_cover(tiling):
    n = len(tiling)
    m = len(tiling[0])
    if ((n + 1) * m * 2 ** m) <= 10000000:
        dp = [[[(n * m if k != 0 or j != 0 or mask != 0 else 0) for mask in range(2 ** m)] for j in range(m)] for k in range(n + 1)]
    for k in range(n):
        for j in range(m):
            for mask in range(2 ** m):
                next_k, next_j = k + ((j + 1) // m), (j + 1) % m
                #  
                if k < n - 1 and tiling[k][j] == '.' and tiling[k + 1][j] == '.' and (mask &  (1 << j)) == 0:
                    dp[next_k][next_j][mask + (1 << j)] = min(dp[next_k][next_j][mask + (1 << j)], dp[k][j][mask])
                #  
                if j < m - 1 and tiling[k][j] == '.' and tiling[k][j + 1] == '.' and (mask & (3 << j)) == 0:
                    dp[next_k][next_j][mask + (2 << j)] = min(dp[next_k][next_j][mask + (2 << j)], dp[k][j][mask])
                #  
                if ((1 << j) &  mask) != 0 or tiling[k][j] != '.':
                    dp[next_k][next_j][(mask | (1 << j)) - (1 << j)] = \
                        min(dp[next_k][next_j][(mask | (1 << j)) - (1 << j)], dp[k][j][mask])
                #   ,    
                else:
                    dp[next_k][next_j][(mask | (1 << j)) - (1 << j)] = \
                        min(dp[next_k][next_j][(mask | (1 << j)) - (1 << j)], dp[k][j][mask] + 1)
    return dp

def cover_maximum_possible(tiling, dp=None):
    if dp is None:
        dp = maxmimum_cover(tiling)
    n = len(dp) - 1
    m = len(dp[0])

    result = [[-1 if tiling[i][j] == '#' else -2 for j in range(m)] for i in range(n)]
    num = 0

    k = n
    j = 0
    mask = 0

    while k > 0 or j > 0:
        #print(k, j, mask)
        prev_j = j - 1
        prev_k = k
        if prev_j == -1:
            prev_j += m
            prev_k -= 1

        #   ,       i, j, mask
        #     ,         
        #   
        for prev_mask in range(2 ** m):
            #    ,        0,  
            # ,       
            if prev_k < n - 1 and tiling[prev_k][prev_j] == '.' and tiling[prev_k + 1][prev_j] == '.' and \
                    (prev_mask & (1 << prev_j)) == 0 and (prev_mask + (1 << prev_j)) == mask and \
                    dp[prev_k][prev_j][prev_mask] == dp[k][j][mask]:
                mask = prev_mask
                result[prev_k][prev_j] = num
                result[prev_k + 1][prev_j] = num
                num += 1
                break
            elif prev_j < m - 1 and tiling[prev_k][prev_j] == '.' and tiling[prev_k][prev_j + 1] == '.' and (prev_mask & (3 << prev_j)) == 0 and \
                    prev_mask + (2 << prev_j) == mask and dp[prev_k][prev_j][prev_mask] == dp[k][j][mask]:
                mask = prev_mask
                result[prev_k][prev_j] = num
                result[prev_k][prev_j + 1] = num
                num += 1
                break
            elif (((1 << prev_j) & prev_mask) != 0 or tiling[prev_k][prev_j] != '.') and \
                    (prev_mask | (1 << prev_j)) - (1 << prev_j) == mask and dp[prev_k][prev_j][prev_mask] == dp[k][j][mask]:
                mask = prev_mask
                break
            elif ((1 << prev_j) & prev_mask) == 0 and tiling[prev_k][prev_j] == '.' and \
                    (prev_mask | (1 << prev_j)) - (1 << prev_j) == mask and dp[prev_k][prev_j][prev_mask] + 1 == dp[k][j][mask]:
                mask = prev_mask
                break

        j = prev_j
        k = prev_k

    return result

tiling_custom=[
    '...####',
    '....###',
    '......#',
    '#.#....',
    '#......',
    '##.....',
    '###...#',
]
filling = cover_maximum_possible(tiling_custom)
draw_filling(filling)

! , . , . $n\times m\times 2^m$ , , $\mathcal{O}(nm2^m)$ , $\mathcal{O}(nm2^{\min(n, m)})$ , $n<m$ . $\mathcal{O}(nm2^m)$ , , $\mathcal{O}(nm)$ , .

- . , — , , . , — . , , — , — , — . , , , . , . , , .

def check_valid(i, j, n, m, tiling):
    return 0 <= i and i < n and 0 <= j and j < m and tiling[i][j] != '#'

def find_augmenting_path(x, y, n, m, visited, matched, tiling):
    if not check_valid(x, y, n, m, tiling):
        return False
    if (x, y) in visited:
        return False
    visited.add((x, y))

    for dx, dy in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
        if not check_valid(x + dx, y + dy, n, m, tiling):
            continue
        if (x + dx, y + dy) not in matched or find_augmenting_path(*matched[(x + dx , y + dy)], n, m, visited, matched, tiling):
            matched[(x + dx, y + dy)] = (x, y)
            return True
    return False

def convert_match(matched, tiling, n, m):
    result = [[-1 if tiling[i][j] == '#' else -2 for j in range(m)] for i in range(n)]
    num = 0
    for x, y in matched:
        _x, _y = matched[(x, y)]
        result[x][y] = num
        result[_x][_y] = num
        num += 1
    return result

def match_with_flow(tiling):
    result_slices = []
    n = len(tiling)
    m = len(tiling[0])

    matched = dict()
    #   
    rows = list(range(n))
    columns = list(range(m))
    random.shuffle(rows)
    random.shuffle(columns)
    result_slices.append(convert_match(matched, tiling, n, m))

    for i in rows:
        for j in columns:
            if (i + j) % 2 == 1:
                continue
            visited = set()
            if find_augmenting_path(i, j, n, m, visited, matched, tiling):
                result_slices.append(convert_match(matched, tiling, n, m))

    return result_slices

sequencial_match = match_with_flow(tiling_custom)

( ) " ". , , , . , , , , , .

UPD. - . , :

( )
, .

, .

! 5 .

. , . , . : , , , . . 4 , . 5 , , ( 5 , )

Painting dominoes in 5 colors

def color_5(filling):
    result = [[i for i in row] for row in filling]
    #  
    domino_tiles = [[] for i in range(max(map(max, filling)) + 1)]
    domino_neighbours = [set() for i in range(max(map(max, filling)) + 1)]
    degree = [0 for i in range(max(map(max, filling)) + 1)]

    n = len(filling)
    m = len(filling[0])

    for i, row in enumerate(filling):
        for j, num in enumerate(row):
            if num >= 0:
                domino_tiles[num].append((i, j))

    for i, tiles in enumerate(domino_tiles):
        for x, y in tiles:
            for dx, dy in [(-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (-1, 1), (1, -1), (1, 1)]:
                a, b = x + dx, y + dy
                if 0 <= a and a < n and 0 <= b and b < m and filling[a][b] >= 0 and filling[a][b] != i \
                        and filling[a][b] not in domino_neighbours[i]:
                    domino_neighbours[i].add(filling[a][b])
                    degree[i] += 1

    #       ,      5  
    # .     ,   ,     
    #           ,     
    active_degrees = [set() for i in range(max(degree) + 1)]
    for i, deg in enumerate(degree):
        active_degrees[deg].add(i)

    reversed_order = []

    for step in range(len(domino_tiles)):
        min_degree = min([i for i, dominoes in enumerate(active_degrees) if len(dominoes) > 0])
        domino = active_degrees[min_degree].pop()

        reversed_order.append(domino)
        for other in domino_neighbours[domino]:
            if other in active_degrees[degree[other]]:
                active_degrees[degree[other]].remove(other)
                degree[other] -= 1
                active_degrees[degree[other]].add(other)                

    #             ,
    #     5 ,      ,
    #    .     
    colors = [-1 for domino in domino_tiles]
    slices = [draw_filling(result)]
    for domino in reversed(reversed_order):
        used_colors = [colors[other] for other in domino_neighbours[domino] if colors[other] != -1]

        domino_color = len(used_colors)
        for i, color in enumerate(sorted(set(used_colors))):
            if i != color:
                domino_color = i
                break

        if domino_color < 5:
            colors[domino] = domino_color
            for x, y in domino_tiles[domino]:
                result[x][y] = domino_color

            slices.append(draw_filling(result))        
            continue

        #      ,           
        c = 0
        other = [other for other in domino_neighbours[domino] if colors[other] == c]
        visited = set([other])
        q = Queue()
        q.put(other)
        domino_was_reached = False
        while not q.empty():
            cur = q.get()
            for other in domino_neighbours[cur]:
                if other == domino:
                    domino_was_reached = True
                    break
                if color[other] == c or color[other] == c + 1 and other not in visited:
                    visited.add(other)
                    q.put(other)

        if not domino_was_reached:
            for other in visited:
                color[other] = color[other] ^ 1
                for x, y in domino_tiles[other]:
                    result[x][y] = color[other]
            color[domino] = c
            for x, y in domino_tiles[domino]:
                result[x][y] = c

            slices.append(draw_filling(result))
            continue

        #     2  3.
        c = 2
        other = [other for other in domino_neighbours[domino] if colors[other] == c]
        visited = set([other])
        q = Queue()
        q.put(other)
        domino_was_reached = False
        while not q.empty():
            cur = q.get()
            for other in domino_neighbours[cur]:
                if other == domino:
                    domino_was_reached = True
                    break
                if color[other] == c or color[other] == c + 1 and other not in visited:
                    visited.add(other)
                    q.put(other)
        for other in visited:
            color[other] = color[other] ^ 1
            for x, y in domino_tiles[other]:
                result[x][y] = color[other]
        color[domino] = c
        for x, y in domino_tiles[domino]:
            result[x][y] = c
        slices.append(draw_filling(result))

    return result, slices

If you are going to use this code, then note that each step is drawn there - this was necessary for the animation, which greatly slows down the algorithm. If you need only the final painting, then remove all the code that uses the slices variable.

Well, finally, try one of the examples that was a little earlier

Dominos paving

Dynamic Programming by Profile

! 5 .

More articles: