// A319900 - a(n) is the number of distinct ways to arrange n copies of each of the numbers 1 through n^2 // inside a fixed n X n X n cube, provided that no number appears twice in the same left-right plane, // front-back plane, or top-bottom plane. // C++ program for a(4) by Bert Dobbelaere // This program was specifically written to find a(4) of the sequence but can reproduce smaller terms. // Consider a (n x n x n) cube and a set of n positions that all have unique x, y and z values within the set. // (coordinates are integers in range [0..n-1]) // By permuting the constant y planes we can obtain unique coordinates per dimension of the form (x,x,z) // By permuting the constant z planes next, can can obtain unique coordinates per dimension of the form (x,x,x) // All (n!)^2 possible sets can be transformed to this form using a unique permutation. // We have (n^2)! possibilities to assign numbers to the x=0 plane, of which (n^2-1)! independent from the permutation. // By demanding that the numbers of the form (x,x,x) are identical and freezing the number assignment, the number of solutions // we find under these restrictions needs to be multiplied by (n^2)! * (n-1)!^2 // The program below counts the restricted solutions i.e. the number of ways to partition the cube into n^2 sets with unique // element coordinates per dimension. Each set is precomputed and represented by a bit pattern of which n bits will be set. // Indicative run time: ~30' on Intel Core i7 4770K. // // a(1) = 1 * 1 = 1 // a(2) = 1 * 24 = 24 // a(3) = 10 * 1451520 = 14515200 // a(4) = 10234349568 * 753220435968000 = 7708721243457872461824000 #include <bits/stdc++.h> // Need a 64 bit unsigned integer #if __cplusplus < 201103L typedef unsigned long long u64; #else typedef uint64_t u64; #endif typedef u64 PositionSet; // Can hold any subset of 4x4x4 positions #define N 4 // Configure "n" #if N==1 #define NPATTERNS 1 // Square of (1-1)! #elif N==2 #define NPATTERNS 1 // Square of (2-1)! #elif N==3 #define NPATTERNS 4 // Square of (3-1)! #elif N==4 #define NPATTERNS 36 // Square of (4-1)! #else #error "n>4 not supported" #endif // Cube position set operation macros. p and q are sets, (x,y,z) coordinates defining an element #define DISJOINT(p,q) (((p)&(q))==0) #define UNION(p,q) ((p)|(q)) #define SINGLETON(x, y, z) (1ULL << ((x)*N*N + (y)*N + (z))) #define CONTAINS(x, y, z, p) (!DISJOINT(p, SINGLETON(x,y,z))) // Each pattern contains a conflict-free set of N positions (represented by single bits). // There are in total (N!)^2 such sets in the cube, here organised by the element they contain from plane x=0 (first index in patterns array) PositionSet patterns[N*N][NPATTERNS]; u64 count; int initcnt; // Recursive helper function for initpatterns void completePatterns(int a, PositionSet p_in, int px) { if(px==N) // We have a set of N conflict free positions contained in p_in { assert(initcnt<NPATTERNS); patterns[a][initcnt++]=p_in; return; } for(int py=0;py<N;py++) for(int pz=0;pz<N;pz++) { bool ok=true; for(int x=0;x<px;x++) for(int y=0;y<N;y++) for(int z=0;z<N;z++) if( CONTAINS(x,y,z, p_in) ) { if(pz == z) ok=false; if(py == y) ok=false; } // Fill next plane and go one level deeper if(ok) completePatterns(a, UNION(p_in, SINGLETON(px, py, pz) ),px+1); } } // Initialize the patterns array void initpatterns() { assert( (8*sizeof(PositionSet)) >= (N*N*N) ); for(int y=0;y<N;y++) for(int z=0;z<N;z++) { int k=N*y+z; initcnt=0; // Collect the N element sets that share the same starting element in plane x=0 completePatterns(k, SINGLETON(0,y,z), 1); assert(initcnt==NPATTERNS); } } // Recursive core routine. Counts the number of ways N^2 disjoint sets of N "conflict free" elements each can fill the cube. void solve(int lvl, PositionSet p_in) { if(lvl==N*N) { count++; return; } for(int k=0;k<NPATTERNS;k++) if(DISJOINT(p_in, patterns[lvl][k])) { solve(lvl+1, UNION(p_in, patterns[lvl][k])); } } int main() { initpatterns(); count=0; solve(1,patterns[0][0]); // Start with set { (k,k,k) for k=0..(N-1) } std::cout << "a(" << N << ") = " << count << " * " << (N*N) << "! * (" << (N-1) << "!)^2" << std::endl; return 0; }