[Question]
  [
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:
![](https://i.stack.imgur.com/vWbjJ.png)
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](https://oeis.org/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](https://oeis.org/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](/questions/tagged/fastest-code "show questions tagged '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.
  
[Answer]
## 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:
![Images of the 10 symmetries](https://i.stack.imgur.com/hPWBZ.png)
The number of diamond tilings of a subset of a triangular lattice can be [calculated as a determinant](https://en.wikipedia.org/wiki/FKT_algorithm), 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](https://en.wikipedia.org/wiki/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](http://cheddarmonk.org/papers/distinct-dimer-hex-tilings.pdf) 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;
    }
}
```
[Answer]
# 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.
[Answer]
# Rust, 0.03 s for n=0..=15
Port of [@Peter Taylor's Java answer](https://codegolf.stackexchange.com/a/51662/110802) in Rust.
Run it on [RustPlayground](https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=d9b5191802f5b1785460efb68e8eecb5)!
---
### `src/main.rs`
```
use std::collections::HashMap;
use num_bigint::BigUint;
use num_traits::{Zero, One, Pow};
fn main() {
    for n in 0..16 {
        println!("{}\t{}", n, count(n));
    }
}
fn count(n: usize) -> BigUint {
    if n == 0 { return BigUint::one() }
    let mut minp = vec![0; 3 * n];
    for index in 2..3 * n {
        minp[index] = (2..=index).find(|p| index % p == 0).unwrap_or(2);
    }
    let e = count_e(n, &minp);
    // println!("n={:?}  e={:?}", n, e);
    let r2 = count_r2(n, &minp);
    // println!("n={:?}  r2={:?}", n, r2);
    let f = count_f(n, &minp);
    // println!("n={:?}  f={:?}", n, f);
    let r3 = count_r3(n, &minp);
    // println!("n={:?}  r3={:?}", n, r3);
    let r = count_r(n, &minp);
    // println!("n={:?}  r={:?}", n, r);
    let fr = count_fr(n, &minp);
    // println!("n={:?}  fr={:?}", n, fr);
    let sum = e
        + r3
        + (r2 + r.clone()) * BigUint::from(2u8)
        + (f + fr) * BigUint::from(3u8);
    sum / BigUint::from(12u8)
}
fn count_e(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];
    for i in 0..n {
        for j in (i + 1)..=(i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }
    // println!("w={:?}  minp={:?}",w, minp);
    power_prod(&mut w, minp)
}
fn count_r2(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];
    for i in 0..n {
        w[3 * i + 2] += 1;
        for j in (3 * i + 1)..=(2 * i + n + 1) { w[j] -= 1 }
        for j in (2 * i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }
    power_prod(&mut w, minp)
}
fn count_f(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];
    for i in 0..n {
        for j in (2 * i + 1)..=(2 * i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }
    power_prod(&mut w, minp)
}
fn count_r3(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let res = count_e(n / 2, minp);
    res.pow(2u8)
}
fn count_r(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m = n / 2;
    // let mut w = vec![0; 3 * m - 1];
    let mut w = vec![0; 3 * m - 1];
    for i in 0..m {
        for j in 1..=(3 * i + 1) { w[j] += 2 }
        for j in 1..=(i + m) { w[j] -= 2 }
    }
    // power_prod(&w, minp)
    power_prod(&mut w, minp)
}
fn count_fr(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m:i32 = n as i32 / 2;
    let mut w = vec![0; 3 * n - 2];
    // println!("n={:?}   m={:?}", n, m);
    for j in 1..=m { w[j as usize] -= 1 }
    for j in (2 * m)..(3 * m) { w[j as usize] += 1 }
    for i in 0..=(2 * m - 3) {
        // println!("start{:?}    end{:?}", (i + 2 * m + 1),  (i + 4 * m));
        for j in (i + 2 * m + 1)..=std::cmp::min(w.len() as i32 - 1, i + 4 * m) { w[j as usize] += 1 }
        // println!("start{:?}    end{:?}", (2 * i + 3),  (2 * i + 2 * m + 2));
        for j in (2 * i + 3)..=std::cmp::min(w.len() as i32 - 1, (2 * i + 2 * m + 2)) { w[j as usize] -= 1 }
        // println!("i={:?}",i);
    }
    power_prod(&mut w, minp)
}
fn power_prod(w: &mut [i32], minp: &[usize]) -> BigUint {
    let mut result = BigUint::one();
    // println!("left={:?}   right={:?}", 2, w.len());
    for x in (2..w.len()).rev() {
        if w[x] != 0 {
            let p = minp[x];
            if p == x {
                result *= BigUint::from(p as u8).pow(w[p] as u32);
            } else {
                w[p] += w[x];
                w[x / p] += w[x];
            }
        }
        // println!("x={:?}", x);
    }
    result
}
```
### `Cargo.toml`
```
[package]
name = "rust_hello"
version = "0.1.0"
edition = "2021"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies]
num-bigint = "0.3"
num-traits = "0.2"
[profile.release]
lto = "y"
panic = "abort"
```
### Build and run
```
$ cargo build --release
$ cargo run --release
    Finished release [optimized] target(s) in 0.03s
     Running `target\release\rust_hello.exe`
0       1
1       1
2       6
3       113
4       20174
5       22306955
6       123222909271
7       3283834214485890
8       421263391026827547540
9       260028731850596651411721718
10      772086476515163830856527013278243
11      11025620741283840573496993339545350520150
12      757129347300072898736973484532998417574513923224
13      249988147151001673705019407088899436686198854480028470000
14      396838668136562149340506906750722486726178589993407032745345886384
15      3028477458789741699983305448780147795991274413836547758038383842702471625120
```
[Answer]
# Scala
A port of [@Peter Taylor's Java answer](https://codegolf.stackexchange.com/a/51662/110802) in Scala.
[TIO](https://tio.run/##zVZbb9MwFH5efsXRNJANLN1SxMNEJ22olXihqPDERShr3NYjcSLbaTJN/e3l2E7SZPQCiNtLm5yc63e@Yx81DeNwveZJlkoNyrz50zSO2VTzVPhJrsObmPlXUoZ31/lsxqTX0U1CvfC/1LLbcFmJrvn8tdBsjvpeenOL7mCcaZ5wxaJXaY6fJLBSMxEpuMoyuPcAMsmXoWawDCUkXGQwgFZcf8bjmASUnFEPlfkMSCjnyo@ZmOsFXMIZtV4ArDgJM/MdBpfGr9CxIFMT1wh9nWJylFJUXwGLFassyRmgCo/h/AX1Z6lk4XRBRNuHOj4Rn/TJvfMl6OrYefFa@UdsBtX3CzCBLmADBxblYpkC0PXA5C2ZzqVoafnjN0OvUTNg1HW@hD48AVHXCohWDAWPrOfv4HK6BjKnjDUB4SJiJbw8hbradgAKN3dwTp1Lp0rRs9Go3ryjo6NerwbEmiZf32l8n5PjZw4PC2Y7bRfIOq0D1QDDPZQIcWUF0OvBCCOBXjBQSYhcVNpgmzAoFhz1I75EPwrKxsQlaxMlAegUSvSOPkgJj@CLA9mfMz2WQ2w2kqiyXFX/FdusG@/H6lu59hj0h2hq@z1EQjTSSVCLJ0FbPqrFo452v9Hud@SNuOOjEY86cpUn@GFYVeWHUUQmfdp5DdwfxdGONc/iO9JiHfrI2XiGANGO1cj9HrDq11aYhu@atE3tPNgzMsO9M2O5fpjnzRiLNsn4hmSEcHhqWI5UsY@Cbub91sx7QW4pnA5Qp2MiumaBDbvN9Glj6jiWpQWTb2UakQJPH5v/bhAMYf4CCoU1sHW4hJtSa3lValC9uuoPI9XV/9Nojf4SZYIdoPxf1OkfvHTw1nlsMhmYOFuung/Dydgq1wMJPSSIj2m4g3NH4N8X17TMnGM28OEmJnCKDr9rZLKrkee2BxuK7@hBc0lUBkY52dXr4NfZ@@@AEwhcsAHOlpnsJ7NlfVIhbMHfz2DsBjoNqi719x7HTqszJM9thP0h2qPZ745m7TE4MKI/17aW3oUD9iN24vOu40dixxTemih6sN9VgFarUbMWSbZkUrHt65HhQ2E3Hbs4khqD1hJb9z2r17aSNnJjnhnbkm7Scg97r/bMzX@BD3Tj7UFQu7pNWMQVLks3OcLFtQ/vF1yBWqR5HAETKpe4xjEIlymPQLB5qPmSKb/lxAQx/S3aiRtxibTe8mnVWebcryvJdm/lrdffAA)
```
import scala.collection.mutable.ArrayBuffer
import scala.math._
import java.math.BigInteger
object OptimisedCounter extends App {
  private var minp = ArrayBuffer.fill(2)(0)
  if (args.length > 0) {
    args.map(arg => println(count(arg.toInt)))
  } else {
    (0 until 16).foreach(n => println(s"$n\t${count(n)}"))
  }
  private def count(n: Int): BigInteger = {
    if (n == 0) return BigInteger.ONE
    if (minp.length < 3 * n) {
      val wider = ArrayBuffer.fill(3 * n)(0)
      for (index <- 0 until (minp.length) by 1) wider(index) = minp(index)
            //println(minp.mkString(","))
      (minp.length until wider.length).foreach { x =>
        // Find the smallest prime which divides x
        wider(x) = (2 to x).find(x % _ == 0).getOrElse(2)
      }
      minp = wider
            //println(minp.mkString(","))
    }
    val E = countE(n)
    val R2 = countR2(n)
    val F = countF(n)
    val R3 = countR3(n)
    val R = countR(n)
    val FR = countFR(n)
    val sum = E
      .add(R3)
      .add(R2.add(R).multiply(BigInteger.valueOf(2)))
      .add(F.add(FR).multiply(BigInteger.valueOf(3)))
    sum.divide(BigInteger.valueOf(12))
  }
  private def countE(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((i + 1) to (i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }
  private def countR2(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      w(3 * i + 2) += 1
      ((3 * i + 1) to (2 * i + n + 1)).foreach(j => w(j) -= 1)
      ((2 * i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }
  private def countF(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((2 * i + 1) to (2 * i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }
  private def countR3(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    countE(n / 2).pow(2)
  }
  private def countR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * m - 1)(0)
    (0 until m).foreach { i =>
      (1 to (3 * i + 1)).foreach(j => w(j) += 2)
      (1 to (i + m)).foreach(j => w(j) -= 2)
    }
    powerProd(w.toArray)
  }
  private def countFR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * n - 2)(0)
    (1 to m).foreach(j => w(j) -= 1)
    (2 * m until 3 * m).foreach(j => w(j) += 1)
    (0 to 2 * m - 3).foreach { i =>
      ((i + 2 * m + 1) to (i + 4 * m)).foreach(j => w(j) += 1)
      ((2 * i + 3) to (2 * i + 2 * m + 2)).foreach(j => w(j) -= 1)
    }
    powerProd(w.toArray)
  }
  private def powerProd(w: Array[Int]): BigInteger = {
    var result = BigInteger.ONE
    (2 until w.length).reverse.foreach { x =>
      if (w(x) == 0) ()
      else {
        val 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)
        }
      }
    }
    result
  }
}
```
]