using System.Diagnostics;
internal class Program
{
    const int n = 16;//number of polycube cells. Need n >= 4 if single threading, or n >= filterDepth >= 5 if multithreading (I think)
    const int filterDepth = 12;
    const int threads = 5;
    private static void Main()
    {
        long count = 0;//enumerate the sum over the order 24 group of the size of the fix of each group element, and divide by 24 (Burnside's lemma)
        {
            Stopwatch sw = Stopwatch.StartNew();
            string[] descriptions = new string[4] { "orthogonal order 2 rotation", "orthogonal order 4 rotation", "short diagonal order 2 rotation", "long diagonal order 3 rotation" };
            int[] autClassSizes = new int[4] { 3, 6, 6, 8 };
            int[][] matrixReps = new int[4][] { new int[9] { -1, 0, 0, 0, -1, 0, 0, 0, 1 }, new int[9] { 0, -1, 0, 1, 0, 0, 0, 0, 1 }, new int[9] { 0, 1, 0, 1, 0, 0, 0, 0, -1 }, new int[9] { 0, 0, 1, 1, 0, 0, 0, 1, 0 } };
            int[][] affine1 = new int[4][] { new int[3] { 1, 0, 0 }, new int[3] { 1, 0, 0 }, new int[3] { 1, -1, 0 }, new int[3] { 1, 0, -1 } };
            int[][] affine2 = new int[4][] { new int[3] { 0, 1, 0 }, new int[3] { 0, 1, 0 }, new int[3] { 0, 0, 1 }, new int[3] { 0, 1, -1 } };
            int[][] biases = new int[4][] { new int[3] { 2 * n, 2 * n, 0 }, new int[3] { 2 * n, 0, 0 }, new int[3] { 0, 0, 2 }, new int[3] { n - 1, 0, 1 - n } };
            for (int sym = 0; sym < 4; sym++)
            {
                long subcount = 0;
                for (int i = 1 - n; i <= n - 1; i++)
                {
                    for (int j = 1 - n + Math.Abs(i); j <= n - 1 - Math.Abs(i); j++)
                    {
                        subcount += CountSymmetricPolycubes(matrixReps[sym], new int[3] { i * affine1[sym][0] + j * affine2[sym][0] + biases[sym][0], i * affine1[sym][1] + j * affine2[sym][1] + biases[sym][1], i * affine1[sym][2] + j * affine2[sym][2] + biases[sym][2] });
                    }
                }
                Console.WriteLine($"{subcount:N0} polycubes fixed under each {descriptions[sym]}");
                Console.WriteLine($"{autClassSizes[sym] * subcount:N0} in total (symmetry occurs {autClassSizes[sym]} times)");
                count += autClassSizes[sym] * subcount;
            }
            Console.WriteLine("total count for nontrivial symmetries is {0:N0} for polycubes with {1} cells\nTook {2:%d} day(s), {2:hh\\:mm\\:ss\\.f} h:m:s", count, n, sw.Elapsed);
        }
        Console.WriteLine();
        {
            Stopwatch sw = Stopwatch.StartNew();
            long subcount = 0;
            Task<long>[] tasks = new Task<long>[threads];//please don't judge my multithreading - this was basically just my first attempt, to see if it helped
            int i = 4 * (n - filterDepth) - 2;//this is the maximum left stack length, which is the value being filtered to separate work
            for (int j = 0; j < tasks.Length; j++)
            {
                int filter = i--;//copy of i, since lambda expression captures the variable
                tasks[j] = Task<long>.Factory.StartNew(() => CountExtensionsSubset(filter));
            }
            while (i >= 0)
            {
                int j = 0;
                while (!tasks[j].IsCompleted)
                {
                    j++;
                    j %= tasks.Length;
                    Thread.Sleep(50);
                }
                subcount += tasks[j].Result;
                int filter = i--;//copy of i, since lambda expression captures the variable
                tasks[j] = Task<long>.Factory.StartNew(() => CountExtensionsSubset(filter));
            }
            subcount += tasks.Sum(t => t.Result);
            count += subcount;
            Console.WriteLine("{0:N0} polycubes with {1} cells (number of polycubes fixed by trivial symmetry)\nTook {2:%d} day(s), {2:hh\\:mm\\:ss\\.f} h:m:s", subcount, n, sw.Elapsed);
        }
        Console.WriteLine();
        count /= 24;
        Console.WriteLine("{0:N0} free polycubes with {1} cells", count, n);
    }
    private static long CountSymmetricPolycubes(int[] linearMap, int[] affineShift)
    {
        byte[,,] adjacencyCounts = new byte[2 * n + 1, 2 * n + 1, n + 2];
        for (int z = 0; z < 2; z++)
        {
            for (int y = 0; y < 2 * n + 1; y++)
            {
                for (int x = 0; x < 2 * n + 1; x++)
                {
                    adjacencyCounts[x, y, z] = 1;
                    if (z == 1 && y == n && x == n) goto loopExit;
                }
            }
        }
    loopExit:
        HashSet<(int, int, int)> requiredCells = new HashSet<(int, int, int)>();
        Stack<(int, int, int)> extensionStack = new Stack<(int, int, int)>();
        extensionStack.Push((n, n, 1));
        Stack<(int, int, int)> recoveryStack = new Stack<(int, int, int)>();
        return CountExtensions(n);

        long CountExtensions(int cellsToAdd)
        {
            cellsToAdd--;
            long count = 0;
            int originalLength = extensionStack.Count;
            while (extensionStack.Count > 0)
            {
                int x, y, z;
                recoveryStack.Push((x, y, z) = extensionStack.Pop());

                bool existingRequirement = requiredCells.Remove((x, y, z));
                if (!existingRequirement)
                {
                    if (cellsToAdd < requiredCells.Count) continue;//number of required cells will only grow, so if already impossible, go to next
                    int tempX = x;
                    int tempY = y;
                    int tempZ = z;
                    while (true)//works for general transformations of finite order
                    {
                        (tempX, tempY, tempZ) = (linearMap[0] * tempX + linearMap[1] * tempY + linearMap[2] * tempZ + affineShift[0],
                                                 linearMap[3] * tempX + linearMap[4] * tempY + linearMap[5] * tempZ + affineShift[1],
                                                 linearMap[6] * tempX + linearMap[7] * tempY + linearMap[8] * tempZ + affineShift[2]);
                        if (x == tempX && y == tempY && z == tempZ) break;
                        requiredCells.Add((tempX, tempY, tempZ));
                    }
                }

                if (cellsToAdd >= requiredCells.Count)//if there are too many required cells, then no valid polycubes are possible
                {
                    if (cellsToAdd == 0) count++;
                    else
                    {
                        int innerOriginalLength = extensionStack.Count;
                        if (adjacencyCounts[x - 1, y, z]++ == 0) extensionStack.Push((x - 1, y, z));
                        if (adjacencyCounts[x, y - 1, z]++ == 0) extensionStack.Push((x, y - 1, z));
                        if (adjacencyCounts[x, y, z - 1]++ == 0) extensionStack.Push((x, y, z - 1));
                        if (adjacencyCounts[x + 1, y, z]++ == 0) extensionStack.Push((x + 1, y, z));
                        if (adjacencyCounts[x, y + 1, z]++ == 0) extensionStack.Push((x, y + 1, z));
                        if (adjacencyCounts[x, y, z + 1]++ == 0) extensionStack.Push((x, y, z + 1));

                        count += CountExtensions(cellsToAdd);

                        --adjacencyCounts[x - 1, y, z];
                        --adjacencyCounts[x, y - 1, z];
                        --adjacencyCounts[x, y, z - 1];
                        --adjacencyCounts[x + 1, y, z];
                        --adjacencyCounts[x, y + 1, z];
                        --adjacencyCounts[x, y, z + 1];
                        while (extensionStack.Count != innerOriginalLength) extensionStack.Pop();//should replace this with custom stack to avoid this unnecessary loop
                    }
                }

                if (existingRequirement)
                {
                    requiredCells.Add((x, y, z));
                    break;//this required cell will no longer be available in the extension stack, so no more valid polycubes are possible in this branch
                }
                else
                {
                    int tempX = x;
                    int tempY = y;
                    int tempZ = z;
                    while (true)
                    {
                        (tempX, tempY, tempZ) = (linearMap[0] * tempX + linearMap[1] * tempY + linearMap[2] * tempZ + affineShift[0],
                                                 linearMap[3] * tempX + linearMap[4] * tempY + linearMap[5] * tempZ + affineShift[1],
                                                 linearMap[6] * tempX + linearMap[7] * tempY + linearMap[8] * tempZ + affineShift[2]);
                        if (x == tempX && y == tempY && z == tempZ) break;
                        requiredCells.Remove((tempX, tempY, tempZ));
                    }
                }
            }
            while (extensionStack.Count != originalLength) extensionStack.Push(recoveryStack.Pop());
            return count;
        }
    }
    private static long CountExtensionsSubset(int filter)
    {
        unsafe
        {
            const int X = (n + 5) / 4 * ((n + 5) / 4 * 3 - 2);//set X<Y<Z such that aX+bY+cZ = 0 implies a = b = c = 0 or |a|+|b|+|c| > n
            const int Y = X + 1;//trivial choice is X = 1, Y = n, Z = n * n. A simple reduction is X = 1, Y = n, Z = n * (n / 2) + (n + 1) / 2
            const int Z = X + (n + 5) / 4 * 3;//minimising Z is memory efficient. Unclear if this noticably affects performance. Z ~ 3/16 n^2 is the best I can find for arbitrary n
            byte* byteBoard = stackalloc byte[(n + 2) * Z];//could use ints or shorts as offsets to save memory, but it's faster to directly store the pointers to avoid adding pointer offsets at every lookup
            byte** refStack = stackalloc byte*[(n - 2) * 4];//total length of the two stacks is at most 4n-9. One stack grows from the left, the other stack grows from the right
            *refStack = byteBoard += Z;//seeded with first index of the byte board as the only allowed extension
            for (byte* i = byteBoard + (n + 1) * Z; --i != byteBoard;) *i = 255;//the first Z + 1 bytes are disallowed extensions; first Z are less than the minimum, last 1 due to edge case of initial polycube having no neighbours

            Console.WriteLine($"started task {filter}");
            long count = CountExtensions(n, refStack + 1, refStack + (n - 2) * 4);
            Console.WriteLine($"finished task {filter} with subcount {count:N0}");
            return count;

            long CountExtensions(int depth, byte** stackTop1, byte** stackTop2)
            {
                long count = 0;
                byte** stackTopOriginal = stackTop1;
                while (stackTop1 != refStack)
                {
                    byte* index = *--stackTop1;
                    byte** stackTopInner = stackTop1;
                    if (++*(index - X) == 0) *stackTopInner++ = index - X;
                    if (++*(index - Y) == 0) *stackTopInner++ = index - Y;
                    if (++*(index - Z) == 0) *stackTopInner++ = index - Z;
                    if (++*(index + X) == 0) *stackTopInner++ = index + X;
                    if (++*(index + Y) == 0) *stackTopInner++ = index + Y;
                    if (++*(index + Z) == 0) *stackTopInner++ = index + Z;
                    if (depth == 4) count += CountFinalExtensions(stackTopInner);
                    else if (depth != filterDepth || stackTop1 - refStack == filter) count += CountExtensions(depth - 1, stackTopInner, stackTop2);//if multithreading is not wanted, remove "if (condition)" from this else statement
                    --*(index - X);
                    --*(index - Y);
                    --*(index - Z);
                    --*(index + X);
                    --*(index + Y);
                    --*(index + Z);
                    *--stackTop2 = index;//doing this push before the recursion would add one extra unnecessary element to the stack at each level of recursion
                }
                while (stackTop1 != stackTopOriginal) *stackTop1++ = *stackTop2++;
                return count;
            }

            int CountFinalExtensions(byte** stackTop)
            {
                const int X2 = X + X, Y2 = Y + Y, Z2 = Z + Z, sYX = Y + X, sZX = Z + X, sZY = Z + Y, dYX = Y - X, dZX = Z - X, dZY = Z - Y;
                int length = (int)(stackTop - refStack);
                int count = length * (length - 1) * (length - 2) / 6;
                byte** stackTopTemp = stackTop;
                while (stackTopTemp != refStack)
                {
                    byte* i = *--stackTopTemp;
                    int neighbours = 0;
                    int subcount = 128;
                    if (*(i - X) > 127) { neighbours++; subcount += *(i - X2) + *(i - sYX) + *(i - sZX) + *(i + dYX) + *(i + dZX); count += --*(i - X); }
                    if (*(i - Y) > 127) { neighbours++; subcount += *(i - Y2) + *(i - sYX) + *(i - sZY) + *(i - dYX) + *(i + dZY); count += --*(i - Y); }
                    if (*(i - Z) > 127) { neighbours++; subcount += *(i - Z2) + *(i - sZX) + *(i - sZY) + *(i - dZX) + *(i - dZY); count += --*(i - Z); }
                    if (*(i + X) > 127) { neighbours++; subcount += *(i + X2) + *(i + sYX) + *(i + sZX) + *(i - dYX) + *(i - dZX); count += --*(i + X); }
                    if (*(i + Y) > 127) { neighbours++; subcount += *(i + Y2) + *(i + sYX) + *(i + sZY) + *(i + dYX) + *(i - dZY); count += --*(i + Y); }
                    if (*(i + Z) > 127) { neighbours++; subcount += *(i + Z2) + *(i + sZX) + *(i + sZY) + *(i + dZX) + *(i + dZY); count += --*(i + Z); }
                    count += (subcount >> 8) + (neighbours * (neighbours + ((length << 1) - 511)) >> 1);
                }
                while (stackTop != refStack)
                {
                    byte* i = *--stackTop;
                    *(i - X) |= (byte)(*(i - X) >> 4);
                    *(i - Y) |= (byte)(*(i - Y) >> 4);
                    *(i - Z) |= (byte)(*(i - Z) >> 4);
                    *(i + X) |= (byte)(*(i + X) >> 4);
                    *(i + Y) |= (byte)(*(i + Y) >> 4);
                    *(i + Z) |= (byte)(*(i + Z) >> 4);
                }
                return count;
            }
        }
    }
}

//written by Stanley Dodds (that's me!)
//I can probably supply a full proof of correctness if needed, but I have nothing written down right now - so currently, source: dude, just trust me (or try reading my code)
//Of course, you can compare up to n = 19 with previous results and see that it's definitely doing something right