# Extending OEIS: Counting Diamond Tilings

I promise, this will be my last challenge about diamong tilings (for a while, anyway). On the bright side, this challenge doesn't have anything to do with ASCII art, and is not a code golf either, so this is actually completely different.

So just as a reminder, every hexagon can be titled with three different diamonds:

An interesting question to ask is how many of these tilings exist for a given hexagon size. It seems these numbers have been studied fairly thoroughly and can be found in OEIS A008793.

However, the problem gets trickier if we ask how many tilings exist up to rotation and reflection. For instance, for side length N = 2, the following 20 tilings exist:

   ____     ____     ____     ____     ____     ____     ____     ____     ____     ____
/\_\_\   /\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
/\/\_\_\ /\/_/\_\ /\/_/_/\ /\/_/\_\ /\_\/\_\ /\_\/_/\ /\/_/_/\ /\_\/\_\ /\_\/_/\ /_/\/\_\
\/\/_/_/ \/\_\/_/ \/\_\_\/ \/_/\/_/ \/\_\/_/ \/\_\_\/ \/_/\_\/ \/_/\/_/ \/_/\_\/ \_\/\/_/
\/_/_/   \/_/_/   \/_/_/   \_\/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/   \_\/_/   \_\/_/
____     ____     ____     ____     ____     ____     ____     ____     ____     ____
/_/_/\   /\_\_\   /_/\_\   /_/_/\   /_/\_\   /_/\_\   /_/_/\   /_/_/\   /_/_/\   /_/_/\
/\_\_\/\ /\/_/_/\ /_/\/_/\ /\_\_\/\ /\_\/_/\ /_/\/_/\ /_/\_\/\ /\_\_\/\ /_/\_\/\ /_/_/\/\
\/\_\_\/ \/_/_/\/ \_\/\_\/ \/_/\_\/ \/_/_/\/ \_\/_/\/ \_\/\_\/ \/_/_/\/ \_\/_/\/ \_\_\/\/
\/_/_/   \_\_\/   \_\/_/   \_\/_/   \_\_\/   \_\_\/   \_\/_/   \_\_\/   \_\_\/   \_\_\/


But many of these are identical under rotation and reflection. If we take these symmetries into account, only 6 distinct tilings remain:

   ____     ____     ____     ____     ____     ____
/\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
/\/\_\_\ /\/_/\_\ /\/_/_/\ /\_\/_/\ /\_\/_/\ /_/\/\_\
\/\/_/_/ \/\_\/_/ \/\_\_\/ \/\_\_\/ \/_/\_\/ \_\/\/_/
\/_/_/   \/_/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/

2        2        6        6        1        3


where the numbers indicate the multiplicity of each tiling. Note that for larger hexagons there are also tilings with multiplicities 4 and 12.

It appears that the number of tilings up to symmetry has been studied less thoroughly. The OEIS entry A066931 only lists the five terms:

1, 1, 6, 113, 20174


where the first term is for side length N = 0 and the last term for side length N = 4.

I'm sure we can do better than that!

Your task is to compute the number of tilings for a given side length.

This is . Your score will be the highest side-length N for which your code produces the correct result within 30 minutes on my machine. In case of a tie, I will accept the submission which produces the result for that N fastest.

As usual, you must not hardcode results you already know to win the tie-breaker. The algorithm that solve N = 3 should be identical to the one that solves N = 5.

Your submission must not use more than 4GB of memory. I will give some leeway on this if you're operating close to that limit, but if you're consistently above that limit, or if you spike significantly beyond it, I will not count that N for your submission.

I will test all submissions on my Windows 8 machine, so make sure your language of choice is freely available on Windows. The only exception to this is Mathematica (because I happen to have a licence for it). Please include instructions for how to compile/run your code.

Of course, feel free to compute more terms in your own time (for science, and for others to check their numbers against), but your answer's score will be determined in those 30 minutes.

• Note that since N = 6 gives an output of more than 10^12, a non-constructive solution is almost certainly necessary to get that far. – Peter Taylor Jun 4 '15 at 12:44
• @PeterTaylor I was hoping that would allow more room for improvement. Maybe a couple of simple constructive answers first that can do N = 5 to gain more insight into the problem, and then potentially hybrid approaches that don't need to construct all tilings but can extrapolate the total number from a few constructed ones... and then maybe something analytic if we're really lucky. :) – Martin Ender Jun 4 '15 at 12:50
• At the risk of stating the obvious, it seems to me that each such tiling corresponds to a projection of an assemblage of unit cubes as viewed from a distant point of view, e.g, from (100,-100,100). I find that this lightens the burden of constructing tilings. – DavidC Jun 4 '15 at 13:47
• @DavidCarraher Indeed. More specifically, such an arrangement of unit cubes is a 3D Young diagram. (Maybe that helps someone.) – Martin Ender Jun 5 '15 at 12:34
• @DavidCarraher If you look hard enough at the big hexagon, you'll see there's 2 different ways to interpret it as a Young diagram. The obvious way (for me at least) is to see a flat area at the top and left with a 2x2x1 cuboid missing from the top left corner. But there's another way of seeing it: an empty zone in that area, with a 2x2x1 cuboid sitting in it. Tilting 60 degrees may help. It hurts my eyes but I think the two young diagrams fit together, possibly by reflection of one of them. OEIS A008793 is very careful with its wording: "number of plane partitions whose young diagrams..." – Level River St Jun 5 '15 at 21:03

## Algebra, graph theory, Möbius inversion, research, and Java

The symmetry group of the hexagon is the dihedral group of order 12, and is generated by a 60 degree rotation and a mirror flip across a diameter. It has 16 subgroups, but some of them are in non-trivial conjugacy groups (the ones which only have reflections have 3 choices of axis), so there are 10 fundamentally different symmetries which a tiling of the hexagon can have:

The number of diamond tilings of a subset of a triangular lattice can be calculated as a determinant, so my initial approach was to set up one determinant for each of the symmetries of the hexagon, to calculate the number of tilings which have at least those symmetries; and then use Möbius inversion in the incidence algebra of their poset (basically a generalisation of the inclusion-exclusion principle) to work out the number of tilings whose symmetry group is exactly each of the 10 cases. However, some of the symmetries have nasty edge conditions, so I was forced to sum over exponentially many determinants. Fortunately, the values obtained for n < 10 gave me enough data to be able to identify relevant sequences in OEIS and piece together a closed form (for some value of "closed" which allows finite products). There's a bit of discussion of the sequences, and references to proofs, in the formal writeup which I prepared to justify the OEIS sequence updates.

Once the double-counting is taken care of, it turns out that four of the ten values cancel out neatly, so we only have to calculate the remaining six and then do a weighted sum.

This code takes under 30 seconds for N=1000 on my machine.

import java.math.BigInteger;

public class OptimisedCounter {
private static int[] minp = new int[2];

public static void main(String[] args) {
if (args.length > 0) {
for (String arg : args) System.out.println(count(Integer.parseInt(arg)));
}
else {
for (int n = 0; n < 16; n++) {
System.out.format("%d\t%s\n", n, count(n));
}
}
}

private static BigInteger count(int n) {
if (n == 0) return BigInteger.ONE;

if (minp.length < 3*n) {
int[] wider = new int[3*n];
System.arraycopy(minp, 0, wider, 0, minp.length);
for (int x = minp.length; x < wider.length; x++) {
// Find the smallest prime which divides x
for (wider[x] = 2; x % wider[x] != 0; wider[x]++) { /* Do nothing */ }
}
minp = wider;
}

BigInteger E = countE(n), R2 = countR2(n), F = countF(n), R3 = countR3(n), R = countR(n), FR = countFR(n);
return sum.divide(BigInteger.valueOf(12));
}

private static BigInteger countE(int n) {
int[] w = new int[3*n];
for (int i = 0; i < n; i++) {
for (int j = i + 1; j <= i + n; j++) w[j]--;
for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
}
return powerProd(w);
}

private static BigInteger countR2(int n) {
int[] w = new int[3*n];
for (int i = 0; i < n; i++) {
w[3*i+2]++;
for (int j = 3*i + 1; j <= 2*i + n + 1; j++) w[j]--;
for (int j = 2*i + n + 1; j <= i + n + n; j++) w[j]++;
}
return powerProd(w);
}

private static BigInteger countF(int n) {
int[] w = new int[3*n];
for (int i = 0; i < n; i++) {
for (int j = 2*i + 1; j <= 2*i + n; j++) w[j]--;
for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
}
return powerProd(w);
}

private static BigInteger countR3(int n) {
if ((n & 1) == 1) return BigInteger.ZERO;
return countE(n / 2).pow(2);
}

private static BigInteger countR(int n) {
if ((n & 1) == 1) return BigInteger.ZERO;
int m = n / 2;
int[] w = new int[3*m-1];
for (int i = 0; i < m; i++) {
for (int j = 1; j <= 3*i+1; j++) w[j] += 2;
for (int j = 1; j <= i + m; j++) w[j] -= 2;
}
return powerProd(w);
}

private static BigInteger countFR(int n) {
if ((n & 1) == 1) return BigInteger.ZERO;
int m = n / 2;
int[] w = new int[3*n-2];
for (int j = 1; j <= m; j++) w[j]--;
for (int j = 2*m; j <= 3*m-1; j++) w[j]++;
for (int i = 0; i <= 2*m-3; i++) {
for (int j = i + 2*m + 1; j <= i + 4*m; j++) w[j]++;
for (int j = 2*i + 3; j <= 2*i + 2*m + 2; j++) w[j]--;
}
return powerProd(w);
}

private static BigInteger powerProd(int[] w) {
BigInteger result = BigInteger.ONE;
for (int x = w.length - 1; x > 1; x--) {
if (w[x] == 0) continue;

int p = minp[x];
if (p == x) result = result.multiply(BigInteger.valueOf(p).pow(w[p]));
else {
// Redistribute it. This should ensure we avoid negatives.
w[p] += w[x];
w[x / p] += w[x];
}
}

return result;
}
}

• You are truly a god amongst mortals. I hope to see your solution published in a prestigious journal. – Alex A. Jun 12 '15 at 23:28
• This is awesome. BTW my (currently not posted) code gives 22306956 for N=5: 22231176(12)+275(4)+75328(6)+352(2), a discrepancy of 1, which is odd. I've no idea you're doing here, is it suitable for a a breakdown by symmetries? For N=4 I'm 16 lower than you and oeis.org/A066931/a066931.txt From that reference it seems I have 16 too many of multiplicity 12, which I need to convert to 32 of multiplicity 6. I'm not too surprised, even N are more difficult for me. But I have no problems with odd N and I get right answers for 0<N<4. Will look for obvious problems and post my code tomorrow. – Level River St Jun 13 '15 at 2:14
• @steveverrill, if I understand the notation, for N=5 I make it 22231176(12) + 75328(6) + 275(4) + 176(2). I think you're failing to divide the index 2 ones by 2. (FWIW for odd numbers they all have an axis of symmetry passing through two vertices and rotational symmetry of order 3). – Peter Taylor Jun 13 '15 at 6:22
• @steveverrill, and for N=4 your discrepancy seems to be a perfect fit for the number which have an axis of symmetry passing through the midpoints of two edges. – Peter Taylor Jun 13 '15 at 6:32
• Impressive that you solved this. I'm hoping you will eventually post an answer that non-mathematicians can follow. – DavidC Jun 13 '15 at 13:26

# C

Introduction

As commented by David Carraher, the simplest way of analysing the hexagon tiling seemed to be to take advantage of its isomorphism with the 3 dimensional Young Diagram, essentially an x,y square filled with integer height bars whose z heights must stay the same or increase as the z axis is approached.

I started with an algorithm for finding the totals which is more amenable to adaptation for symmetry counting than the published algorithm, which is based on a bias to one of the three cartesian axes.

Algorithm

I start by filling the cells of the x,y, and z planes with 1's, while the rest of the area contains zeroes. Once that is done, I build up the pattern layer by layer, with each layer containing the cells which have a common 3D manhattan distance from the origin. A cell can only contain a 1 if the three cells below it also contain a 1. if any of them contains a 0, then the cell must be a 0.

The advantage of building the pattern up in this way is that each layer is symmetrical about the x=y=z line. This means that each layer can be checked independently for symmetry.

Symmetry checking

The symmetries of the solid are as follows: 3 fold rotation about the x=y=z line --> 3 fold rotation about the hexagon centre; and 3 x reflections about the 3 planes containing the x=y=z line and each of the axes x,y,z --> reflection about the lines through the hexagon corners.

This adds up only to 6 fold symmetry. To get the full symmetry of the hexagon, another kind of symmetry must be considered. Each solid (built from 1's) has a complementary solid (built from 0's). Where N is odd, the complementary solid must be different from the original solid (because it is not possible for them to have the same number of cubes). Yet when the complementary solid is turned round, it will be found that its 2D representation as a diamond tiling is identical (except for a 2-fold symmetry operation) to the original solid. Where N is even, it is possible for the solid to be self-inverse.

This can be seen in the examples for N=2 in the question. If viewed from the left, the first hexagon looks like a solid cube with 8 small cubes, while the last hexagon looks like an empty shell with 0 small cubes. If viewed from the right, the reverse is true. The 3rd, 4th and 5th hexagons and the 16th, 17th and 18th hexagons look like they contain either 2 or 6 cubes, and thus they complement each other in 3 dimensions. They are related to each other in 2 dimensions by a 2-fold symmetry operation (2-fold rotation, or reflection about an axis through the hexagon's edges.) On the other hand the 9th,10th,11th and 12th hexagons show 3D patterns which are their own complements, and therefore have a higher symmetry (these are therefore the only patterns with odd multiplicity).

Note that having (N^3)/2 cubes is a necessary condition to be self complement, but in general it is not a sufficient condition if N>2. The result of all this is that for odd N, tilings always occur in pairs (N^3)/2 cubes must be carefully inspected.

Current code (generates the right total for N=1,2,3,5. Error as discussed for N=4.)

int n;                     //side length

char t[11][11][11];        //grid sized for N up to 10

int q[29][192], r[29];     //tables of coordinates for up to 10*3-2=28 layers

int c[9];                  //counts arrangements found by symmetry class. c[8] contains total.

//recursive layer counting function. m= manhattan distance, e= number of cells in previous layers, s=symmetry class.
void f(int m,int e,int s){

int u[64], v[64], w[64]; //shortlists for x,y,z coordinates of cells in this layer
int j=0;
int x,y,z;

for (int i=r[m]*3; i; i-=3){
// get a set of coordinates for a cell in the current layer.
x=q[m][i-3]; y= q[m][i-2]; z= q[m][i-1];
// if the three cells in the previous layer are filled, add it to the shortlist u[],v[],w[]. j indicates the length of the shortlist.
if (t[x][y][z-1] && t[x][y-1][z] && t[x-1][y][z]) u[j]=x, v[j]=y, w[j++]=z ;
}

// there are 1<<j possible arrangements for this layer.
for (int i = 1 << j; i--;) {

int d = 0;

// for each value of i, set the 1's bits of t[] to the 1's bits of i. Count the number of 1's into d as we go.
for (int k = j; k--;) d+=(t[u[k]][v[k]][w[k]]=(i>>k)&1);

// we have no interest in i=0 as it is the empty layer and therefore the same as the previous recursion step.
// Still we loop through it to ensure t[] is properly cleared.

if(i>0){
int s1=s;    //local copy of symmetry class. 1's bit for 3 fold rotation, 2's bit for reflection in y axis.
int sc=0;    //symmetry of self-complement.

//if previous layers were symmetrical, test if the symmetry has been reduced by the current layer
if (s1) for (int k = j; k--;) s1 &= (t[u[k]][v[k]][w[k]]==t[w[k]][u[k]][v[k]]) | (t[u[k]][v[k]][w[k]]==t[w[k]][v[k]][u[k]])<<1;

//if exactly half the cells are filled, test for self complement
if ((e+d)*2==n*n*n){
sc=1;
for(int A=1; A<=(n>>1); A++)for(int B=1; B<=n; B++)for(int C=1; C<=n; C++) sc&=t[A][B][C]^t[n+1-A][n+1-B][n+1-C];
}

//increment counters for total and for symmetry class.
c[8]++; c[s1+(sc<<2)]++;

//uncomment for graphic display of each block stacking with metadata. not recommended for n>3.
//printf("m=%d  j=%d  i=%d c1=%d-2*%d=%d c3=%d cy=%d(cs=%d) c3v=%d ctot=%d\n",m,j,i,c[0],c[2],c[0]-2*c[2],c[1],c[2],c[2]*3,c[3],c[8]);
//printf("m=%d  j=%d  i=%d C1=%d-2*%d=%d C3=%d CY=%d(CS=%d) C3V=%d ctot=%d\n",m,j,i,c[4],c[6],c[4]-2*c[6],c[5],c[6],c[6]*3,c[7],c[8]);
//for (int A = 0; A<4; A++, puts(""))for (int B = 0; B<4; B++, printf(" "))for (int C = 0; C<4; C++) printf("%c",34+t[A][B][C]);

//recurse to next level.
if(m<n*3-2)f(m + 1,e+d,s1);

}
}
}

main()
{
scanf("%d",&n);

int x,y,z;

// Fill x,y and z planes of t[] with 1's
for (int a=0; a<9; a++) for (int b=0; b<9; b++) t[a][b][0]= t[0][a][b]= t[b][0][a]= 1;

// Build table of coordinates for each manhattan layer
for (int m=1; m < n*3-1; m++){
printf("m=%d : ",m);
int j=0;
for (x = 1; x <= n; x++) for (y = 1; y <= n; y++) {
z=m+2-x-y;
if (z>0 && z <= n) q[m][j++] = x, q[m][j++] = y, q[m][j++]=z, printf(" %d%d%d ",x,y,z);
r[m]=j/3;
}
printf(" : r=%d\n",r[m]);
}

// Set count to 1 representing the empty box (symmetry c3v)
c[8]=1; c[3]=1;

// Start searching at f=1, with 0 cells occupied and symmetry 3=c3v
f(1,0,3);

// c[2 and 6] only contain reflections in y axis, therefore must be multiplied by 3.
// Similarly the reflections in x and z axis must be subtracted from c[0] and c[4].
c[0]-=c[2]*2; c[2]*=3;
c[4]-=c[6]*2; c[6]*=3;

int cr[9];cr[8]=0;
printf("non self-complement                   self-complement\n");
printf("c1  %9d/12=%9d           C1  %9d/6=%9d\n",   c[0], cr[0]=c[0]/12,     c[4], cr[4]=c[4]/6);
if(cr[0]*12!=c[0])puts("c1 division error");if(cr[4]*6!=c[4])puts("C1 division error");

printf("c3  %9d/4 =%9d           C3  %9d/2=%9d\n",   c[1], cr[1]=c[1]/4,      c[5], cr[5]=c[5]/2);
if(cr[1]*4!=c[1])puts("c3 division error");if(cr[5]*2!=c[5])puts("C3 division error");

printf("cs  %9d/6 =%9d           CS  %9d/3=%9d\n",   c[2], cr[2]=c[2]/6,      c[6], cr[6]=c[6]/3);
if(cr[2]*6!=c[2])puts("cs division error");if(cr[6]*3!=c[6])puts("CS division error");

printf("c3v %9d/2 =%9d           C3V %9d/1=%9d\n",   c[3], cr[3]=c[3]/2,      c[7], cr[7]=c[7]);
if(cr[3]*2!=c[3])puts("c3v division error");

for(int i=8;i--;)cr[8]+=cr[i];
printf("total =%d unique =%d",c[8],cr[8]);
}


Output

The program generates an output table of 8 entries, in accordance with the 8 symmetries of the solid. The solid can have any of 4 symmetries as follows (Schoenflies notation)

c1: no symmetry
c3: 3-fold axis of rotation (produces 3-fold axis of rotation in hexagon tiling)
cs: plane of reflection (produces line of reflection in hexagon tiling)
c3v both of the above (produces 3-fold axis of rotation and three lines of reflection through the hexagon corners)


Additionally, when the solid has exactly half the cells with 1's and half with 0's, there exists the possibility of flipping all the 1's and 0's, then inverting the coordinates through the centre of the cube space. This is what I am calling self-complement, but a more mathematical term would be "antisymmetric with respect to a centre of inversion."

This symmetry operation gives a 2-fold axis of rotation in the hexagon tiling.

Patterns which have this symmetry are listed in a separate column. They only occur where N is even.

My count seems to be slightly off for N=4. In discussion with Peter Taylor it appears that I am not detecting tilings which only have a symmetry of a line through the hexagon edges. This is presumably because I have not tested for self complement (antisymmetry) for operations other than (inversion)x(identity.) Testing for self complement for the operatons (inversion)x(reflection) and (inversion)x(3-fold rotation) may uncover the missing symmetries. I would then expect the first line of the data for N=4 to look like this (16 less in c1 and 32 more in C1):

c1   224064/12=18672          C1  534/6=89


This would bring the totals in line with Peter's answer and https://oeis.org/A066931/a066931.txt

current output is as follows.

N=1
non self-complement     self-complement
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs      0/6 = 0           CS  0/3= 0
c3v     2/2 = 1           C3V 0/1= 0
total =2 unique =1

non self-complement     self-complement
N=2
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs     12/6 = 2           CS  3/3= 1
c3v     4/2 = 2           C3V 1/1= 1
total =20 unique =6

N=3
non self-complement     self-complement
c1    672/12=56           C1  0/6= 0
c3      4/4 = 1           C3  0/2= 0
cs    288/6 =48           CS  0/3= 0
c3v    16/2 = 8           C3V 0/1= 0
total =980 unique =113

N=4 (errors as discussed)
non self-complement     self-complement
c1   224256/12=18688          C1  342/6=57
c3       64/4 =16             C3  2/2= 1
cs     8064/6 =1344           CS  54/3=18
c3v      64/2 =32             C3V 2/1= 2
total =232848 unique =20158

N=5
non self-complement     self-complement
c1  266774112/12=22231176        C1  0/6= 0
c3       1100/4 =275             C3  0/2= 0
cs     451968/6 =75328           CS  0/3= 0
c3v       352/2 =176             C3V 0/1= 0
total =267227532 unique =22306955


To-do list (updated)

Tidy up current code.

Done, more or less

Implement symmetry checking for current layer, and pass a parameter for the symmetry of the previous layer (no point in checking if the last layer was asymmetrical.)

Done, results for odd N agree with published data

Add an option to suppress counting asymmetrical figures (should run much faster)

This can be done by adding another condition to the recursion call: if(s1 && m<n*3-2)f(m + 1,e+d,s1) It reduces the run time for N=5 from 5 minutes down to about a second. As a result the first line of the output becomes total garbage (as do the overall totals) but if the total is already known from OEIS the number of asymmetric tilings can be reconstituted, at least for odd N.

But for even N, the number of asymmetric (according to c3v symmetries) solids that are self-complement would be lost. For this case, a separate program dedicated to solids with exactly (N**3)/2 cells with a 1 may be useful. With this available (and counting correctly) it may be possible to try N=6, but it will take a long time to run.

Implement counting of cells to reduce search to up to (N^3)/2 cubes.

Not done, savings expected to be marginal

Implement symmetry (complementary solid) checking for patterns containing exactly (N^3)/2 cubes.

Done, but seems to have omissions, see N=4.

Find a way to pick the lexically lowest figure from an asymetrical one.

Savings are not expected to be that great. Supressing asymmetrical figures eliminates most of this. The only reflection that is checked is the plane through the y axis (x and z are calculated later by multiplying by 3.) Figures with only rotational symmetry are counted in both their enantiomeric forms. Perhaps it would run nearly twice as fast if only one was counted.

To facilitate this, possibly improve the way the coordinates in each layer are listed (they form degenerate groups of 6 or 3, with possibly a group of 1 in the exact centre of the layer.)

Interesting but probably there are other questions on the site to explore.