// 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;
}