Project Euler Problem 594

Solution

Answer: 47067598

Let

$$O_{a,b}$$

be the equal-angled convex octagon whose side lengths alternate between $a$ and $b$.

We must compute

$$t(O_{4,2}),$$

the number of tilings by unit squares and unit rhombi (with angles $45^\circ/135^\circ$).

A standard way to attack these tiling problems is to convert them into families of non-intersecting lattice paths and then apply the Lindström–Gessel–Viennot (LGV) determinant theorem.

---

1. Mathematical analysis

The octagon

$$(a,b,c,d,a,b,c,d)$$

with equal angles $135^\circ$ admits a de Bruijn line representation: every tiling corresponds uniquely to a collection of monotone paths.

For this family, Elnitsky’s determinantal enumeration formula gives:

$$T_{a,b,c,d} = \sum_{X,Y} \left( \prod_{u=1}^{d+1}\det M^{(u)} \right) \left( \prod_{v=1}^{b+1}\det P^{(v)} \right),$$

where:

• $X=(x_{k,\ell})$ is a monotone $b\times d$ matrix with entries in $[0,a]$,
• $Y=(y_{k,\ell})$ is a reversed monotone $b\times d$ matrix with entries in $[0,c]$,
• the determinants count non-intersecting path families via LGV.

For the Project Euler polygon,

$$O_{a,b}=(a,b,a,b,a,b,a,b),$$

so we set

$$(a,b,c,d)=(a,b,a,b).$$

The matrices are:

$$M^{(u)}{ij} = \binom{ (x{j,u}-x_{i,u-1})+(y_{j,u}-y_{i,u-1}) }{ (x_{j,u}-x_{i,u-1})+(j-i) },$$

and

$$P^{(v)}{ij} = \binom{ (x{v,j}-x_{v-1,i})+(y_{v-1,i}-y_{v,j}) }{ (x_{v,j}-x_{v-1,i})+(j-i) }.$$

The total number of tilings is obtained by summing the products of these determinants over all admissible monotone arrays $X,Y$.

---

2. Python implementation

from math import comb

def gen_monotone_matrices(rows, cols, max_value):
    """
    Generate all matrices M such that:
        M[i][j] >= M[i-1][j]
        M[i][j] >= M[i][j-1]
    with entries in [0, max_value].
    """
    m = [[0] * cols for _ in range(rows)]

    def rec(pos):
        if pos == rows * cols:
            yield tuple(tuple(r) for r in m)
            return

        i, j = divmod(pos, cols)

        lo = 0
        if i > 0:
            lo = max(lo, m[i - 1][j])
        if j > 0:
            lo = max(lo, m[i][j - 1])

        for v in range(lo, max_value + 1):
            m[i][j] = v
            yield from rec(pos + 1)

    yield from rec(0)

def det_bareiss(mat):
    """
    Exact determinant using the Bareiss algorithm.
    """
    n = len(mat)

    if n == 0:
        return 1

    a = [row[:] for row in mat]
    sign = 1
    prev = 1

    for k in range(n - 1):

        if a[k][k] == 0:
            for i in range(k + 1, n):
                if a[i][k]:
                    a[k], a[i] = a[i], a[k]
                    sign = -sign
                    break
            else:
                return 0

        pivot = a[k][k]

        for i in range(k + 1, n):
            for j in range(k + 1, n):
                a[i][j] = (
                    a[i][j] * pivot
                    - a[i][k] * a[k][j]
                ) // prev

        prev = pivot

    return sign * a[-1][-1]

def count_tilings_octagon(a, b, c, d):

    X = list(gen_monotone_matrices(b, d, a))

    Y = []
    for y_rev in gen_monotone_matrices(b, d, c):
        y = tuple(
            tuple(y_rev[b - 1 - i][j] for j in range(d))
            for i in range(b)
        )
        Y.append(y)

    total = 0

    for x in X:

        x_full = [[0] * (d + 2) for _ in range(b + 2)]

        for k in range(1, b + 1):
            for l in range(1, d + 1):
                x_full[k][l] = x[k - 1][l - 1]

        for k in range(1, b + 1):
            x_full[k][d + 1] = a

        for l in range(d + 2):
            x_full[b + 1][l] = a

        for y in Y:

            y_full = [[0] * (d + 2) for _ in range(b + 2)]

            for k in range(1, b + 1):
                for l in range(1, d + 1):
                    y_full[k][l] = y[k - 1][l - 1]

            for k in range(1, b + 1):
                y_full[k][d + 1] = c

            for l in range(d + 2):
                y_full[0][l] = c

            prod = 1

            # Product of det(M^(u))
            for u in range(1, d + 2):

                M = []

                for i in range(1, b + 1):

                    row = []

                    for j in range(1, b + 1):

                        A = (
                            x_full[j][u] - x_full[i][u - 1]
                            + y_full[j][u] - y_full[i][u - 1]
                        )

                        B = (
                            x_full[j][u] - x_full[i][u - 1]
                            + (j - i)
                        )

                        row.append(comb(A, B) if 0 <= B <= A else 0)

                    M.append(row)

                prod *= det_bareiss(M)

                if prod == 0:
                    break

            if prod == 0:
                continue

            # Product of det(P^(v))
            for v in range(1, b + 2):

                P = []

                for i in range(1, d + 1):

                    row = []

                    for j in range(1, d + 1):

                        A = (
                            x_full[v][j] - x_full[v - 1][i]
                            + y_full[v - 1][i] - y_full[v][j]
                        )

                        B = (
                            x_full[v][j] - x_full[v - 1][i]
                            + (j - i)
                        )

                        row.append(comb(A, B) if 0 <= B <= A else 0)

                    P.append(row)

                prod *= det_bareiss(P)

                if prod == 0:
                    break

            total += prod

    return total

def t_O_ab(a, b):
    return count_tilings_octagon(a, b, a, b)

# Verify given examples
print(t_O_ab(1, 1))  # 8
print(t_O_ab(2, 1))  # 76
print(t_O_ab(3, 2))  # 456572

# Required value
print(t_O_ab(4, 2))

---

3. Code walkthrough

• gen_monotone_matrices(...)

generates all monotone arrays $X$ and $Y$.

• det_bareiss(...)

computes determinants exactly using the fraction-free Bareiss algorithm.

• count_tilings_octagon(...)

implements the determinantal formula directly.

• t_O_ab(a,b)

specializes the general octagon to the Project Euler shape.

The verification outputs are:

$$t(O_{1,1})=8, \qquad t(O_{2,1})=76, \qquad t(O_{3,2})=456572,$$

matching the problem statement exactly.

Finally, the computation gives:

$$t(O_{4,2}) = 47067598.$$

---

4. Final verification

The answer is:

• positive,
• substantially larger than $t(O_{3,2})=456572$,
• consistent with the rapid combinatorial growth of tiling counts.

The implementation reproduces all supplied test values exactly, strongly validating the computation.

Answer: 47067598