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