Square triangulations

,

Project Euler problem 270 asks:

A square piece of paper with integer dimensions \(n×n\) is placed with a corner at the origin and two of its sides along the \(x\)- and \(y\)-axes. Then, we cut it up respecting the following rules:

Counting any reflections or rotations as distinct, we call \(C(n)\) the number of ways to cut an \(n×n\) square. What is \(C(30) \bmod 10^8\)?

Spoilers below.

The question asks us to compute the number of ways to triangulate an oriented \(n×n\) square using only lines joining the \(4n\) integer points around the edge of the square. For example, here are four ways in which the \(3×3\) square can be triangulated according to these rules:

To solve this, we need to find a recurrence relation, ideally one expressing \(C(n)\) in terms of \(C(m)\) for \(m < n\). But a quick look at the shapes that you can get after removing a few triangles shows that this is not going to be as simple as that. These smaller shapes in the recurrence can have more sides than the original, and sides of different lengths at that:

So I am going to describe a convex polygon by a list of integers giving the number of interior vertices on each side. So the 3×3 square is described by (2, 2, 2, 2), and the five-sided figure on the right is described by (0, 1, 2, 2, 1):

Why count the interior vertices on each side, and not the edges, or all the vertices? Well, the choice is arbitrary, but in the representation I’ve chosen the exceptional cases arise when a side has no interior vertices, and 0 feels more comfortable as an exceptional value than 1 or 2.

Now to work out the recurrence. If you read my post on counting the triangulations of a convex polygon, then you’ll recall that one of the difficulties in developing a recurrence is avoiding double counting, and my solution there was to pick a distinguished edge and consider the unique triangle in the triangulation containing this edge. So let’s try the same thing here and consider the cases.

  1. The edge might be removed as part of an ear, leaving the remainder in one piece. This can happen in six ways, depending on the length of the side containing the distinguished edge, and on the lengths of the adjacent sides.

    Each of these diagrams shows an illustrative example that stands for a general case. For example, diagram A(i) stands for the general case \(a, \ldots, z → 0, a − 1, \ldots, z − 1\).

  2. The edge might be removed as part of a triangle whose third vertex is an interior vertex on one of the other sides of the figure, leaving the remainder in two pieces. This can happen in two ways, depending on whether the edge is the only edge on its side or not:

  3. The edge might be removed as part of a triangle whose third vertex is a corner vertex, leaving the remainder in two pieces. This can happen in two ways, depending on whether the edge is the only edge on its side or not:

    Note that case C(ii) can only happen when there are at least five sides in the figure, and so it would have been easy to miss!

Ten cases! Let’s code it up, using the @memoized decorator from my earlier post:

# Base case: exactly one way to triangulate a triangle.
@memoized({(0, 0, 0): 1})
def C(c):
    """Return the number of ways to triangulate an oriented convex
    polygon, using only lines joining points on the edges of the
    polygon. The tuple c gives the number of interior points on each
    side of the polygon.

    """
    assert(len(c) >= 3)

    # Number of triangulations found.
    t = 0

    # Case A: distinguished edge removed as part of an ear.
    if c[0] > 0 and c[-1] > 0:                     # case A(i)
        t += C((0, c[0] - 1) + c[1:-1] + (c[-1] - 1,))
    elif c[0] > 0 and c[-1] == 0:                  # case A(ii)
        t += C((0, c[0] - 1) + c[1:-1])
    elif c[0] == 0 and c[-1] > 0:                  # case A(iv)
        t += C((0,) + c[1:-1] + (c[-1] - 1,))
    elif c[0] == 0 and c[-1] == 0 and len(c) >= 4: # case A(vi)
        t += C((0,) + c[1:-1])
    if c[0] == 0 and c[1] > 0:                     # case A(iii)
        t += C((0, c[1] - 1) + c[2:])
    elif c[0] == 0 and c[1] == 0 and len(c) >= 4:  # case A(v)
        t += C((0,) + c[2:])

    # Case B: distinguished edge removed as part of triangle whose
    # third vertex is the j'th interior vertex on side i.
    for i in range(1, len(c) - 1):
        for j in range(c[i]):
            u = C((0, c[i] - j - 1) + c[i + 1:])
            if c[0] > 0:                           # case B(i)
                t += u * C((0,) + (c[0] - 1,) + c[1:i] + (j,))
            elif c[0] == 0 and i > 1:              # case B(ii)
                t += u * C((0,) + c[1:i] + (j,))

    # Case C: distinguished edge removed as part of triangle whose
    # third vertex is the corner between sides i and i + 1.
    for i in range(1, len(c) - 2):
        u = C((0,) + c[i+1:])
        if c[0] > 0:                               # case C(i)
            t += u * C((0,) + (c[0] - 1,) + c[1:i+1])
        elif c[0] == 0 and i > 1:                  # case C(ii)
            t += u * C((0,) + c[1:i+1])

    return t

When you have code with a lot of cases, it’s hard to be sure that you haven’t made a mistake in one or two of them. What tests can we try? Well, there are a couple of sample values in the Project Euler problem:

>>> C((0,) * 4)
2
>>> C((1,) * 4)
30

Also, we already have a program that counts the number of triangulations of a convex polygon, so we can compare with that:

>>> all(T(i) == C((0,) * i) for i in range(3, 100))
True

And the final test is whether it runs quickly enough to get the answer to the Project Euler problem in a reasonable amount of time:

>>> timeit(lambda:C((29,) * 4), number=1)
0.2441568149952218

An interesting wrinkle in this problem is that there is a bit of flexibility in how to implement the recurrence, because each figure has multiple representations. For example, when we cut an ear off the square (2, 2, 2, 2) we might represent the resulting pentagon as (0, 1, 2, 2, 1) or (1, 2, 2, 1, 0) or (2, 2, 1, 0, 1) or (2, 1, 0, 1, 2) or (1, 0, 1, 2, 2). Any of these would produce the right answer. But some of them would take longer to do so, because they would visit a wider range of figures along the way. In the implementation above, I’ve arranged always to recurse onto figures whose representation starts with a zero. This means that the recursion tends to visit figures with shorter representations, because the side represented by the zero comes first and so gets completely removed in cases A(v), A(vi), B(ii) and C(ii).

As a concrete example, consider the computation of \(C(1,1,1,1) = 30\). There are, I believe, 28 possible figures1 that might need to be visited in the course of this computation, including one hexagon (that you’d get if you started by cutting off diagonally opposite ears). But the implementation above visits only nine of them, a fact we can verify by inspecting the table used by the @memoized decorator to cache the arguments and corresponding results:2

>>> C((1,) * 4)
30
>>> from inspect import getclosurevars
>>> sorted(getclosurevars(C).nonlocals['table'])
[(0, 0, 0), (0, 0, 1), (0, 0, 1, 0), (0, 0, 1, 1), (0, 0, 1, 1, 0), (0, 1, 0),
 (0, 1, 1), (0, 1, 1, 0), (1, 1, 1, 1)]

If I rewrite the implementation so that it sometimes puts the zero at the beginning and sometimes at the end, the computation needs to visit fifteen figures:

>>> sorted(getclosurevars(D).nonlocals['table'].keys())
[(0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 0, 0),
 (0, 1, 0, 0, 0), (0, 1, 1), (0, 1, 1, 0), (0, 1, 1, 0, 0), (1, 0, 0),
 (1, 0, 0, 0), (1, 1, 0), (1, 1, 0, 0), (1, 1, 1, 1)]

And that inefficiency makes a big difference to the computation time for larger examples:

>>> timeit(lambda:D((29,) * 4), number=1)
0.5947780360002071

Update 2013-08-09. In this article I show how to revise and extend this code to draw the triangulations.


  1.  000, 001, 010, 100, 011, 101, 110, 0000, 0001, 0010, 0100, 1000, 0011, 0110, 1100, 1001, 00000, 00011, 00110, 01100, 11000, 10001, 00001, 00010, 00100, 01000, 10000, and 000000.

  2.  You’ll need Python 3.3 in order to use inspect.getclosurevars.