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 fastest-code. 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.
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;
    }
}