/*
Filename:	Squares.cpp

Purpose:		To calculate all possible square tilings of a square lattice bounded by different sizes of rectangle.
				To calculate all possible cubic tilings of a cubic lattice bounded by different sizes of rectangular cuboid.
				To calculate the sets of tilings reduced for symmetry and the set of partitions.
				To determine the sets of tilings in each symmetry group.
				To count the numbers of tilings with various properties e.g.
				.	different numbers of parts
				.	different numbers of isolated nodes
				.	...

Development history:
				Developed by Chris Gribble in 2013.
*/

#include <iomanip>
#include <iostream>
#include <fstream>
#include <malloc.h>
#include <stdlib.h>
#include <string>

using namespace std;

const int maxNodes     = 400;								/* Maximum number of nodes allowed in the bounded lattice. */
const int maxNeighbors = 6;									/* Maximum number of neighbors a node can have. */
const int maxPartSets  = 200;								/* Maximum numbers of distinct sets of square parts. */
const int maxTilings   = 10000000;							/* Maximum number of tilings. */

int		count[10], zeroCount;								/* Counters. */
int		nodeNeighbors[maxNodes][maxNeighbors];				/* List of neighbors of each node in the grid. */
int		numberOfNodeNeighbors[maxNodes];					/* Number of neighbors of each node. */
int		numberOfZeros[maxNodes];							/* Histogram of number of 0s in each unrestricted tiling. */
int		numberOfZerosInequivalent[maxNodes];				/* Histogram of number of 0s in each inequivalent tiling. */
int		numberOfZerosPartSet[maxNodes];						/* Histogram of number of 0s in each part set. */
int		partSetCount[maxPartSets];							/* Number of times each set of parts occurs in the set of unrestricted tilings. */
int		partSetCountInequivalent[maxPartSets];				/* Number of times each set of parts occurs in the set of inequivalent tilings. */
int		partSets[maxPartSets][12];							/* Distinct sets of parts. */

int		a, b, c;											/* Dimensions of lattice boundary in x , y and z directions, respectively. */
int		h, i, j, k, l, m, n, s;								/* Loop counters. */
int		firstTiling, lastTiling;							/* Pointers to first and last tilings generated in an iteration in the list of tilings. */
int		numberOfTilings;									/* Number of tilings for display. */
int		tilingIndex;										/* Index to list of tilings. */

short	parts[maxNodes];									/* Number of squares of each size in a tiling. */
short	rotate090[maxNodes];								/* Rotation of nodal array by 90 degs. */
short	rotate180[maxNodes];								/* Rotation of nodal array by 180 degs. */
short	rotate270[maxNodes];								/* Rotation of nodal array by 270 degs. */
short	reflectHoriz[maxNodes];								/* Reflection of nodal array about a horizontal axis. */
short	reflectVert[maxNodes];								/* Reflection of nodal array about a vertical axis. */
short	refHorizRot90[maxNodes];							/* Reflection of nodal array about a horizontal axis then a rotation by 90 degs. */
short	refVertRot90[maxNodes];								/* Reflection of nodal array about a vertical axis then a rotation by 90 degs. */
short	sum;												/* Accumulative sum of terms. */
short	tilingInd[maxTilings];								/* -1 = Tiling is not distinct, i.e. is a rotation or reflection of a tiling with lower index
															    0 = Tiling not yet examined
															  > 0 = Tiling is distinct and has 1, 2, 4 or 8 different transformations. */
short	tilingPartSet[maxTilings];							/* Part set index for each tiling. */

short	latticeArea;										/* Area of bounded lattice. */
short	latticeDimension;									/* User-specified lattice dimension (2 or 3). */
short	latticeVolume;										/* Volume of bounded lattice. */
short	maxCubes;											/* Maximum number of cubes that can be incorporated into tiling. */
short	maxCubeSize;										/* Maximum size of cube that will fit in rectangular cuboid. */
short	maxSquares;											/* Maximum number of squares that can be incorporated into tiling. */
short	maxSquareSize;										/* Maximum size of square that will fit in rectangle. */
short	node;												/* Node number. */
short	numberOfCubes;										/* Number of cubes to be incorporated into tiling. */
short	numberOfNodes;										/* Number of nodes in the bounded lattice. */
short	numberOfPartSets;									/* Number of distinct sets of parts. */
short	numberOfSquares;									/* Number of squares to be incorporated into tiling. */
short	partsMatch;											/* Number of parts that match in 2 tilings. */

bool    finished;											/* Indicates whether the program run is finished or not. */
bool	duplicate;											/* Indicates whether a duplicate tiling has been found or not. */
bool	stop;												/* Indicates whether nested loop operations should be terminated or not. */
bool	tiling[maxNodes];									/* Tiling consists of a value of 0 or 1 for each node:
															    0 = Node has no neighbors and is inside a square with dimension 2x2 <= square <= axa
																1 = Node has neighbors and is on the boundary of a square with dimension 1x1 <= square <= axa. */

bool	**tilingList;										/* List of tilings including different rotations and reflections. */

int main ()
{
	/*
	This program calculates all possible tilings with squares of different sizes of a square lattice bounded by a user-specified size of rectangle.
	It also calculates all possible tilings with cubes of different sizes of a cubic lattice bounded by a user-specified size of rectangular cuboid.
	In particular, it determines all equivalent and inequivalent tilings under rotation and reflection.
	It partitions equivalent and inequivalent tilings into groups of order 1, 2, 4 and 8.
	It counts the number of squares/cubes of each size in each tiling.

	In the 2D mode, the user is prompted to supply the dimensions of the rectangle, a and b.
	In the 3D mode, the user is prompted to supply the dimensions of the rectangular cuboid, a, b and c.

	The first part of the program determines the set of neighbors for each node in the relevant lattice.

	The node numbering convention for a selection of 2D lattices is:

		a = b = 2		0  1
						2  3

		a = 2, b = 3	0  1
						2  3
						4  5

		a = 3, b = 2	0  1  2
						3  4  5

		a = b = 3		0  1  2
						3  4  5
						6  7  8

		a = b = 4		 0  1  2  3
						 4  5  6  7
						 8  9 10 11
						12 13 14 15

		a = b = 5		 0  1  2  3  4
						 5  6  7  8  9
						10 11 12 13 14
						15 16 17 18 19
						20 21 22 23 24

	The node numbering convention for a selection of 3D lattices is:

		a = b = c = 2		 0  1
						   2  3
							 4  5
						   6  7

		a = b = c = 3		 0  1  2
						   3  4  5
						 6  7  8
							 9 10 11
						  12 13 14
						15 16 17
							18 19 20
						  21 22 23
						24 25 26

		a = b = c = 4		 0  1  2  3
						   4  5  6  7
						 8  9 10 11
					   12 13 14 15
							16 17 18 19
						   20 21 22 23
						  24 25 26 27
						28 29 30 31
							32 33 34 35
						   36 37 38 39
						  40 41 42 43
						 44 45 46 47
							48 49 50 51
						   52 53 54 55
						  56 57 58 59
						 60 61 62 63

	a = 5, b = 4, c = 3		  0  1  2  3  4
						    5  6  7  8  9
						  10 11 12 13 14
						15 16 17 18 19
							20 21 22 23 24
						   25 26 27 28 29
						  30 31 32 33 34
						35 36 37 38 39
							40 41 42 43 44
						   45 46 47 48 49
						  50 51 52 53 54
						 55 56 57 58 59
	*/

	/* Create output file. */

	fstream file;
	file.open("Squares.txt", ios_base::out | ios_base::trunc);

	/* Initialise variables. */

	numberOfTilings = 0;										/* Set the lnumber of tilings generated to zero. */

	/* Get lattice type. */

	cout << "Enter dimension of lattice [2 = 2D / 3 = 3D] ";
	cin  >> latticeDimension;

	if (latticeDimension != 2 && latticeDimension != 3)
		goto error;

	if (latticeDimension == 2)
	{
		/* Get dimensions of lattice. */

		cout << "Enter length of rectangle > 1, a = ";			/* Get the length dimension   a a a a */
		cin  >> a;

		if (a < 2)
			goto error;

		cout << "Enter width of rectangle > 1,  b = ";			/* Get the width dimension   b */
		cin  >> b;												/*                           b */
																/*                           b */
																/*                           b */
		if (b < 2)
			goto error;

		numberOfNodes = a * b;

		if (numberOfNodes > maxNodes)
			goto error;

		file << "Lattice has 2 dimensions" << endl;
		file << "Length of rectangle, a = " << a << endl;
		file << "Width of rectangle,  b = " << b << endl;

		/* The first part of the program creates the neighbor array.  This part of the program has 3 parts. */

		/* Make sure the array is clean. */

		for (i = 0; i < numberOfNodes; i++)
		{
			for (j = 0; j < maxNeighbors; j++)
			{
				nodeNeighbors[i][j] = 0;
			}
		}

		/* Part one identifies the corner points of the lattice, each of which has 2 neighbors, */
		/* and identifies them. */

		/* Part One */

		/* Top left hand corner. */
		i = 0;
		numberOfNodeNeighbors[i] = 2;
		nodeNeighbors[i][0]      = 1;
		nodeNeighbors[i][1]      = a;

		/* Top right hand corner. */
		i = a - 1;
		numberOfNodeNeighbors[i] = 2;
	    nodeNeighbors[i][0]      = a - 2;
	    nodeNeighbors[i][1]      = 2 * a - 1;

		/* Bottom left hand corner. */
		i = a * b - a;
		numberOfNodeNeighbors[i] = 2;
	    nodeNeighbors[i][0]      = a * b - 2 * a;
	    nodeNeighbors[i][1]      = a * b - a + 1;

		/* Bottom right hand corner. */
		i = a * b - 1;
		numberOfNodeNeighbors[i] = 2;
	    nodeNeighbors[i][0]      = a * b - a - 1;
	    nodeNeighbors[i][1]      = a * b - 2;

		/* Part Two figures out the lattice points that have 3 neighbors. */
		/* These points reside along the 4 edges of the rectangle and are not corner points. */

		/* Part Two */

		/* Edges in the length [a] direction.  This will not be done if a = 2 since there are no edge nodes then. */

		if (a > 2)
		{
			/* Top edge. */
			for (i = 1; i <= a - 2; i++)
			{
				numberOfNodeNeighbors[i] = 3;
		        nodeNeighbors[i][0]      = i - 1;
		        nodeNeighbors[i][1]      = i + 1;
		        nodeNeighbors[i][2]      = i + a;
			}

			/* Bottom edge. */
		    for (i = a * b - a + 1; i <= a * b - 2; i++)
			{
				numberOfNodeNeighbors[i] = 3;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + 1;
			}
		}

		/* Edges in the length [b] direction. */

		if (b > 2)
		{
			/* Left hand edge. */
		    for (i = a; i <= a * b - 2 * a; i += a)
			{
				numberOfNodeNeighbors[i] = 3;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i + 1;
		        nodeNeighbors[i][2]      = i + a;
			}

			/* Right hand edge. */
		    for (i = 2 * a - 1; i <= a * b - a - 1; i += a)
			{
				numberOfNodeNeighbors[i] = 3;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + a;
			}
		}

		/* Part Three identifies the lattice points that have 4 neighbors. */
	    /* These points all lie within the rectangle. */

		/* Part Three */

		for (i = 0; i < numberOfNodes; i++)
		{
		    if (numberOfNodeNeighbors[i] == 0)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0] = i - a;
		        nodeNeighbors[i][1] = i - 1;
		        nodeNeighbors[i][2] = i + 1;
		        nodeNeighbors[i][3] = i + a;
			}
		}
	}
	else if (latticeDimension == 3)
	{
		/* Get dimensions of lattice. */

		cout << "Enter length of cuboid, a = ";					/* Get the length dimension   a a a a */
		cin  >> a;

		if (a < 2)
			goto error;

		cout << "Enter width of cuboid,  b = ";					/* Get the width dimension   b */
		cin  >> b;												/*                          b  */
																/*                         b   */
																/*                        b    */
		if (b < 2)
			goto error;

		cout << "Enter height of cuboid, c = ";					/* Get the height dimension  c */
		cin  >> c;												/*                           c */
																/*                           c */
																/*                           c */
		if (c < 2)
			goto error;

		numberOfNodes = a * b * c;

		if (numberOfNodes > maxNodes)
			goto error;

		file << "Lattice has 3 dimensions" << endl;
		file << "Length of cuboid, a = " << a << endl;
		file << "Width of cuboid,  b = " << b << endl;
		file << "Height of cuboid, c = " << c << endl;

		/* The first part of the program creates the neighbor array.  This part of the program has 4 parts. */

		/* Make sure the array is clean. */

		for (i = 0; i < numberOfNodes; i++)
		{
			for (j = 0; j < maxNeighbors; j++)
			{
				nodeNeighbors[i][j] = 0;
			}
		}

		/* Part one finds the corner points of the lattice and assigns them the number 3 since they have three neighbors. */
		/* Then, for each corner, it figures out the neighbors. */

		/* Part One */

		/* Upper face, top left hand corner. */
		i = 0;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i + 1;
	    nodeNeighbors[i][1]      = i + a;
	    nodeNeighbors[i][2]      = i + a * b;

		/* Upper face, top right hand corner. */
		i = a - 1;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - 1;
	    nodeNeighbors[i][1]      = i + a;
	    nodeNeighbors[i][2]      = i + a * b;

		/* Upper face, bottom left hand corner. */
		i = a * b - a;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a;
	    nodeNeighbors[i][1]      = i + 1;
	    nodeNeighbors[i][2]      = i + a * b;

		/* Upper face, bottom right hand corner. */
		i = a * b - 1;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a;
	    nodeNeighbors[i][1]      = i - 1;
	    nodeNeighbors[i][2]      = i + a * b;

		/* Lower face, top left hand corner. */
		i = a * b * c - a * b;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a * b;
	    nodeNeighbors[i][1]      = i + 1;
	    nodeNeighbors[i][2]      = i + a;

		/* Lower face, top right hand corner. */
		i = a * b * c - a * b + a - 1;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a * b;
	    nodeNeighbors[i][1]      = i - 1;
	    nodeNeighbors[i][2]      = i + a;

		/* Lower face, bottom left hand corner. */
		i = a * b * c - a;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a * b;
	    nodeNeighbors[i][1]      = i - a;
	    nodeNeighbors[i][2]      = i + 1;

		/* Lower face, bottom right hand corner. */
		i = a * b * c - 1;
		numberOfNodeNeighbors[i] = 3;
	    nodeNeighbors[i][0]      = i - a * b;
	    nodeNeighbors[i][1]      = i - a;
	    nodeNeighbors[i][2]      = i - 1;

		/* Part Two figures out the lattice points that have 4 neighbors. */
		/* These points reside along the edges of the cube. */
		/* Since there are twelve edges, we need twelve loops. */

		/*Part Two*/
		/*Top and bottom face edges in the length [a] direction.  This will not be done if a = 2 since there are no edge nodes then. */

		if (a > 2)
		{
			for (i = 1; i <= a - 2; i++)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - 1;
		        nodeNeighbors[i][1]      = i + 1;
		        nodeNeighbors[i][2]      = i + a;
		        nodeNeighbors[i][3]      = i + a * b;
			}

			for (i = a * b - a + 1; i <= a * b - 2; i++)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + 1;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = a * b * c - a * b + 1; i <= a * b * c - a * b + a - 2; i++)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + 1;
		        nodeNeighbors[i][3]      = i + a;
			}

		    for (i = a * b * c - a + 1; i <= a * b * c - 2; i++)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i - 1;
		        nodeNeighbors[i][3]      = i + 1;
			}
		}

		/* Top and bottom face edges in the length [b] direction. */

		if (b > 2)
		{
		    for (i = a; i <= a * b - 2 * a; i += a)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i + 1;
		        nodeNeighbors[i][2]      = i + a;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = 2 * a - 1; i <= a * b - a - 1; i += a)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + a;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = a * b * c - a * b + a; i <= a * b * c - 2 * a; i += a)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i + 1;
		        nodeNeighbors[i][3]      = i + a;
			}

		    for (i = a * b * c - a * b + 2 * a - 1; i <= a * b * c - a - 1; i += a)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i - 1;
		        nodeNeighbors[i][3]      = i + a;
			}
		}

		/* Lateral edges. */

		if (c > 2)
		{
		    for (i = a * b; i <= a * b * c - 2 * a * b; i += a * b)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i + 1;
		        nodeNeighbors[i][2]      = i + a;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = a * b + a - 1; i <= a * b * c - 2 * a * b + a - 1; i += a * b)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - 1;
		        nodeNeighbors[i][2]      = i + a;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = 2 * a * b - a; i <= a * b * c - a * b - a; i += a * b)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i + 1;
		        nodeNeighbors[i][3]      = i + a * b;
			}

		    for (i = 2 * a * b - 1; i <= a * b * c - a * b - 1; i += a * b)
			{
				numberOfNodeNeighbors[i] = 4;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i - 1;
		        nodeNeighbors[i][3]      = i + a * b;
			}
		}

		/* Part Three */
		/* Each lattice point that is on the inside of a face has 5 neighbors. */

		if (a > 2 && b > 2)
		{
			/* Top face. */

			for (i = 1; i <= b - 2; i++)
			{
		        for (j = 1; j <= a - 2; j++)
				{
					k = i * a + j;
					numberOfNodeNeighbors[k] = 5;
		            nodeNeighbors[k][0]      = k - a;
		            nodeNeighbors[k][1]      = k - 1;
		            nodeNeighbors[k][2]      = k + 1;
		            nodeNeighbors[k][3]      = k + a;
		            nodeNeighbors[k][4]      = k + a * b;
				}
			}

			/* Bottom face. */

			for (i = 1; i <= b - 2; i++)
			{
		        for (j = 1; j <= a - 2; j++)
				{
					k = i * a + j + a * b * (c - 1);
					numberOfNodeNeighbors[k] = 5;
					nodeNeighbors[k][0]      = k - a * b;
		            nodeNeighbors[k][1]      = k - a;
		            nodeNeighbors[k][2]      = k - 1;
		            nodeNeighbors[k][3]      = k + 1;
		            nodeNeighbors[k][4]      = k + a;
				}
			}
		}

		if (a > 2 && c > 2)
		{
			/* Front face. */

			for (i = 1; i <= c - 2; i++)
			{
		        for (j = 1; j <= a - 2; j++)
				{
					k = (i + 1) * a * b + j - a;
					numberOfNodeNeighbors[k] = 5;
		            nodeNeighbors[k][0]      = k - a * b;
		            nodeNeighbors[k][1]      = k - a;
		            nodeNeighbors[k][2]      = k - 1;
		            nodeNeighbors[k][3]      = k + 1;
		            nodeNeighbors[k][4]      = k + a * b;
				}
			}

			/* Back face. */

			for (i = 1; i <= c - 2; i++)
			{
		        for (j = 1; j <= a - 2; j++)
				{
					k = i * a * b + j;
					numberOfNodeNeighbors[k] = 5;
		            nodeNeighbors[k][0]      = k - a * b;
		            nodeNeighbors[k][1]      = k - 1;
		            nodeNeighbors[k][2]      = k + 1;
		            nodeNeighbors[k][3]      = k + a;
		            nodeNeighbors[k][4]      = k + a * b;
				}
			}
		}

		if (b > 2 && c > 2)
		{
			/* Right face. */

			for (i = 1; i <= c - 2; i++)
			{
		        for (j = 1; j <= b - 2; j++)
				{
					k = i * a * b + a * j + a - 1;
					numberOfNodeNeighbors[k] = 5;
		            nodeNeighbors[k][0]      = k - a * b;
		            nodeNeighbors[k][1]      = k - a;
		            nodeNeighbors[k][2]      = k - 1;
		            nodeNeighbors[k][3]      = k + a;
		            nodeNeighbors[k][4]      = k + a * b;
				}
			}

			/* Left face. */

			for (i = 1; i <= c - 2; i++)
			{
		        for (j = 1; j <= b - 2; j++)
				{
					k = i * a * b + j * a;
					numberOfNodeNeighbors[k] = 5;
		            nodeNeighbors[k][0]      = k - a * b;
		            nodeNeighbors[k][1]      = k - a;
		            nodeNeighbors[k][2]      = k + 1;
		            nodeNeighbors[k][3]      = k + a;
		            nodeNeighbors[k][4]      = k + a * b;
				}
			}
		}

		/* Part Four */
		/* Every other node must be inside the cuboid and has 6 neighbors. */

		for (i = 0; i < numberOfNodes; i++)
		{
		    if (numberOfNodeNeighbors[i] == 0)
			{
				numberOfNodeNeighbors[i] = 6;
		        nodeNeighbors[i][0]      = i - a * b;
		        nodeNeighbors[i][1]      = i - a;
		        nodeNeighbors[i][2]      = i - 1;
		        nodeNeighbors[i][3]      = i + 1;
		        nodeNeighbors[i][4]      = i + a;
		        nodeNeighbors[i][5]      = i + a * b;
			}
		}
	}

	file << endl;

	if (latticeDimension == 2)
	{
		/* Generate all tilings of rectangle. */

		latticeArea = (a - 1) * (b - 1);

		if (a <= b)
		{
			maxSquareSize = a;
		}
		else
		{
			maxSquareSize = b;
		}

		/* Allocate space for tiling arrays. */

		tilingList = (bool **) malloc (maxTilings * sizeof(bool *));

		for (i = 0; i < maxTilings; i++)
		{
			tilingList[i] = (bool *) malloc (numberOfNodes * sizeof(bool));
		}

		/* Construct initial tiling with all edges present. */

		for (i = 0; i < numberOfNodes; i++)
		{
			tiling[i] = 1;
		}

		/* Specify square parts. */

		for (i = 0; i <= maxSquareSize; i++)
		{
			parts[i] = 0;
		}

		parts[1] = latticeArea;

		/* Record tiling and parts. */

		numberOfTilings = 1;
		cout << numberOfTilings << endl;

		for (i = 0; i < numberOfNodes; i++)
		{
			tilingList[numberOfTilings-1][i] = tiling[i];
		}

		numberOfPartSets = 1;

		for (i = 0; i <= maxSquareSize; i++)
		{
			partSets[numberOfPartSets - 1][i] = parts[i];
		}

		tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

		/* Output tiling. */

		file << "Tiling " << numberOfTilings << endl;

		for (i = 0; i < numberOfNodes; i++)
		{
			file << tiling[i] << " ";

			if (i % a == a - 1) file << endl;
		}

		file << endl;

		for (s = 2; s < 10; s++)
		{
			if (maxSquareSize > s)
			{
				/* Determine all tilings involving ((a - 1)(b - 1) - s^2) 1x1 squares and 1 sxs square. */

				/* Specify square parts. */

				for (j = 0; j <= maxSquareSize; j++)
				{
					parts[j] = 0;
				}

				parts[1] = latticeArea - s * s;
				parts[s] = 1;

				numberOfPartSets++;

				for (j = 0; j <= maxSquareSize; j++)
				{
					partSets[numberOfPartSets - 1][j] = parts[j];
				}

				for (i = 0; i < numberOfNodes; i++)
				{
					sum  = 0;
					stop = false;

					for (k = 0; k < s - 1; k++)
					{
						for (m = 0; m < s - 1; m++)
						{
							if (numberOfNodeNeighbors[i + k * a + m] == 4)
							{
								sum++;
							}
							else
							{
								stop = true;
								break;
							}
						}

						if (stop == true) break;
					}

					if (sum == (s - 1) * (s - 1))
					{
						/* Start constructing tiling with all edges present. */

						for (j = 0; j < numberOfNodes; j++)
						{
							tiling[j] = 1;
						}

						/* Eliminate connections of internal nodes to make an sxs square. */

						for (k = 0; k < s - 1; k++)
						{
							for (m = 0; m < s - 1; m++)
							{
								tiling[i + k * a + m] = 0;
							}
						}

						/* Record tiling. */

						numberOfTilings++;
						cout << numberOfTilings << endl;

						for (j = 0; j < numberOfNodes; j++)
						{
							tilingList[numberOfTilings-1][j] = tiling[j];
						}

						tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

						/* Output tiling. */

						file << "Tiling " << numberOfTilings << endl;

						for (j = 0; j < numberOfNodes; j++)
						{
							file << tiling[j] << " ";

							if (j % a == a - 1) file << endl;
						}

						file << endl;
					}
				}

				/* Determine all tilings involving 1x1 .. sxs squares. */

				maxSquares = ((a - 1) / s) * ((b - 1) / s);

				firstTiling = 1;
				lastTiling  = numberOfTilings - 1;

				for (numberOfSquares = 1; numberOfSquares <= maxSquares; numberOfSquares++)
				{
					for (tilingIndex = firstTiling; tilingIndex <= lastTiling; tilingIndex++)
					{
						for (i = 0; i < numberOfNodes; i++)
						{
							sum  = 0;
							stop = false;

							for (k = 0; k < s - 1; k++)
							{
								for (m = 0; m < s - 1; m++)
								{
									if (numberOfNodeNeighbors[i + k * a + m] == 4)
									{
										sum++;
									}
									else
									{
										stop = true;
										break;
									}
								}

								if (stop == true) break;
							}

							if (sum == (s - 1) * (s - 1))
							{
								/* Get previously generated tiling and its square parts. */

								for (j = 0; j < numberOfNodes; j++)
								{
									tiling[j] = tilingList[tilingIndex][j];
								}

								for (j = 0; j <= maxSquareSize; j++)
								{
									parts[j] = partSets[tilingPartSet[tilingIndex]][j];
								}

								sum  = 0;
								stop = false;

								for (k = -1; k < s; k++)
								{
									for (m = -1; m < s; m++)
									{
										if (tiling[i + k * a + m] == 1)
										{
											sum++;
										}
										else
										{
											stop = true;
											break;
										}
									}

									if (stop == true) break;
								}

								if (sum == (s + 1) * (s + 1))
								{
									/* Eliminate connections of internal nodes to make an sxs square. */

									for (k = 0; k < s - 1; k++)
									{
										for (m = 0; m < s - 1; m++)
										{
											tiling[i + k * a + m] = 0;
										}
									}

									/* Compare tiling with relevant others to prevent duplication. */

									duplicate = false;

									for (k = lastTiling; k <= numberOfTilings; k++)
									{
										n = 0;

										for (l = 0; l < numberOfNodes; l++)
										{
											if (tiling[l] == tilingList[k][l])
											{
												n = n + 1;
											}
										}

										if (n == numberOfNodes)
										{
											duplicate = true;
										}
									}

									if (duplicate == false)
									{
										/* Record tiling. */

										numberOfTilings++;
										cout << numberOfTilings << endl;

										for (j = 0; j < numberOfNodes; j++)
										{
											tilingList[numberOfTilings-1][j] = tiling[j];
										}

										/* Form square parts. */

										parts[1] = parts[1] - s * s;
										parts[s] = parts[s] + 1;

										for (j = 0; j <= maxSquareSize; j++)
										{
											if (parts[j] != partSets[numberOfPartSets - 1][j])
											{
												break;
											}
										}

										if (j < maxSquareSize)
										{
											numberOfPartSets++;

											for (j = 0; j <= maxSquareSize; j++)
											{
												partSets[numberOfPartSets - 1][j] = parts[j];
											}
										}

										tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

										/* Output tiling. */

										file << "Tiling " << numberOfTilings << endl;

										for (j = 0; j < numberOfNodes; j++)
										{
											file << tiling[j] << " ";

											if (j % a == a - 1) file << endl;
										}

										file << endl;
									}
								}
							}
						}
					}

					/* Update lower and upper tiling indices. */

					firstTiling = lastTiling + 1;
					lastTiling  = numberOfTilings - 1;
				}
			}
		}

		/* Output number of tilings. */

		file << endl;
		file << "Number of tilings = " << numberOfTilings << endl << endl << endl;

		/* Determine the inequivalent tilings. */

		/* Initialise indicator array. */

		for (i = 0; i < numberOfTilings; i++)
		{
			tilingInd[i] = 0;
		}

		/* Calculate rotation and reflection arrays. */

		if (a == b)
		{
			/* Rotate left by 90 degs. */

			for (m = 0; m < b; m++)
			{
				for (n = 0; n < a; n++)
				{
					rotate090[m * a + n] = (a - n - 1) * a + m;
				}
			}
		}

		/* Rotate left by 180 degs. */

		for (m = 0; m < b; m++)
		{
			for (n = 0; n < a; n++)
			{
				rotate180[m * a + n] = (b - m) * a - (n + 1);
			}
		}

		if (a == b)
		{
			/* Rotate left by 270 degs. */

			for (m = 0; m < b; m++)
			{
				for (n = 0; n < a; n++)
				{
					rotate270[m * a + n] = (n + 1) * a - (m + 1);
				}
			}
		}

		/* Reflect horizontally. */

		for (m = 0; m < b; m++)
		{
			for (n = 0; n < a; n++)
			{
				reflectHoriz[m * a + n] = (b - m - 1) * a + n;
			}
		}

		/* Reflect vertically. */

		for (m = 0; m < b; m++)
		{
			for (n = 0; n < a; n++)
			{
				reflectVert[m * a + n] = (m + 1) * a - (n + 1);
			}
		}

		if (a == b)
		{
			/* Reflect horizontally then rotate by 90 degs = Reflect diagonally. */

			for (m = 0; m < b; m++)
			{
				for (n = 0; n < a; n++)
				{
					refHorizRot90[m * a + n] = m + a * n;
				}
			}

			/* Reflect vertically then rotate by 90 degs = Reflect antidiagonally. */

			for (m = 0; m < b; m++)
			{
				for (n = 0; n < a; n++)
				{
					refVertRot90[m * a + n] = (a - n) * a - (m + 1);
				}
			}
		}

		h = 0;

		while (h < numberOfTilings)
		{
			if (tilingInd[h] == 0)
			{
				tilingInd[h]++;

				if (a == b)
				{
					/* Rotate tiling by 90 degs. */

					for (j = 0; j < numberOfNodes; j++)
					{
						tiling[j] = tilingList[h][rotate090[j]];
					}

					/* Find this tiling in the list. */

					for (i = h; i < numberOfTilings; i++)
					{
						for (j = 0; j < numberOfNodes; j++)
						{
							if (tilingList[i][j] != tiling[j]) break;
						}

						if ((j == numberOfNodes) &&
							(tilingInd[i] == 0))
						{
							tilingInd[i] = -1;
							tilingInd[h]++;
							break;
						}
					}
				}

				/* Rotate tiling by 180 degs. */

				for (j = 0; j < numberOfNodes; j++)
				{
					tiling[j] = tilingList[h][rotate180[j]];
				}

				/* Find this tiling in the list. */

				for (i = h; i < numberOfTilings; i++)
				{
					for (j = 0; j < numberOfNodes; j++)
					{
						if (tilingList[i][j] != tiling[j]) break;
					}

					if ((j == numberOfNodes) &&
						(tilingInd[i] == 0))
					{
						tilingInd[i] = -1;
						tilingInd[h]++;
						break;
					}
				}

				if (a == b)
				{
					/* Rotate tiling by 270 degs. */

					for (j = 0; j < numberOfNodes; j++)
					{
						tiling[j] = tilingList[h][rotate270[j]];
					}

					/* Find this tiling in the list. */

					for (i = h; i < numberOfTilings; i++)
					{
						for (j = 0; j < numberOfNodes; j++)
						{
							if (tilingList[i][j] != tiling[j]) break;
						}

						if ((j == numberOfNodes) &&
							(tilingInd[i] == 0))
						{
							tilingInd[i] = -1;
							tilingInd[h]++;
							break;
						}
					}
				}

				/* Reflect tiling horizontally. */

				for (j = 0; j < numberOfNodes; j++)
				{
					tiling[j] = tilingList[h][reflectHoriz[j]];
				}

				/* Find this tiling in the list. */

				for (i = h; i < numberOfTilings; i++)
				{
					for (j = 0; j < numberOfNodes; j++)
					{
						if (tilingList[i][j] != tiling[j]) break;
					}

					if ((j == numberOfNodes) &&
						(tilingInd[i] == 0))
					{
						tilingInd[i] = -1;
						tilingInd[h]++;
						break;
					}
				}

				/* Reflect tiling vertically. */

				for (j = 0; j < numberOfNodes; j++)
				{
					tiling[j] = tilingList[h][reflectVert[j]];
				}

				/* Find this tiling in the list. */

				for (i = h; i < numberOfTilings; i++)
				{
					for (j = 0; j < numberOfNodes; j++)
					{
						if (tilingList[i][j] != tiling[j]) break;
					}

					if ((j == numberOfNodes) &&
						(tilingInd[i] == 0))
					{
						tilingInd[i] = -1;
						tilingInd[h]++;
						break;
					}
				}

				if (a == b)
				{
					/* Reflect tiling horizontally then rotate by 90 degs = Reflect diagonally. */

					for (j = 0; j < numberOfNodes; j++)
					{
						tiling[j] = tilingList[h][refHorizRot90[j]];
					}

					/* Find this tiling in the list. */

					for (i = h; i < numberOfTilings; i++)
					{
						for (j = 0; j < numberOfNodes; j++)
						{
							if (tilingList[i][j] != tiling[j]) break;
						}

						if ((j == numberOfNodes) &&
							(tilingInd[i] == 0))
						{
							tilingInd[i] = -1;
							tilingInd[h]++;
							break;
						}
					}

					/* Reflect tiling vertically then rotate by 90 degs = Reflect antidiagonally. */

					for (j = 0; j < numberOfNodes; j++)
					{
						tiling[j] = tilingList[h][refVertRot90[j]];
					}

					/* Find this tiling in the list. */

					for (i = h; i < numberOfTilings; i++)
					{
						for (j = 0; j < numberOfNodes; j++)
						{
							if (tilingList[i][j] != tiling[j]) break;
						}

						if ((j == numberOfNodes) &&
							(tilingInd[i] == 0))
						{
							tilingInd[i] = -1;
							tilingInd[h]++;
							break;
						}
					}
				}
			}

			h++;
		}

		/* Output tilings. */

		for (i = 0; i < numberOfTilings; i++)
		{
			for (j = 0; j <  numberOfNodes; j++)
			{
				file << tilingList[i][j];
			}

			file << endl;
		}

		file << endl << endl;

		/* Output counts of tilings with each degree of symmetry. */

		for (i = 0; i < 10; i++)
		{
			count[i] = 0;
		}

		for (h = 0; h < numberOfTilings; h++)
		{
			count[tilingInd[h] + 1]++;
		}

		file << "Numbers of tilings with different numbers of images" << endl << endl;

		for (i = -1; i < 9; i++)
		{
			file << i << "    ";
		}

		file << endl;

		for (i = 0; i < 10; i++)
		{
			file << count[i] << "    ";
		}

		file << endl << endl << endl;

		/* Output number of inequivalent tilings. */

		n = 0;

		for (i = 2; i < 10; i++)
		{
			n = n + count[i];
		}

		file << "Number of inequivalent tilings = " << n << endl << endl << endl;

		/* Form histograms of the number of nodes in each unrestricted tiling, inequivalent tiling and part set not connected to any neighbors. */

		for (i = 0; i <= numberOfNodes; i++)
		{
			numberOfZeros[i] = 0;					/* Histogram of number of 0s in each unrestricted tiling. */
			numberOfZerosInequivalent[i] = 0;		/* Histogram of number of 0s in each inequivalent tiling. */
			numberOfZerosPartSet[i] = 0;			/* Histogram of number of 0s in each part set. */
		}

		for (i = 0; i < numberOfTilings; i++)
		{
			zeroCount = 0;

			for (j = 0; j <  numberOfNodes; j++)
			{
				if (tilingList[i][j] == 0)
				{
					zeroCount++;
				}
			}

			numberOfZeros[zeroCount]++;

			if (tilingInd[i] >= 0)
			{
				numberOfZerosInequivalent[zeroCount]++;
			}
		}

		for (i = 0; i < numberOfPartSets; i++)
		{
			zeroCount = 0;

			for (j = 1; j <= a; j++)
			{
				zeroCount = zeroCount + partSets[i][j] * (j - 1) * (j - 1);
			}

			numberOfZerosPartSet[zeroCount]++;
		}

		file << "Numbers of part sets (P), inequivalent (I) and unrestricted (U) tilings containing different numbers of zeros (Z)" << endl << endl;
		file << "  Z     P     I     U" << endl << endl;

		for (i = 0; i <= (a - 2) * (a - 2); i++)
		{
			file << "  " << i <<  "     " << numberOfZerosPartSet[i] << "     " << numberOfZerosInequivalent[i] << "     " << numberOfZeros[i] << endl;
		}

		file << endl << endl;

		/* Generate square parts stats. */

		for (i = 0; i < maxPartSets; i++)
		{
			partSetCount[i] = 0;
			partSetCountInequivalent[i] = 0;
		}

		for (i = 0; i < numberOfTilings; i++)
		{
			partSetCount[tilingPartSet[i]]++;

			if (tilingInd[i] >= 0)
			{
				partSetCountInequivalent[tilingPartSet[i]]++;
			}
		}

		file << "Numbers of inequivalent (I) and unrestricted (U) tilings with each different set of square parts" << endl << endl;
		file << "Square side   ";

		for (i = 1; i <= a; i++)
		{
			file << i << " ";
		}

		file << endl;
		file << "    Counts" << endl;
		file << "    I    U" << endl;

		for (i = 0; i < numberOfPartSets; i++)
		{
			file << "    " << partSetCountInequivalent[i] << "    " << partSetCount[i] << "    ";

			for (j = 1; j <= a; j++)
			{
				file << partSets[i][j] << " ";
			}

			file << endl;
		}

		file << endl << endl;

		/* Output number of part sets. */

		file << "Number of distinct part sets = " << numberOfPartSets << endl << endl;

		/* Release allocated memory. */

		for (i = 0; i < maxTilings; i++)
		{
			free (tilingList[i]);
		}

		free (tilingList);
	}
	else if (latticeDimension == 3)
	{
		/* Generate all tilings of rectangular cuboid. */

		latticeVolume = (a - 1) * (b - 1) * (c - 1);

		if ((a <= b) && (a <= c))
		{
			maxCubeSize = a;
		}

		if ((b <= a) && (b <= c))
		{
			maxCubeSize = b;
		}

		if ((c <= a) && (c <= b))
		{
			maxCubeSize = c;
		}

		/* Allocate space for tiling arrays. */

		tilingList = (bool **) malloc (maxTilings * sizeof(bool *));

		for (i = 0; i < maxTilings; i++)
		{
			tilingList[i] = (bool *) malloc (numberOfNodes * sizeof(bool));
		}

		/* Construct initial tiling with all edges present. */

		for (i = 0; i < numberOfNodes; i++)
		{
			tiling[i] = 1;
		}

		/* Specify cubic parts. */

		for (i = 0; i <= maxCubeSize; i++)
		{
			parts[i] = 0;
		}

		parts[1] = latticeVolume;

		/* Record tiling and parts. */

		numberOfTilings = 1;
		cout << numberOfTilings << endl;

		for (i = 0; i < numberOfNodes; i++)
		{
			tilingList[numberOfTilings-1][i] = tiling[i];
		}

		numberOfPartSets = 1;

		for (i = 0; i <= maxCubeSize; i++)
		{
			partSets[numberOfPartSets - 1][i] = parts[i];
		}

		tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

		/* Output 3D tiling using c 2D tilings of size a x b. */
/*
		file << "Tiling " << numberOfTilings << endl << endl;

		m = 0;

		for (k = 0; k < numberOfNodes; k++)
		{
			file << tiling[k] << " ";

			if (k % a == a - 1)
			{
				file << endl;

				if (m % b == b - 1)
				{
					file << endl;
					m = 0;
				}
				else
				{
					m++;
				}
			}
		}

		file << endl;
*/
		for (s = 2; s < 10; s++)
		{
			if (maxCubeSize > s)
			{
				/* Determine all tilings involving ((a - 1)(b - 1)(c - 1) - s^3) 1x1x1 cubes and 1 sxsxs cube. */

				/* Specify cubic parts. */

				for (j = 0; j <= maxCubeSize; j++)
				{
					parts[j] = 0;
				}

				parts[1] = latticeVolume - s * s * s;
				parts[s] = 1;

				numberOfPartSets++;

				for (j = 0; j <= maxCubeSize; j++)
				{
					partSets[numberOfPartSets - 1][j] = parts[j];
				}

				for (i = 0; i < numberOfNodes; i++)
				{
					sum  = 0;
					stop = false;

					for (j = 0; j < s - 1; j++)
					{
						for (k = 0; k < s - 1; k++)
						{
							for (m = 0; m < s - 1; m++)
							{
								if (numberOfNodeNeighbors[i + j*a*b + k*a + m] == 6)
								{
									sum++;
								}
								else
								{
									stop = true;
									break;
								}
							}

							if (stop == true) break;
						}

						if (stop == true) break;
					}

					if (sum == (s - 1) * (s - 1) * (s - 1))
					{
						/* Start constructing tiling with all edges present. */

						for (j = 0; j < numberOfNodes; j++)
						{
							tiling[j] = 1;
						}

						/* Eliminate connections of internal nodes to make an sxsxs cube. */

						for (j = 0; j < s - 1; j++)
						{
							for (k = 0; k < s - 1; k++)
							{
								for (m = 0; m < s - 1; m++)
								{
									tiling[i + j*a*b + k*a + m] = 0;
								}
							}
						}

						/* Record tiling. */

						numberOfTilings++;
						cout << numberOfTilings << endl;

						for (j = 0; j < numberOfNodes; j++)
						{
							tilingList[numberOfTilings-1][j] = tiling[j];
						}

						tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

						/* Output 3D tiling using c 2D tilings of size a X b. */
/*
						file << "Tiling " << numberOfTilings << endl << endl;

						m = 0;

						for (k = 0; k < numberOfNodes; k++)
						{
							file << tiling[k] << " ";

							if (k % a == a - 1)
							{
								file << endl;

								if (m % b == b - 1)
								{
									file << endl;
									m = 0;
								}
								else
								{
									m++;
								}
							}
						}

						file << endl;
*/
					}
				}

				/* Determine all tilings involving 1x1x1 .. sxsxs cubes. */

				maxCubes = ((a - 1) / s) * ((b - 1) / s) * ((c - 1) / s);

				firstTiling = 1;
				lastTiling  = numberOfTilings - 1;

				for (numberOfCubes = 1; numberOfCubes <= maxCubes; numberOfCubes++)
				{
					for (tilingIndex = firstTiling; tilingIndex <= lastTiling; tilingIndex++)
					{
						for (i = 0; i < numberOfNodes; i++)
						{
							sum  = 0;
							stop = false;

							for (j = 0; j < s - 1; j++)
							{
								for (k = 0; k < s - 1; k++)
								{
									for (m = 0; m < s - 1; m++)
									{
										if (numberOfNodeNeighbors[i + j*a*b + k*a + m] == 6)
										{
											sum++;
										}
										else
										{
											stop = true;
											break;
										}
									}

									if (stop == true) break;
								}

								if (stop == true) break;
							}

							if (sum == (s - 1) * (s - 1) * (s - 1))
							{
								/* Get previously generated tiling and its cubic parts. */

								for (j = 0; j < numberOfNodes; j++)
								{
									tiling[j] = tilingList[tilingIndex][j];
								}

								for (j = 0; j <= maxCubeSize; j++)
								{
									parts[j] = partSets[tilingPartSet[tilingIndex]][j];
								}

								sum  = 0;
								stop = false;

								for (j = -1; j < s; j++)
								{
									for (k = -1; k < s; k++)
									{
										for (m = -1; m < s; m++)
										{
											if (tiling[i + j*a*b + k*a + m] == 1)
											{
												sum++;
											}
											else
											{
												stop = true;
												break;
											}
										}

										if (stop == true) break;
									}

									if (stop == true) break;
								}

								if (sum == (s + 1) * (s + 1) * (s + 1))
								{
									/* Eliminate connections of internal nodes to make an sxsxs cube. */

									for (j = 0; j < s - 1; j++)
									{
										for (k = 0; k < s - 1; k++)
										{
											for (m = 0; m < s - 1; m++)
											{
												tiling[i + j*a*b + k*a + m] = 0;
											}
										}
									}

									/* Compare tiling with relevant others to prevent duplication. */

									duplicate = false;

									for (k = lastTiling; k <= numberOfTilings; k++)
									{
										n = 0;

										for (l = 0; l < numberOfNodes; l++)
										{
											if (tiling[l] == tilingList[k][l])
											{
												n = n + 1;
											}
										}

										if (n == numberOfNodes)
										{
											duplicate = true;
										}
									}

									if (duplicate == false)
									{
										/* Record tiling. */

										numberOfTilings++;
										cout << numberOfTilings << endl;

										for (j = 0; j < numberOfNodes; j++)
										{
											tilingList[numberOfTilings-1][j] = tiling[j];
										}

										/* Form cubic parts. */

										parts[1] = parts[1] - s * s * s;
										parts[s] = parts[s] + 1;

										for (j = 0; j <= maxCubeSize; j++)
										{
											if (parts[j] != partSets[numberOfPartSets - 1][j])
											{
												break;
											}
										}

										if (j < maxCubeSize)
										{
											numberOfPartSets++;

											for (j = 0; j <= maxCubeSize; j++)
											{
												partSets[numberOfPartSets - 1][j] = parts[j];
											}
										}

										tilingPartSet[numberOfTilings-1] = numberOfPartSets - 1;

										/* Output 3D tiling using c 2D tilings of size a X b. */
/*
										file << "Tiling " << numberOfTilings << endl << endl;

										m = 0;

										for (k = 0; k < numberOfNodes; k++)
										{
											file << tiling[k] << " ";

											if (k % a == a - 1)
											{
												file << endl;

												if (m % b == b - 1)
												{
													file << endl;
													m = 0;
												}
												else
												{
													m++;
												}
											}
										}

										file << endl;
*/
									}
								}
							}
						}
					}

					/* Update lower and upper tiling indices. */

					firstTiling = lastTiling + 1;
					lastTiling  = numberOfTilings - 1;
				}
			}
		}

		/* Output number of tilings. */

		file << endl;
		file << "Number of tilings = " << numberOfTilings << endl << endl << endl;

		/* Output tilings. */
/*
		for (i = 0; i < numberOfTilings; i++)
		{
			for (j = 0; j <  numberOfNodes; j++)
			{
				file << tilingList[i][j];
			}

			file << endl;
		}

		file << endl << endl;
*/
		/* Form histograms of the number of nodes in each unrestricted tiling and part set not connected to any neighbors. */

		for (i = 0; i <= numberOfNodes; i++)
		{
			numberOfZeros[i] = 0;					/* Histogram of number of 0s in each unrestricted tiling. */
			numberOfZerosPartSet[i] = 0;			/* Histogram of number of 0s in each part set. */
		}

		for (i = 0; i < numberOfTilings; i++)
		{
			zeroCount = 0;

			for (j = 0; j <  numberOfNodes; j++)
			{
				if (tilingList[i][j] == 0)
				{
					zeroCount++;
				}
			}

			numberOfZeros[zeroCount]++;
		}

		for (i = 0; i < numberOfPartSets; i++)
		{
			zeroCount = 0;

			for (j = 1; j <= a; j++)
			{
				zeroCount = zeroCount + partSets[i][j] * (j - 1) * (j - 1) * (j - 1);
			}

			numberOfZerosPartSet[zeroCount]++;
		}

		file << "Numbers of part sets (P) and unrestricted (U) tilings containing different numbers of zeros (Z)" << endl << endl;
		file << "  Z     P     U" << endl << endl;

		for (i = 0; i <= (a - 2) * (a - 2) * (a - 2); i++)
		{
			file << "  " << i <<  "     " << numberOfZerosPartSet[i] << "     " << numberOfZeros[i] << endl;
		}

		file << endl << endl;

		/* Generate cubic parts stats. */

		for (i = 0; i < maxPartSets; i++)
		{
			partSetCount[i] = 0;
		}

		for (i = 0; i < numberOfTilings; i++)
		{
			partSetCount[tilingPartSet[i]]++;
		}

		file << "Number of unrestricted (U) tilings with each different set of cubic parts" << endl;
		file << "U    Cube side" << endl;
		file << "     ";

		for (i = 1; i <= a; i++)
		{
			file << i << " ";
		}

		file << endl;

		for (i = 0; i < numberOfPartSets; i++)
		{
			file << partSetCount[i] << "    ";

			for (j = 1; j <= a; j++)
			{
				file << partSets[i][j] << " ";
			}

			file << endl;
		}

		file << endl << endl;

		/* Output number of part sets. */

		file << "Number of distinct part sets = " << numberOfPartSets << endl << endl;

		/* Release allocated memory. */

		for (i = 0; i < maxTilings; i++)
		{
			free (tilingList[i]);
		}

		free (tilingList);
	}

	cout << endl << "The program is done" << endl;
	file << endl << "The program is done" << endl;
	goto done;

error:

	cout << endl << "Input error: Program terminated" << endl;
	file << endl << "Input error: Program terminated" << endl;

done:

	file.close();

	return 0;
}