Extending OEIS: Counting Diamond Tilings
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);
BigInteger sum = E.add(R3);
sum = sum.add(R2.add(R).multiply(BigInteger.valueOf(2)));
sum = sum.add(F.add(FR).multiply(BigInteger.valueOf(3)));
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;
}
}
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.