Two algorithms for randomly generating aperiodic tilings

[Simon Tatham, 2023-04-10]

• Introduction
• Substitution systems
  • Penrose tilings: half-tile triangles
  • Hat tiling: a system of four metatiles
• Choosing a random patch from a fixed expansion
• Combinatorial coordinates
  • Penrose tilings: indexing triangle edges
  • Hats: indexing kite edges, and dealing with non-uniqueness
  • Making it up as you go along
  • Generating the output tiling
  • Advantages
• Conclusion
• Appendix: exact coordinate systems
  • Complex roots of unity
  • Separate x and y coordinates using √5
• Appendix: interleaving the Penrose tiling types
• References

Introduction

In the 1970s, Roger Penrose discovered several sets of polygons
which will tile the plane, but only aperiodically,
without the tiling repeating in a fixed pattern.

Samples of his two best-known tiling types, called P2 and P3,
are shown below. P2 is composed of quadrilaterals usually
referred to as ‘kite’ and ‘dart’; P3 is composed of two sizes of
rhombus, known as ‘rhombs’ in this context for some reason:

(Strictly speaking, any one of these polygons can tile
the plane periodically, because so can any quadrilateral at all.
The tiling is only forced to be aperiodic if you restrict the
ways the tiles can connect. You could do this by putting
jigsaw-piece protrusions and holes on the edges, but that’s
ugly, so people normally prefer to show the plain quadrilaterals
and let the matching rules be implicit. In this article I’ll
ignore this detail just like everyone else.)

In 2023, a long-open question was answered: can you force
aperiodicity in the same way using only one type of
tile, instead of Penrose’s two? You can! The shape in question
is described as a ‘hat’ by its discoverers Smith, Myers, Kaplan
and Goodman-Strauss, and a sample of its tiling is shown
below:

(There’s a quibble with this one too: the hat shape is
asymmetric, and in order to tile the whole plane, you need to
use it in both handednesses. For some purposes this still means
you need two types of tile! But from a mathematical perspective
it’s more natural to regard reflections as not a fundamentally
different shape.)

One of the games in my
puzzle
collection, namely “Loopy” (also known as Slitherlink), can
be played on a lot of different tilings of the plane. Most of
them are periodic, but all three of these aperiodic tilings are
also perfectly feasible to play Loopy on. So the question
arises: how do you best get a computer to generate pieces of
these tilings?

One really obvious option would be to pre-generate a fixed
patch of each tiling, and set Loopy puzzles on that by just
varying the puzzle solution (the closed loop you’re trying to
reconstruct from clues). But that’s boring – surely the whole
point of a tiling that does not repeat is that it
shouldn’t be the same every time! So when you request a new
game, you should get a different patch of the tiling
itself, as well as a different loop to find in it.

So our mission goal is: given an output region and a source
of random numbers, generate a randomly selected patch of each
of these three aperiodic tilings.

In this article I’ll discuss two very different algorithms for
doing this, and say which one I prefer. (Spoiler: there will be
a very definite one that I prefer.)

Substitution systems

All three of the tilings we’ll discuss here can be generated
recursively by the same basic method.

For each of these tilings, there’s a system for starting with
an existing tiling, and then substituting each tile for smaller
tiles according to a fixed rule. This generates a new tiling,
using the same set of tiles, but more of them at a smaller
scale.

Because the output tiling is of the same type, it follows that
you can substitute that in the same way, and so on. So if you
repeat this many times, you can start from a single initial tile,
and expand it into as large a region of the plane as you like.

One slight complication, in every case, is that the most
convenient substitution system is not based on the actual tiles
you want to end up with. In each case, you have to substitute
a different system of tiles, and after you’ve done that
enough, convert the result into the final tiles you want as
output.

Penrose tilings: half-tile triangles

The substitution systems for the Penrose P2 and P3 tilings work
in very similar ways, so I’ll describe both of them
together.

In both cases, the most convenient way to think of the
substitution system is to imagine each of the starting tiles
being divided in half, to make two isosceles triangles. That
way, the substitution can replace each of those triangles with
either two or three smaller triangles, exactly filling the same
area as the original triangle. (If you didn’t do this, then most
smaller tiles would half-overlap two larger tiles, and it would
be more awkward to keep track of which tiles you’d expanded from
where.)

The divisions for each tile type look like this. For the P2
tiling, the kite and dart are each divided down their line of
symmetry. For the P3 tiling, the thin rhomb is cut along its
short diagonal, and the thick one along its long diagonal:

(Therefore, in the P2 tiling, the acute isosceles triangles are
larger than the obtuse ones, whereas in the P3 tiling, vice
versa.)

The two halves of each original tile are shown in different
colours, because they will be treated differently by the
substitution rules. Specifically, they will be mirror images of
each other.

Here are the substitution rules for each of those types of
triangle:

To see the complete system in action, select the radio buttons
below to see the effect of running a different number of
iterations, starting from a single obtuse triangle of each
tiling. In each step, you can see that every large triangle
transforms into smaller ones in accordance with the rules
above.

P2:
0
1
2
3
4
5
6
7
8

P3:
0
1
2
3
4
5
6
7
8

So one way to generate Penrose tilings is to expand a starting
triangle like this, and when you’ve expanded it enough times,
glue each pair of adjacent half-tiles back together into a
single kite, dart or rhomb.

Hat tiling: a system of four metatiles

The hat tiling can also be generated by a substitution system.
In this case, however, the system is based on a set of four
tiles that have no obvious relation to the actual hat shape!

The paper refers to these as “metatiles”, to avoid confusion
with the hat itself (the single tile). The four
metatiles are shown below, with their substitution rules:

The arrows on three of the metatile types indicate orientation:
those three metatiles are symmetric, but in the final conversion
to hats, each one will generate an entirely asymmetric pattern
of hats. So the substitution process has to remember what every
metatile’s orientation should be, in order not to get the last
step wrong.

(The remaining metatile is already asymmetric, so it doesn’t
need an arrow.)

In these substitution diagrams, the area of the original large
metatile is completely covered by the smaller ones, and the
smaller ones also overlap the edges. This means that when two
metatiles are placed next to each other, their expansions will
overlap too. In fact, they will always overlap in a way that
causes some of the resulting metatiles to coincide.

Here’s an illustration of the expansion in action, expanding a
single starting metatile for a few levels of recursion:

0
1
2
3

So, similarly to the previous section, if you want to generate
a large patch of hat tiling, you can start with a single
metatile and expand it repeatedly until you think it’s big
enough. Finally, each metatile is converted into either 1, 2 or
4 hats according to the following rules:

(Incidentally, the hat in the centre of the cluster of four
expanded from the hexagonal metatile, shown here in the darkest
red colour, is the only one of these hats that is reflected
compared to the rest. So you can readily identify the rare
reflected hats in this tiling, if you’re generating it using
this metatile system.)

However, beware! The final expansion step isn’t quite as easy
as it looks, because the expansion into hats
is distorting.

As I’ve drawn the metatiles here, the substitution of one set
of metatiles for a smaller set is geometrically precise: each
large metatile from one iteration of the algorithm can be
overlaid on the next iteration and the vertices will align
exactly as shown in the diagrams above. But when each metatile
expands into a cluster of hats, the clusters aren’t quite
exactly the same sizes as the metatiles they replace. So the
precise positions in the plane will vary in a nonlinear way, and
there’s no simple formula that matches a metatile’s
coordinates to the coordinates of its hats.

The paper offers an alternative version of this scheme in which
the smallest level of metatiles do match up
geometrically to the hats they turn into – but in that version,
the downside is that each size of metatiles is not precisely the
same shape as the next size. Instead, the metatiles themselves
gradually distort as they expand.

Either way, this makes it more difficult to perform the overall
procedure. With the fixed-shape metatiles as shown here, in the
final conversion to hats, you can’t just iterate over each
metatile and independently write out the coordinates of its
hats; with the alternative version, you have the same problem
during the expansion of metatiles into smaller ones. In each
case, you have to do some kind of a graph search (breadth-first
or depth-first, whichever you find easiest), choosing one
metatile to start with and processing the rest in an order that
guarantees that each metatile you process is adjacent to one
you’ve already placed. That way you can figure out the location
of each output hat or metatile, by knowing of an existing
adjacent thing it has to fit together with.

Here’s an example, which also illustrates the distortion. This
shows the same set of metatiles as in step 2 of the iteration
above, and its expansion into hats. The top left corners of the
two images are aligned – but you can see that at the other end,
the corresponding pieces of tiling are quite offset from each
other.

Metatiles
Hats

Choosing a random patch from a fixed expansion

In the previous section, I’ve shown how all three of our
aperiodic tiling types (P2, P3, hats) can be generated by
repeated subdivision of a starting tile, followed by a final
postprocessing step. (In the two Penrose cases, gluing
half-tiles back together into whole ones; in the hat case,
substituting the four metatiles for their hat clusters.)

This is all very well if you want to generate a fixed
piece of tiling. But what if you want a different piece every
time?

One obvious answer is: make a fixed patch of tiles much
larger than you need, and select a small region from it at
random.

This works basically OK in principle. But if you implement it
naïvely, it’s very slow, because you spend a lot of effort on
generating all the fine detail in huge regions of the tiling
that you’re just going to throw away. And the more uniform a
probability distribution you want for your output region, the
larger the fixed region you have to generate, and the more
pointless effort you waste.

We can save a lot of that effort by not being quite so naïve
about the procedure. Suppose that, instead of generating a patch
of tiling and then deciding which part of it to return,
we instead decide in advance where our output region is
going to be within the starting tile. Then, in each
iteration, we can identify a tile that’s completely outside
the output region and won’t contribute any tiles to it, and
discard it early, before we waste any more time expanding it
into lots of subtiles that will all be thrown away.

For example, here’s an illustration of how much work this might
save when selecting a small rectangle from the same P3 example
pictured above:

0
1
2
3
4
5
6
7
8

Without pruning

At every stage, we identify triangles – large or small – that
are out of our region, and stop bothering to expand them any
further. This bounds the wasted effort at a much more acceptable
level. Compare the number of triangles in the final pruned
output with the number in the unpruned version!

This technique is easy to apply for Penrose tilings, because
the expansion process is geometrically exact. At every stage of
the process, the outline of each triangle is exactly the
boundary of the region that will contain all its eventual output
tiles. So it’s easy to identify whether a triangle
intersects the target region.

But in the hats tiling, where the expansion process is
geometrically distorting, it’s not so easy. In order to work out
whether a given metatile is going to contribute to your target
region, first you have to work out what coordinates in the
current layer of metatile expansion correspond to the
target region in the final hat tiling. This is computationally
tricky: even if you could find an exact formula that translates
each corner of your target region into its coordinates
at the current iteration, then you’re not done, because it’s
also not guaranteed that the image of a straight line between
two of those corners will still be straight after the distortion
is applied. In other words, if you’re trying to construct a
patch of hats to fill a specific rectangle, you might find the
region of metatiling you need to keep isn’t even
rectangular!

It might be possible to get this technique to work regardless,
by computing a conservative approximation to the target region.
That way maybe you do a little more work than you need,
computing a few metatiles that never contribute hats to the
output but your approximation couldn’t quite prove it; but the
effort saving would still be large compared to the naïve
approach. But this is likely to be tricky, and has plenty of
room for subtle bugs. Err in one direction, and you do more work
than you need; err in the other direction and you generate wrong
output.

This entire approach has another downside. As I mentioned
above, if you want your random patch of tiling to be
selected uniformly from the limiting probability
distribution of finite pieces of tiling over the whole plane,
then you can’t achieve it this way, in an exact manner. You can
only approximate, by running a large enough number of iterations
that the huge patch you’re selecting a rectangle from is
close enough to that limiting distribution. And the
closer you want the distribution, the more you have to expand,
and the more work you have to do.

More subtly, this means you have to figure out in
advance how many iterations you want to run, based on the
size of your target region and the error threshold you’re
prepared to tolerate in probability distribution. That
involves a fiddly formula which is easy to get wrong,
introducing subtle bugs you’ll probably never spot.

In the next section, I’ll describe a completely different
approach that avoids all these problems, and works for
hats as easily as it does for Penrose tiles.

Combinatorial coordinates

The alternative approach begins with the following observation.

Suppose that, in any of these recursive substitution systems,
whenever we substitute one larger tile for a number of smaller
ones, we assign a different label to each smaller tile.

For example, in the P2 tiling, an acute triangle expands into
three smaller triangles; so we could number them 0, 1, 2 in some
arbitrary but consistent way. And an obtuse triangle expands
into just two, which we could call 0 and 1. Or, perhaps more
elegantly, we can observe that no two triangles expanded from
the same larger triangle are the same type (counting reflections
as different), so we could simply label each triangle by one of
the four types (two types of acute, two types of obtuse).

For the metatiles in the hat system, each metatile has between
7 and 13 children, and some of them are the same type as each
other. So the problem is larger – but not more difficult. We can
still just assign each child a numeric label in an arbitrary
way, so that the hexagonal metatile has children labelled from 0
to 12 inclusive, the triangle from 0 to 6, etc.

Then, after you’ve run the expansion process some particular
number of times, each tile of the output can be uniquely
identified by the sequence of labels it generated as it was
expanded. You could identify some particular tile by saying
something like: “From the starting tile, take the child
with this label; from there, take the child
with that label, then the child with the other
label, …” and after you’ve listed the label at every step of
the expansion, your listener knows exactly how to find the same
tile you were thinking of.

Now, supposing I tell you the coordinates of one particular
output tile, specified in this way. Can you figure out the
coordinates of one of its neighbours, sharing an edge
with it?

Assuming the tile you want is not on the very edge of the giant
starting tile, yes, you can. And you can do it entirely by
manipulating the string of labels, without ever having to think
about geometry.

Penrose tilings: indexing triangle edges

We’ll tackle the Penrose tilings first. In this system, the
smallest unit we ever deal with is a single half-tile triangle,
so we’ll assign coordinates to each of those.

I’ll label the tile types with letters. For the P2 tiling, the
two handednesses of acute triangle are called A and B; for the P3
they’re C and D. The obtuse triangles in P2 will be U and V; the
ones for P3 will be X and Y.

Here’s an illustration that shows the coordinates of each
output triangle for the first few layers of expansion, in both
tiling types. The coordinates are written with the label of the
smallest triangle first, followed by its parents in increasing
order of size. (This matches the natural way to store them in an
array, because that way, the fixed array index 0 corresponds to
the real output triangles.) So whenever a triangle is divided
into two or three smaller triangles, the coordinates for each
smaller triangle are derived from the original triangle by
appending a new type label to the front.

P2:
0
1
2
3

P3:
0
1
2
3

We’ll also need to index the edges of each triangle in
a unique way, so as to describe unambiguously which neighbour of
a triangle we’re going to ask for. For example, we could say
that we always assign index 0 to the base of the isosceles
triangle, and the two equal sides are assigned indexes 1 and 2,
in anticlockwise order from the base. So when you want to know
about a particular neighbour of your starting tile, you’ll
identify which neighbour by giving the index of an edge: “What
triangle is on the other side of edge {0/1/2} of my current
one?”

Then, for each substitution rule, we can make a map of which
smaller triangles border on which others, and which of their
edges correspond:

From these maps, with a little recursion, you can work out
everything you need.

For example, suppose you’re in a type B triangle, and you want
to see what’s on the other side of its edge #0. Look at the type
of its parent triangle. If the parent is type A or type U, then
the maps above for those triangle types show that edge #0 of the
type-B sub-triangle is also edge #1 of a type-U triangle. So you
can adjust the lowest-order label in your coordinate list to
specify a different sub-triangle, and return success.

More symbolically: to move across edge #0 of any triangle
whose coordinate string begins with BA or BU, just replace the
initial B with a U.

On the other hand, suppose your type-B triangle is expanded
from another type-B triangle. In that case, the map for
the parent triangle says it doesn’t know what’s on the far side
of your edge #0 – that’s off the edge of the map.

But that’s all right: recurse further up to find a map on a
larger scale! We know that going out of edge #0 of the smaller
type-B triangle means going out of edge #1 of the larger type-B
triangle. So ask the same type of question one layer further up
(perhaps recursing again, if necessary).

When the recursive call returns, it will give you an entire
updated set of coordinates identifying the second
smallest triangle you’re going to end up in, and which edge of
it is on the other side of this one’s edge #1. So, finally, you
can append the lowest-order coordinate to that, by figuring out
which of the smaller triangles of the new large
triangle is the one you need to end up at.

In our example case, there’s one final step. Going out of edge
#0 of our original B meant going out of edge #1 of the larger B.
But that edge is divided into two segments, each belonging to a
different sub-triangle. So we need to remember which
segment of the larger edge we were crossing: were we to the left
or the right of the division?

We’ll always expect to find that the incoming edge of the new
large triangle is subdivided into segments in the same way, and
the map for that triangle will let us find which sub-triangle
corresponds to each segment of the outer edge. So you can still
figure out which sub-triangle you end up in.

For example, suppose we had coordinates ending BB, and our
recursive call (“we’re going out of edge #1 of a B, what happens
next?”) rewrote the second B to an A (perhaps modifying further
labels above that) and told us we were coming in along edge #2
of the A. Then we consult the A map, and we see that its edge #2
is indeed divided in two, and coming in via the right-hand one
of those segments leads us in the A child. So we end up with AA
as our new lowest-order coordinates (plus whatever rewrites
were made at higher orders by the recursion).

(You can check these examples in the expansion shown above. In
the 3-iteration P2 diagram, there are triangles labelled BABU
and BUUU, and as described above, the edge #0 of each one – that
is, its short base edge – connects to UABU or UUUU respectively,
with only the lowest-order label differing. But the triangle
labelled BBBU is a more difficult case, requiring a rewrite at
the second layer: its edge #0 connects to a triangle labelled
AABU.)

Hats: indexing kite edges, and dealing with non-uniqueness

For the hats tiling, indexing the edges of the tiles isn’t
quite so convenient. The hat polygon itself has 13 edges, and
the maps for which ones border on each neighbouring hat are
pretty complicated. Not only that, but because the metatile
expansions overlap each other, every tile potentially has
multiple different, equally legal, lists of coordinates. And the
metatiles don’t necessarily meet edge-to-edge, in the sense that
one edge of this metatile might join up with the edges
of two other metatiles.

To solve the first of those problems, I found the easiest thing
is to stop considering a whole hat at a time, and use a smaller
and more regular unit. The hat tiling can be made to align
neatly with a periodic tiling of the plane with kites, shaped so
that six of them joined at the pointy ends make a regular
hexagon and three joined at the blunt ends make an equilateral
triangle:

If you align the hats to this tiling, then every hat occupies
exactly eight whole kites, as shown above. And the kites are
simply shaped, with just four edges. So I found it’s easier to
consider each individual kite to have a set of
combinatorial coordinates, and to devise an algorithm that will
tell you the coordinates of a neighbouring kite, given the
coordinates of the starting kite and which of the four edges you
want to head out of.

The coordinates of an individual kite look something like this:

• “I am the kth kite in a hat …
• which is the hth hat expanded from a first-order
  metatile of type m …
• which is the cth child of a second-order
  metatile of type m₂ …
• which is the c₂th child of a third-order
  metatile of type m₃ …”

and so on. At the highest-order end, this terminates with you
knowing the outermost metatile type, but not anything
about its parent, or which child of that parent it might be.

To navigate this coordinate system, we’ll make a set of giant
lookup tables. But first, we’ll have to assign numeric indexes
to all three of the expansion processes involved: metatiles to
smaller metatiles, then to hats, then to kites. It doesn’t
matter how we do this, as long as we do it consistently.

The first type of lookup table we’ll need, I call
a kitemap. For each type of second-order metatile (i.e.
each possible value of m₂), you expand it into
its set of first-order metatiles, and then into hats, and then
into kites, and you keep the coordinate labels you generated at
each stage. To illustrate this, here’s the kitemap for the
triangular metatile. (It make a manageable example, because it’s
the smallest. The others are similar but larger.)

Starting metatile
Expand metatiles
Expand to hats
Expand to kites

From this visual map, you can read off the coordinates of each
kite adjacent to a starting kite. For example, consider the kite
labelled 7.3.0, at the top right of the central hat in this
diagram: you can see that the four kites bordering it are 0.0.0,
6,3,0, 3.1.1 and 4.1.1, and you can identify which edge of the
kite takes you to each of those neighbours. So if the
coordinates of your current kite started with
(k, h, c) = (7, 3, 0), then you would know
how to generate the coordinates of each neighbouring kite,
simply by rewriting those three low-order labels.

The version of the kitemap used by the algorithm is not a
visual map like this: it’s a lookup table, or rather one for
each of the four metatile types. Each table is indexed by the
triple (k, h, c) and a particular kite
edge, and it tells you the new values of
(k, h, c) corresponding to the kite on the
far side of that edge.

So, to compute the coordinates of a neighbouring kite, you
start by choosing the kitemap lookup table corresponding to the
metatile type m₂, and look up
(k, h, c, kite edge) in it. If it has an
answer, you’re done – just replace the three low-order indices
in your input coordinates with the ones you got out of the
kitemap lookup table, and you’ve got a set of coordinates for
the new kite.

(Of course, when you rewrite c in the coordinates, you
must also rewrite m to match the type of the new
first-order metatile. To do this, you just need a much simpler
lookup table of which type of metatile each child is.)

But sometimes this doesn’t work. For example, suppose you’re in
the kite labelled 0.0.4 in the above map, right at the bottom.
Then traversing one of the four kite edges will take you inwards
to 1.0.4, but for the other three, this map has no answer. In
cases like this, the entry in the kitemap lookup table for “what
happens if I head in this direction
from here?” will be a special value meaning “don’t
know, that’s off the edge of the map.”

By analogy with the Penrose case, an obvious thing to try here
might be to specify (in some way) which edge of the
kitemap you’re stepping off, and then recurse to a higher layer
to find out where you’ve come back in to some other kitemap. If
we had made smaller kitemaps based only on a
single first-order metatile (i.e. just containing 1, 2
or 4 hats), then we would have no choice but to do exactly this,
because all first-order metatiles’ hat expansions are disjoint.
But hat expansions of metatiles have extremely twiddly and
complicated edges, and that sounded like a very tricky thing to
get right.

And we don’t have to! There’s a nicer way, based on the fact
that the metatile expansions overlap each other. If we’re about
to head off the outermost border of the kitemap, that must mean
that we’re near the edge of the expansion of our second-order
metatile m₂. So the first-order
metatile we’re in must be one of the outermost metatiles in that
expansion – which means it overlaps with the expansion of
another second-order metatile, or maybe even two of them. And
among the second-order metatiles whose expansions contain this
particular first-order one, there must be at least one in which
we’re about to step inwards towards the centre, not
outwards towards the edge.

In other words, the fact that there are multiple valid
coordinates for the same kite is not a problem after all. It’s
an opportunity! If we can rewrite two layers
of metatile child index in our coordinate system, then
we can find an equivalent representation of our current
kite, which places it in a kitemap that we’re not about to step
off the edge of. And now we’ve completely avoided having to deal
with those long crinkly complicated edges of all the maps.

To implement this step, we make a second type of map, which I
call a metamap. This time, take
each third-order metatile type m₃,
and expand it twice, into second- and then into first-order
metatiles. As you do that, keep track of all the different ways
that each first-order metatile is generated (there will be more
than one if it’s in the overlap between the expansions of two or
three second-order metatiles). So we end up with a map in which
some metatiles have multiple coordinates:

Starting metatile
Expand once
Expand twice

Finally, we translate that map into a lookup table indexed by
any of the coordinate pairs in this diagram, giving all the
other coordinate pairs equivalent to it. This lookup table
allows you to rewrite the part of the coordinate system that
says

“… metatile type m, which is
child c of type m₂, which is
child c₂ of m₃ …”

by selecting the metamap for tile type m₃,
and looking up the pair (c, c₂) in it.
The metamap might return you one or two alternative
(c, c₂) pairs, and for each one, you
can substitute those in your coordinates (which doesn’t change
which kite you’re referring to). This will also give you a new
value of m₂, so you’re representing your
existing location relative to a different kitemap which overlaps
the previous one. Now you can try looking up your new
(k, h, c) triple in it (using the rewritten
value of c), and this time, perhaps you’re not on the
edge of the map any more, and can find the coordinates of the
neighbouring kite successfully.

If even that doesn’t work, it’s because you’re on the
outer edge of the metamap, meaning you’re near the edge of the
expansion of the third-order metatile m₃.
And now every layer above this looks the same – so you
can try applying the same kind of metamap rewrite one layer up,
by looking up (c₂, c₃) in
the metamap for m₄ and seeing if that permits
a more useful rewrite one layer down. If even that doesn’t help,
try further up still, and so on.

Making it up as you go along

In both of the previous sections, I’ve described a technique
for computing the coordinates of each smallest-size element of
the tiling (a Penrose half-tile, or a single kite), using a
recursive algorithm which looks at higher- and higher-order
labels in the coordinate system as necessary.

Of course, there’s one obvious way this can fail. You can only
store a finite number of coordinates. What if the recursion tries
to go all the way off the top, and look at a higher-order
coordinate than you have at all?

One obvious option is: report failure. If for some
reason you had decided in advance that you were only
interested in some specific fixed patch of tiling, and that
nobody should ever go beyond the edges of that at all, then a
failure of this type should be interpreted as “Clonk! You have
run into the wall, and can’t go that way.”

(Perhaps if you had decided to set a text adventure game in the
7-level expansion of a type-A Penrose half-tile, or something
along those lines, then this might be what you’d actually want!
Each individual triangle would correspond to a room, and not all
rooms would have exits in all three directions.)

But in the main situation I’m discussing here, we don’t ever
want to report failure. Our aim is to generate a randomly
selected patch of tiling to cover a particular target rectangle.
Just because that rectangle turns out to protrude over the edge
of the highest-order tile we currently know about doesn’t mean
there’s no possible way to extend the tiling in that
direction; it just means we haven’t yet made a decision
about which way to extend it.

If you imagine an infinite tiling of the entire plane
in any of our tiling systems, then in principle, any given
element has a set of coordinates that go on forever to
higher- and higher-order, larger and larger, supertiles. If we
had started by inventing an infinite sequence of
coordinates for our starting element, then there would be no
problem with running out of coordinates during the recursion –
we could recurse as high as necessary.

Of course, you can’t invent infinitely many things up front in
a computer. But what you can do is to invent them
lazily, by inventing the first few, and being prepared
to generate a new one as and when it’s requested.

So this is what we do. In our algorithm, we invent some random
coordinates of our original starting element, up to a certain
point. And then, if the recursion ever tries to go beyond the
level we know about, we simply invent another layer, on demand,
append it to the coordinates of the starting element,
and pretend it was there all along, and that our
starting element has always had that extra coordinate.
So if any other coordinate list lying around elsewhere in the
software is shorter than the current coordinates of the starting
position, then we can extend it by appending the extra elements
that have been appended to the starting position since it was
generated. (Because it has the semantics “When we last got here,
we didn’t know what the next few levels looked like, but now we
do know”.)

Of course, when you make up a parent coordinate at random, you
have to make it consistent with the rules about which tiles can
be parents of which other tiles. In the Penrose system, for
example, you can see from the maps in a previous section that
the type-B triangle can be a child of A or U or another B, but
it can’t be a child of a V. So if your current topmost tile type
is B, then you must select the next one up by randomly choosing
from the set {A, U, B}, and not from all four tile types.
Similarly in the hats system, not all metatiles can appear as
children of other metatiles: the triangular
metatile only appears in the expansion of the hexagonal
one, so if your current topmost metatile is a triangle, then you
only have one choice for what the next larger one will be. And
even if your current metatile occurs as a child of all
the metatile types, you also have to decide what its child index
will be relative to its parent, and that index must be chosen in
a way that matches the tile types.

Better still, once you’ve finished generating your output patch
of tiling, you can retrieve the full set of coordinates you
ended up having to invent for the starting position – and those
can act as a short and convenient identifier that allows another
person running the same algorithm to generate the identical
piece of tiling, because they should never find
themselves needing a coordinate you hadn’t already generated.
That would only happen if they tried to extend your tiling into
a larger region than you had tried yourself.

Another thing you can do here is to choose the next parent
using a deliberately biased probability distribution, to match
the limiting distribution over the whole plane.

In each of these tiling systems, there’s a system of linear
equations that tells you how many tiles of each type you have
after an expansion, if you know how many tiles of each type you
started with. By eigenvalue/eigenvector analysis, you can
process these equations to determine the proportions of the
various tile types occurring overall – or, more precisely, the
limit of the proportions in a very large finite patch, as the
patch size tends to ∞.

This overall distribution of tile types, in turn, can be
processed into a set of conditional probabilities, each
answering a question of the form “Of all the tiles of
type t in the plane, what proportion are child #n
of tile type u?” So if you work out all those
probabilities, then you can select the next coordinate pair
(n, u) in a way that matches that limiting
distribution.

(Also, when you randomly select your initial tile, you
should do it according to the original, unconditional, version
of the limiting distribution.)

The effect of this refinement is that our piece of tiling will
be generated as if the coordinates of the starting element had
been chosen uniformly at random from the whole plane!
That’s normally a meaningless concept, but in this case it (just
about) makes sense, because the number of output patches of
tiling you could possibly return is finite, and the distribution
of those converges to a limit if you choose uniformly
from larger and larger regions of the plane.

Compare this to the recursive-expansion technique for
generating tilings. In that technique, your output distribution
is – by construction – the one you’d get by choosing from
an actual fixed patch of tiling of a finite (if large)
size, because that’s exactly what you are doing. If you
start from a larger and larger patch, your output
distribution approaches the idealised limit, but in
order to get it closer and closer, you have to do more and more
work. Here, we can directly generate a tiling
from precisely the limiting distribution, and we didn’t
have to do any extra work at all to get there!

Generating the output tiling

So, now we have a method for moving around a tiling and
determining a string of coordinates describing each
smallest-sized element we come to. What do we do with
that?

Our ultimate aim is to actually generate a set of tiles that
cover the specified target area. For this, we don’t really need
the full set of higher- and higher-order coordinates at all. We
only need to know the lowest-order tile type labels: enough to
know which kind of thing appears next in the tiling.
Those low-order labels are the output of the algorithm.
The higher-order parts of the coordinate system are just an
internal detail of how the algorithm keeps track of where it had
got to. (Also, as I mention in the previous section, they can
also act as an identifier to make a particular run
repeatable.)

So, to use this technique to actually generate a Penrose
tiling, a simple approach is to use a graph search (say,
breadth-first). Pick an initial starting triangle, and place it
somewhere in your target area (the centre is an obvious choice,
but anywhere will do). Then, repeatedly, find some triangle edge
that you don’t know what’s on the other side of yet, and compute
the combinatorial coordinates of the triangle on the other side,
by starting from the coordinates of the triangle you already
have. The lowest-order label in the new triangle’s coordinates
will tell you which shape of triangle it is (acute or obtuse),
and which of its edges corresponds to the edge you just
traversed. That’s enough information to calculate the
coordinates of the new triangle’s third vertex, starting from
the two endpoints of that edge. Now add the two new edges of
that triangle to your queue of edges to try later. Then go back
and pick another edge to explore, and so on.

At every stage, if you generate a triangle that’s completely
outside the target area, there’s no need to store it or to queue
up its other edges. If you discard out-of-bounds triangles as
you get to them, then this search process will eventually
terminate, because the queue of edges that still need exploring
has run out, and you’ve covered the whole target region in
half-tile triangles.

Every time you generate the coordinates of a triangle, you find
out which of the four triangle types it is: not just
whether it’s acute or obtuse, but also which handedness
it is. So you know which half of an output tile each triangle
is. Therefore, as you go along, you can also generate the full
output Penrose tile corresponding to every triangle. This will
generate most of them twice (any tile with both halves inside
the target area), but as long as you notice that and
de-duplicate them, that’s fine.

For generating the hat tiling, you could do the same
kind of thing – but there’s an even easier approach that doesn’t
need a graph search at all, using the fact that the whole tiling
lives on a fixed grid of kites.

For the hat tiling, the easiest approach is to specify your
target region as a connected set of kites, and simply iterate
over that set in some fixed way. For example, to fill a
rectangular region, you might process the kites in the region in
raster order, iterating along each row from left to right and
then going on to the next row. Or you might find it easier to
process alternate rows in opposite orders. As long as you
process the kites in an order that means every kite (after the
first one) is adjacent to a kite you’ve already found the
coordinates of, anything will work. At the end of the iteration,
you know the combinatorial coordinates of every kite in your
region.

Once you’ve done that, you can generate the output in whichever
way is easiest. For every kite in the region, you know which
kite it is out of the 8 making up its containing hat, and that
(together with knowing whether the hat is a reflected one) is
enough to generate its complete outline. So you could do that
for every kite in the region, and discard duplicates.

But an even easier way than that is to only generate an output
hat when the iteration finds a kite which was #0 in its hat, and
even then, discard it if that whole hat doesn’t fit in the
output area. This guarantees to generate every valid output
hat exactly once, so no de-duplication is even
needed.

(If you need to generate hats that cover your output
area, instead of the subset of hats that fit entirely inside it,
then that algorithm can still be used – just expand the output
region by the maximum width of a hat on each side, so that any
hat that intersects the true region must fall entirely within
the expanded one.)

As long as you’re only looking at the three coordinates
(k, h, m) of each kite (which kite it is,
in which hat, of which type of first-order metatile),
it doesn’t matter that multiple different coordinates can refer
to the same location, because they’ll all agree on those three
points. The lowest-order coordinate index that can differ
is c, the question of which child of which second-order
metatile the first one is. So the information you really need
(kite index, and whether the hat is reflected) is easy to read
off reliably.

Advantages

I’ve mentioned in passing a few ways that this combinatorial
coordinate technique is superior to recursive expansion. But
let’s bring them all together:

Uniform distribution. By generating each
coordinate with the correct probability distribution, you can
arrange that your output patch of tiling is chosen
from precisely the overall limiting distribution, and
not just some approximation to it based on a large finite
region.

No need to engage with geometry of higher-order
tiles. The question of what points in the
plane might be the vertices of any higher-order tile just
never comes up, in this system. The only thing we ever track
about the higher-order tile layers is
their combinatorial nature: what tiles exist at all,
which ones are children of which, which ones are adjacent to
which, which ones overlap which. This isn’t such a huge
advantage in the Penrose tilings, where the higher-order
geometry was still fairly easy. But in the hats system, the
distortion of the metatiles is made completely irrelevant by
treating them purely combinatorially, so we never have to worry
about it at all.

No risk of overflow. In the recursive
expansion technique, we need to generate actual coordinates for
all the tiles at every layer. If we’re choosing a small output
region from a really big patch, this might mean we have to worry
about floating-point precision loss in the largest scales – or
alternatively integer overflow, if we use exact integer-based
coordinates (see the appendix below). But in this system, the
only output points we ever generate are the ones in the actual
output region, so we only need to ensure that our number
representation is accurate enough for that. The much
larger tiles surrounding the output region are still described
in this algorithm, but instead of being described by very large
single numbers (representing vertices in the plane), they’re
described by a long string of small numbers (the labels in the
combinatorial coordinate list). It’s as if we had transferred
the working state from machine integers or floats into a bignum
representation – just a slightly weird two-dimensional
mixed-base one.

No manual tuning needed. In the recursive
expansion technique, you have to think about the size of output
region you want, and choose up front how many recursion levels
you’re going to perform (also based on what error you’re
prepared to tolerate from the ideal probability distribution).
But in this system, you find out after the algorithm
finishes how many levels of coordinates it turned out to need to
invent. You don’t have to make a decision up front at all, which
means you don’t need to debug the code that does it.

(In particular, you might not end up with the same
number of coordinate labels, in two runs on the same target area
with different random numbers. It will depend on whether the
starting element’s randomly chosen coordinates happened to put
it close to a high-order tile boundary.)

Close to linear time. The number of coordinate
labels we end up dealing with is variable. But on average it
will be proportional to the log of the number of output tiles.
For each output tile we generate, we perform a bounded number of
operations that step from the coordinates of one element to an
adjacent one. Each of those steps will
take O(log n) time, or
maybe O(log² n) for the hats system
which might need to recurse back and forth rewriting metatile
labels. But that’s only an upper bound: most coordinate
rewrites will be able to succeed by adjusting only the
lowest-order labels, and the occasions of recursing higher up
will become progressively rarer as you ascend the layers.

I don’t have a formal proof that the average overall running
time is any better
than O(n log² n). But it seems
intuitively clear that it will in practice be much closer to
linear time than that expression makes it look!

Only needs to use logarithmic memory. If you’re
generating a hats tiling to cover a rectangular region by
iterating over the kites in that region, you can arrange to keep
a very small number of coordinate lists live at one time – just
two or three is enough to guarantee that every kite you process
is adjacent to one of a small number of past kites that you’re
still keeping the coordinates for. There’s no need to store the
full coordinates of every kite for the duration of the
algorithm: you can emit each hat as you come to it in a
streaming manner, and forget nearly everything about the ones
you’ve already generated.

So this system is suitable for generating absolutely
enormous patches of hat tiling, with only a modest cost
in memory usage. The running time is close to linear in the
number of output hats, and the average memory usage is
just O(log n). Probably whatever is printing out
the hats after you generate them will have higher cost than
that!

I think it should be possible to work this way in the
Penrose case too, even without a convenient underlying grid to
iterate over. But it would be fiddly. I think you could do it by
iterating along a sequence of horizontal lines of the output
region, at a spacing small enough to guarantee any triangle must
intersect at least one line. Then, instead of adding
both of each new triangle’s outgoing edges to a queue,
instead choose just one of them to traverse next, by picking the
one that intersects the horizontal line you’re currently
following. Meanwhile, a similar iteration in the perpendicular
direction is finding all the triangles on the left edge of the
region, to begin each horizontal iteration at. To avoid emitting
any triangle twice, check to find how far above the lowest point
of the triangle the current line hits it; if that’s greater than
the line spacing then you know a previous sweep must have
already caught it.

I’m fairly sure this would work, but I haven’t tried generating
Penrose tilings in this low-memory mode. I’ve only tried the
more obvious breadth-first search. Also, this would be slower by
a constant factor in payment for the memory saving – you’d visit
most triangles more than once.

Conclusion

I’ve presented two algorithms for generating aperiodic tilings
at random. The first algorithm, picking a random region from a
large fixed area, is efficient enough for Penrose tilings, but
works rather badly for hats because of the distortion problem.
The second, combinatorial coordinates, is efficient for both
types of tiling, and especially so for hats, where the
other algorithm performed worse. Also, it has many other
assorted advantages, listed in the previous section.

The combinatorial coordinates system described here is in use
in the live version of Loopy for generating hat tilings, and is
working very nicely.

The Penrose tilings in Loopy are generated by the other
algorithm, picking a random region from a fixed patch. At the
time Penrose tilings were implemented in Loopy, neither I nor
the implementor had thought of the combinatorial coordinates
system.

Now that I have thought of it, it’s tempting to rewrite Loopy’s
Penrose tiling code using the shiny new algorithm. But on
balance I’m not sure that I want to, because its textual game
descriptions are based on specifying the coordinates of the
random sub-region. So if I threw out the existing code, then
existing game descriptions and save files would stop
working.

(I could keep the old code, solely for parsing old
game descriptions, and switch to combinatorial coordinates for
new ones. But that seems wasteful in another sense – I’d end up
with two pieces of complicated code where only one is
really needed!)

However, I have implemented a combinatorial-coordinate Penrose
generator as prototype code, and it’s pretty convenient. If I
were starting from scratch in Loopy with no backwards
compatibility considerations, I’d certainly choose that
algorithm now!

Appendix: exact coordinate systems

In both the Penrose and hat tilings, it’s useful to have a
representation for storing the coordinates of vertices of the
tiling which allows you to compute without floating-point
rounding errors. That way, you can easily check whether a vertex
of the Penrose tiling is the same one you’ve seen before, or
keep track of the coordinates of the current kite in the
underlying grid of the hats tiling, without having to deal with
the awkwardness that when you get back to something that
logically should be the same tile or vertex, its
coordinates have changed a tiny bit.

Complex roots of unity

For the Penrose tiling, one neat approach is as follows.
Consider the coordinates to be complex numbers, rather
than just 2-element vectors. Then let t be the complex
number cos(π/5) + i sin(π/5), which has
modulus 1 and argument π/5. Then, multiplying any other
number by t has the effect of rotating it by exactly a
tenth of a turn around the origin.

Therefore, we know that t⁵ = −1,
because rotating by a tenth of a turn five times must rotate by
half a turn, exactly negating the number you started with.
So t is a zero of the
polynomial z⁵ + 1. This polynomial
factorises as
(z + 1)(z⁴ − z³ + z² − z + 1),
and t is not a zero of the first factor (or else it would
just be −1). So it must be a zero of the second factor, meaning
that t⁴ = t³ − t² + t − 1.

So, suppose you have two complex numbers, each expressed as a
linear combination of the first four powers of t,
say a + bt + ct² + dt³
and w + xt + yt² + zt³.
Then you can multiply those polynomials together, multiply out
the product, and reduce every term
containing t⁴ or higher by applying the
identity t⁴ = t³ − t² + t − 1.
This gives you the product of the two numbers, also in
the form of a linear combination of the first four powers
of t. Better still, if the the two values you’re
multiplying have all four coefficients integers, then so does
the product.

Also, t has further useful properties. It turns out
that t + 1/t = φ, the golden ratio. This is
also the ratio of the two side lengths of the triangles making
up a Penrose tiling. And 1/t is
just t⁹ = −t⁴ = 1 − t + t² − t³.
So φ itself has a representation in this system, with all
four coordinates
integers: φ = 1 + t² − t³.
So does its inverse, because
1/φ = φ − 1 = t² − t³.
So you can scale a vector length up or down by a factor
of φ, by simply multiplying by an appropriate tuple of
four integers, and still leave all four of its
coordinates as integers in this system.

This gives you everything you need to calculate the vertices of
Penrose tilings without rounding error, using either of
the methods described in this article (recursive subdivision of
a large starting tile, or breadth-first search via combinatorial
coordinates). The coordinates of a single triangle can all be
computed from each other by a combination of rotating
by t and scaling up or down by φ; once you have
one triangle, you can similarly compute the coordinates of the
triangle next to it; and in the subdivision process the
different sizes of triangle are scaled by φ each time.
And every coordinate you need will be described by a tuple
(a, b, c, d) of four integers, so
they’re easy to check for equality!

The hat tiling is simpler, because it lives entirely on the
underlying kite tiling. There are lots of ways to represent the
coordinates of vertices of that tiling, and they need not be
mathematically clever; any old ad hoc approach will be
good enough. However, a similar trick to the above is a
particularly nice approach:
let s = cos(π/3) + i sin(π/3)
represent a sixth of a turn around the origin, in which
case s² = s − 1 for similar reasons to
above. Then you can represent all your vertices as integer
combinations of 1 and s, and compute products in the same
way as before, so you can rotate a sixth of a turn by
multiplying by s or by
1/s = s⁵ = −s² = −s + 1.

Separate x and y coordinates using √5

One thing that isn’t convenient, in the above Penrose
representation using t, is checking whether a vertex is
inside or outside your intended target rectangle. The four
coordinates in the tuple
(a, b, c, d) are pretty abstract,
and have no obvious relation to x and y
coordinates in the plane.

Eventually, of course, you’ll have to convert them back to
actual complex numbers, to plot your tiling on the screen or the
page. This is done by these formulae, which give you separate
real and imaginary parts (i.e. x and y
coordinates) for all the numbers involved:

t = cos(π/5) + i sin(π/5)


t² = cos(2π/5) + i sin(2π/5)


t³ = cos(3π/5) + i sin(3π/5)

So you could just use those formulae to convert all your
coordinates into floating point as you generate the tiling, and
use the floating-point coordinates to decide whether any part of
a tile is inside your output rectangle.

But if you’d rather do even that check without
rounding error, there’s a more exact method available. It
happens that cos(π/5) = φ/2 = (1 + √5)/4, and
cos(2π/5) = (1 − √5)/4. The sines appearing in the
imaginary parts don’t have nearly such nice algebraic
representations (in fact they’re roots of a horrible quartic),
but their ratio does: the ratio
sin(2π/5)/sin(π/5) turns out to be just φ
again.

(The t³ term can be reduced to the existing
numbers by observing that sin(3π/5) = sin(2π/5),
and cos(3π/5) = −cos(2π/5).)

So if we define X = 1/4 and Y = sin(π/5)/2
to be our basic units of distance in the x and y
directions respectively, then the formulae above reduce to

t = (1 + √5)X + 2iY


t² = (1 − √5)X + (1 + √5)iY


t³ = (−1 + √5)X + (1 + √5)iY

in which each of the x and y coordinates is
expressed as an integer combination of 1 and √5 (times some
fixed scale unit which we mostly ignore).

In this representation, you can exactly compare two
numbers of that form to determine which is greater. Start by
subtracting the two, giving a single number of the
form a + b√5; then the problem reduces to testing
its sign. If a and b have the same sign as each
other, then the answer is obvious; if they have opposite sign,
then you just need to know which of a and b√5 has
greater magnitude, which is the same as asking which
of a² and 5b² is bigger. And
you can do that calculation in integers!

So if you scale your output rectangle just once to write its
width as a multiple of X and its height as a multiple
of Y, then you can do the entire process of tiling
generation, including bounds checking, in exact integer
arithmetic, and convert back to floating point only for the
final job of actually plotting the points.

(You could also choose to use these √5-based coordinates for
the whole job, instead of converting from the t
representation. It’s slightly more inconvenient to do rotation
and scaling by φ in this system, because some of the
operations you need will involve integer division, whereas in
the t representation, it’s all just multiplication and
addition. But the divisions will be by small constant integers,
and if you generate correct coordinates then they’ll never leave
a remainder. So you might decide that this is less hassle
overall than having two complicated coordinate systems
and a conversion in between. It’s up to you.)

Appendix: interleaving the Penrose tiling types

Here’s an interesting thing about the P2 and P3 tilings.

Both tiling types have substitution systems that involve the
same two types of isosceles triangle, each in two mirror-image
forms: one with an acute apex angle, and one obtuse.

In P2, the acute triangle is the larger one, and the
substitution step divides it into three smaller triangles – but
two of them are the same two triangles that the obtuse
triangle divides into. So you could imagine dividing the acute
triangle first into a smaller acute triangle and an obtuse
triangle of the current size, and then completing the procedure
by dividing all the obtuse triangles.

In P3, the situation is exactly reversed. The obtuse triangle
is larger, and its three substituted tiles include two that
could also have been obtained by first making an acute triangle
and then subdividing that.

It turns out that, by separating each tiling’s subdivision into
these two phases, they can be interleaved! Here’s a single set
of substitution rules that divide both acute and obtuse
triangles into just two pieces:

The rule with this substitution system is that, in each
iteration, you only use the rules for one type of triangle –
whichever one is currently larger. If you have larger obtuse
triangles than acute, then you only divide up the obtuse
triangles; this makes the obtuse triangles smaller and leaves
the acute ones the same size (including the extra ones you just
made). That means, in the next iteration, the acute triangles
are larger, so this time you only subdivide those ones. So in
each iteration, you use whichever rule you didn’t use the
previous time.

Here’s an illustration of the resulting interleaved
substitution system. In each iteration, the appropriate
boundaries between mirror-image triangles are drawn fainter, so
you can see that in alternate iterations they combine into kites
and darts, or into rhombs, depending on which tile type is
currently the larger one:

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

This isn’t of practical use while constructing the tilings, but
isn’t it pretty?

References

Some useful links if you’d like to know more:

• Wikipedia on Penrose tilings. (I don’t have a primary source about these.)
• Website about the hat tiling by its discoverers.
• arXiv location of the paper about the hat tiling.