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 &lt; 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 &gt; 0) {
            for (String arg : args) System.out.println(count(Integer.parseInt(arg)));
        }
        else {
            for (int n = 0; n &lt; 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 &lt; 3*n) {
            int[] wider = new int[3*n];
            System.arraycopy(minp, 0, wider, 0, minp.length);
            for (int x = minp.length; x &lt; 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 &lt; n; i++) {
            for (int j = i + 1; j &lt;= i + n; j++) w[j]--;
            for (int j = i + n + 1; j &lt;= 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 &lt; n; i++) {
            w[3*i+2]++;
            for (int j = 3*i + 1; j &lt;= 2*i + n + 1; j++) w[j]--;
            for (int j = 2*i + n + 1; j &lt;= 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 &lt; n; i++) {
            for (int j = 2*i + 1; j &lt;= 2*i + n; j++) w[j]--;
            for (int j = i + n + 1; j &lt;= i + 2*n; j++) w[j]++;
        }
        return powerProd(w);
    }

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

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

    private static BigInteger countFR(int n) {
        if ((n &amp; 1) == 1) return BigInteger.ZERO;
        int m = n / 2;
        int[] w = new int[3*n-2];
        for (int j = 1; j &lt;= m; j++) w[j]--;
        for (int j = 2*m; j &lt;= 3*m-1; j++) w[j]++;
        for (int i = 0; i &lt;= 2*m-3; i++) {
            for (int j = i + 2*m + 1; j &lt;= i + 4*m; j++) w[j]++;
            for (int j = 2*i + 3; j &lt;= 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 &gt; 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;
    }
}