/*
#
#  Copyright (c) 2008    David Albert Bagley, bagleyd@tux.org
#
#                   All Rights Reserved
#
#  Permission to use, copy, modify, and distribute this software and
#  its documentation for any purpose and without fee is hereby granted,
#  provided that the above copyright notice appear in all copies and
#  that both that copyright notice and this permission notice appear in
#  supporting documentation, and that the name of the author not be
#  used in advertising or publicity pertaining to distribution of the
#  software without specific, written prior permission.
#
#  This program is distributed in the hope that it will be "usable",
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#
#  NOTE: you may THRASH your hard drive if you run out of RAM...
#  IF YOU HAVE LOTS OF DISK ACTIVITY AFTER EXECUTING THIS, TERMINATE IT!
cc -O algorithme.c -lm -o me
nohup nice me 9 2>&1 > me.out &
*/

/*
  System Configuration: Intel Core 2 Duo CPU
  (T7300)
  System clock frequency: 2 GHz
  Memory size: 2 GB
  Using 1 of 2 CPUs
*/

/* Coded after reading "Analysis of the Algorithme 6 Puzzle               */
/* and its Generalisations" by Dick Hess, Cubism For Fun 76               */
/* July 2008 pp 8-13                                                      */

/* This program calculates the algorithme puzzle problem for              */
/* constant height tiles (in example 'N' is 3).                           */

/* Input:                                                                 */
/* # of purple or yellow tiles                                            */
/* # of nodes you expect it to take to search for a solution.             */
/*   (See minSizes table below.)                                          */

/* Output:                                                                */
/* The optimum route to exchange the tiles, but must obey rule            */
/* no tile can go through another tile except top tiles on 1 or 3 can     */
/* pass "around" stack 2 if it is full.  Longer tiles can not be          */
/* on shorter tiles.  Tiles can only move from one stack to               */
/* another along the top.  The width between stacks is one tile width.    */
/* If it takes 2^22 nodes to search for a solution then it will print     */
/* the node and its current move (algorithme order 6 and above).  As it   */
/* stands now a full solution will not be printed for order 7 or more.    */

/* int minMoves[9] = {3, 16, 39, 89, 190, 395, 808, 1637, 3298};          */
/* int minStartMoves[9] = {0, 2, 7, 19, 44, 96, 201, 413, 838};           */
/* int minChurnMoves[9] = {3, 12, 24, 48, 192, 384, 768, 1536};           */
/* int minFinishMoves[9] = {0, 2, 8, 22, 50, 107, 223, 456, 924};         */
/* Say n=tiles                                                            */
/* minStartMoves(n+1) = minStartMoves(n) + floor((3n+1)/2) for n > 3      */
/* minChurnMoves(n+1) = 2*minChurnMoves(n) for n > 2                      */
/* minFinishMoves(n+1) = minFinishMoves(n) + floor(3n/2) for n > 4        */
/* minMoves(n) = minStartMoves(n) + minChurnMoves(n) + minFinishMoves(n)  */

/* minTime[9] = {0., 0., 0., 0., 0., 1sec, 18sec, 7min, 65min};           */

/* To get this to program to run for order 7, the memory had to be cut    */
/* down but this throws away the path to the answer and just leaves you   */
/* with the number of moves.                                              */

/* Farthest moves do not lead to anything very interesting.  It always    */
/* includes the starting point.  The furthest positions have some tiles   */
/* on the opposite side (from top to second from bottom could be          */
/* swapped).  In level 3, a top tile could be in center stack.            */

/*
4x savings from mail by Nick Baxter

The paper also describes a 4x savings, taking advantage of
the symmetry of the puzzle, and by attacking the solution
from both ends. This makes the "solution" check
a bit more complicated (for each new node x, you must
compute and check for x', x* and x*' as well as just x). But
this saves a lot of space.

Depending on the problem, the number of states being analyzed by
the Dykstra algorithm could grow rapidly as the number of moves
increases. The number of states could also remain constant or go
down, but this is less common. I think you characterized the
growth for this puzzle as increasing, but modest. Let's say the average
growth rate for new nodes each turn is p (p>1), and the minimal
solution is 2n. That gives us (p^(2n)-1)/(p-1) total nodes if we just
go from start to finish. If on the other hand we solve simultaneously
from either end, then the number of nodes visited is just 2(p^n-1)/(p-1).
The latter approach saves (p^n-1)^2 nodes! I think I'm +/-1 off
on the exponents, but you get the idea.

Remember that p is growth rate of NEW nodes. If p is just +1%,
then this is 10^19 nodes for level 7! So solving from both ends
can save a lot of search space.

Panex has both color and mirror symmetries. Solving from
the goal is really the same search, but with the puzzle mirrored
OR the colors swapped. Their notation was x* is a color swap,
and x' is a mirror image. If x is reachable from the start in n moves,
then x' and x* are reachable from the finish in n moves. So every time you
add a new position to the list, you check to see if x' or x* are already
there. If so, then you've found a solution for the entire puzzle that
is potentially minimal.

So, as you are building a new n-list for n+1 moves, and you find
a matching x' or x*, it will corresponding to an n or n+1 move path
from the end, giving a 2n+1 or 2n+2 solution. Because of this, I
believe you have to continue building the complete n-list because a
2n+2 solution could still be improved by a 2n+1 solution. It's
conceivable that a 1 move savings could be overlooked if this
isn't done carefully.

There is also the case of x'* (x*' is the same). If x is n moves from
the start, then so is x'*. Thus if x'* is already on a list, you do not
need to add x. This gives additional space savings.

The cost is that you have to compute x', x*, and x'* each time
you need to check if a node exists yet. The assumption is that
the space savings is worth this work.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define FALSE		0
#define TRUE		1
#define LEFT		0
#define MIDDLE		1
#define RIGHT		2
#define MAXSTACKS	3  /* maximum number of stacks */
#define KINDS		2  /* number of different colors of tiles */
#define POSSIBILITIES	6  /* maximum number of possible moves */
#define CHAINS		3  /* for deleting nodes from hash table */
#define ARRAYTWEAKS	10 /* if <= tiles it will consult array for tweaks */
#define FREQBITS	22 /* sets frequency of what node & move its at now */
#define NOFULLSOLUTION	7  /* lower this if you do not have the memory */
/*
Note:
 For full solution size can not be set lower than the maximum number
 of nodes it takes to make "CHAINS" moves anywhere in the search...
 Also if the size is greater than minSizes, a solution path can not be
 generated.
*/

int debug = 0;
int mode = 0;
int arch = 8;
int start = 0, finish = 3;
#define MAX(i, j) (((i) > (j)) ? (i) : (j))
#define MIN(i, j) (((i) < (j)) ? (i) : (j))
#define CEIL(i) ((int) (i + 0.999999)) /* rough ceiling */
/* min bits needed to store a move */
#define MINBITS() CEIL(log2(MAXSTACKS * (tiles + blanks)))
/* max tiles to store in a unit of memory */
#define MAXNUMINUNIT() MIN(((int) (arch / MINBITS())), KINDS * tiles)
/* min units to store all tiles */
#define MINUNITS() CEIL((double) (KINDS * tiles) / MAXNUMINUNIT())
#define MEET_AT_MIDDLE 0
#define STRAIGHT 1
#define FURTHER 2
#define MODES 3

/*#define ARCH128*/
/*#define ARCH64*/

#ifdef ARCH128
#include <bitset>
using namespace std;
const int len = 128;
/* does not work past algorithme 4 due to to_ulong() conversion, also slow */
#define TO_LONG(x) (x).to_ulong()
#define LONG bitset<len>
#else
#define TO_LONG(x) (x)
#define LONG unsigned long long
#endif

#define LOOKUP(x) ((units>1)?(TO_LONG(*x+((*(x+1))*253))%hashSize):\
	(TO_LONG(*x)%hashSize))

static int MEMSTACKS = (sizeof (char)) * MAXSTACKS;
static int EMPTY;

/* already assumes l is not EMPTY */
#define tileNumber(l) ((l < tiles) ? l : l - tiles)
#define tileNum(l) ((l < tiles) ? ((l) << 1) : (((l - tiles) << 1) + 1))
#define numTile(l) ((l >> 1) + ((l & 1) ? tiles : 0))
#define tileChar(l) ((l == EMPTY) ? '-' : ((l < tiles) ? l + 'A' :\
	l + 'a' - tiles))
#define SETBIT(t,p) (((LONG) (p)) << (bits * (t % numInUnit)))
#define setStart(t,p) startTiles[t / numInUnit] |= SETBIT(t,p)
#define setFinish(t,p) finishTiles[t / numInUnit] |= SETBIT(t,p)
#define setStartMid(t,p) startMidTiles[t / numInUnit] |= SETBIT(t,p)
#define setFinishMid(t,p) finishMidTiles[t / numInUnit] |= SETBIT(t,p)

static int hashSizes[ARRAYTWEAKS] = {31, 253, 1021, 65533, 131071,
  262141, 1048573, 8388607, 259000001, 259000001};

/* Maximum number of nodes for a move.  Seems to grow by a factor of 7.77  */
/* This is only for finding the number of moves... a full path is not      */
/* generated.                                                              */
static int moveSizes[MODES][ARRAYTWEAKS] = {
  {1, 6, 32, 214, 1556, 9216, 66004, 435992, 3097980, /*?*/ 83000000},
  {1, 6, 64, 436, 3052, 18044, 124048, 849216, /*?*/ 6200000, /*?*/ 83000000},
  {1, 6, 64, 436, 3052, 18044, 124048, 849216, /*?*/ 6200000, /*?*/ 83000000}};
/* The number of nodes seems to go up roughly by a factor of 18.8 - 20.0.  */
static long long minSizes[MODES][ARRAYTWEAKS] = {
  {4, 23, 318, 3292, 36736, 426128, 5002564, 59392976, 709589548LL, /*?*/ 50290380020LL},
  {4, 43, 783, 10160, 123899, 1491767, 17913061, 214984397LL, /*?*/ 69105921539LL, /*?*/ 1306101917087LL},
  {4, 44, 784, 10160, 123904, 1491776, 17913088, 214984448LL, /*?*/ 69798594048LL, /*?*/ 1312213568102LL}};

static LONG *nodeLog = NULL;
static char *backPoss = NULL;
static int *move = NULL, *backPath = NULL;
static int newNode = 0, currNode = 0;

/* hashless, debug level 2 check */
/* slow but good to check the hashes are working */
static int newBottom = 0;
static int unitMem;

static LONG *startTiles = NULL, *finishTiles = NULL;
static LONG *startMidTiles = NULL, *finishMidTiles = NULL;

typedef struct _hashtype {
	int node;
	struct _hashtype *next;
} hashtype;

static hashtype **hash;
static hashtype *chain[CHAINS];
static int hashNewIndex;
static long long c /* , nc, fc, fnc */;

static char *tileOfPosition = NULL; /* easier to see */
static char *positionOfTile = NULL; /* takes up less memory */
static LONG *shrunkPositionOfTile = NULL; /* takes even less memory */
static LONG *colorPositionOfTile = NULL; /* color swap of previous */
static LONG *mirrorColorPositionOfTile = NULL; /* mirror of previous */
static LONG *mirrorPositionOfTile = NULL; /* mirror of sPOT */

static int tiles; /* # of tiles of the same kind */
static int blanks; /* tiles - 2 */
static int height; /* tiles + blanks */
static int areaPositions; /* MAXSTACKS * (tiles + blanks) */
static int areaTiles; /* KINDS * tiles */
static int size, hashSize;
static long long inputSize;
static int bits, units, numInUnit;
static LONG bitmask;
static int freqmask;
static int loop;

/* chosen by efficiency */
/*static int fromPoss[POSSIBILITIES] = {0, 0, 1, 1, 2, 2};*/
static int toPoss[POSSIBILITIES] = {1, 2, 2, 0, 0, 1};
#define fromPoss(n) (n >> 1)
/*#define toPoss(n) (((n + MAXSTACKS) >> 1) % MAXSTACKS) */

static void
setColorSwap() {
	LONG *pcPOT = colorPositionOfTile;
	char *ppOT = positionOfTile;
	int j, shift, n = 0;

	for (;;) {
		shift = 0;
		*pcPOT = 0;
		for (j = 0; j < numInUnit; j++) {
			if (n & 1)
				ppOT--;
			else
				ppOT++;
			*pcPOT |= (((LONG) *ppOT) << shift);
			if (n & 1) {
				ppOT++;
				ppOT++;
			}
			shift += bits;
			n++;
			if (n >= areaTiles) {
				return;
			}
		}
		pcPOT++;
	}
}

static void
setMirrorImage(LONG *psPOT, LONG *pmPOT) {
	int j, shift, n = 0;
	unsigned long value;

	for (;;) {
		shift = 0;
		*pmPOT = 0;
		for (j = 0; j < numInUnit; j++) {
			value = TO_LONG((*psPOT >> shift) & bitmask);
			value = (value / MAXSTACKS) * MAXSTACKS + MAXSTACKS -
				1 - (value % MAXSTACKS);
			*pmPOT |= (((LONG) value) << shift);
			shift += bits;
			n++;
			if (n >= areaTiles) {
				return;
			}
		}
		psPOT++;
		pmPOT++;
	}
}

static void freeHash(int index) {
	hashtype **oldChain = &(chain[index]);
	hashtype *ctemp = *oldChain;
	hashtype **oldHashPt;
	hashtype *htemp, *temp;

	while (ctemp) {
		oldHashPt = &(hash[LOOKUP(&nodeLog[ctemp->node * units])]);
		htemp = *oldHashPt;

		if ((*oldHashPt)->node == ctemp->node) {
			*oldHashPt = htemp->next;
			(void) free((void *) htemp);
			/* fnc++; */
		} else {
			while (htemp->next->node != ctemp->node)
				htemp = htemp->next;
			temp = htemp->next;
			htemp->next = temp->next;
			(void) free((void *) temp);
			/* fc++; */
		}
		*oldChain = (*oldChain)->next;
		(void) free((void *) ctemp);
		ctemp = *oldChain;
	}
}

#ifdef DEBUG
static void printHash(){
	int i;

	for (i = 0; i < hashSize; i++) {
		hashtype *temp = hash[i];

		while (temp != NULL) {
			(void) printf("%d: %d %ld\n", i, temp->node,
				nodeLog[temp->node * units]);
			temp = temp->next;
		}
	}
}
#endif

static int inHash() {
	LONG *psPOT = shrunkPositionOfTile;
	LONG *tempsPOT, *pnodeLog;
	hashtype *temp;
	int i;

	temp = hash[LOOKUP(psPOT)];
	while (temp != NULL) {
		tempsPOT = psPOT;
		pnodeLog = &(nodeLog[temp->node * units]);
		for (i = 0; i < units; i++) {
			if (*pnodeLog != *tempsPOT)
				goto HASHCONTINUE;
			pnodeLog++;
			tempsPOT++;
		}
		return TRUE;
HASHCONTINUE:	temp = temp->next;
	}
	setColorSwap();
	setMirrorImage(colorPositionOfTile, mirrorColorPositionOfTile);
	psPOT = mirrorColorPositionOfTile;
	temp = hash[LOOKUP(psPOT)];
	while (temp != NULL) {
		tempsPOT = psPOT;
		pnodeLog = &(nodeLog[temp->node * units]);
		for (i = 0; i < units; i++) {
			if (*pnodeLog != *tempsPOT)
				goto CSMICONTINUE;
			pnodeLog++;
			tempsPOT++;
		}
		return TRUE;
CSMICONTINUE:	temp = temp->next;
	}
	return FALSE;
}

/* MEET_AT_MIDDLE */
#define COLORSWAP
#define MIRRORIMAGE
/* The same as previous except that it returns node with shortest path */
static int checkInHash(
	LONG *pcPOT,
	LONG *pmPOT) {
	LONG *tempsPOT, *pnodeLog;
	hashtype *temp;
	int i, flag = -1;

#ifdef COLORSWAP
	temp = hash[LOOKUP(pcPOT)];
	while (temp != NULL) {
		tempsPOT = pcPOT;
		pnodeLog = &(nodeLog[temp->node * units]);
		for (i = 0; i < units; i++) {
			if (*pnodeLog != *tempsPOT)
				goto CSCONTINUE;
			pnodeLog++;
			tempsPOT++;
		}
		if (debug)
			(void) printf("colorswap %d %d\n",
				move[temp->node], move[newNode]);
		if ((flag == -1) || (move[flag] > move[temp->node]))
			flag = temp->node;
CSCONTINUE:	temp = temp->next;
	}
#endif
#ifdef MIRRORIMAGE
	temp = hash[LOOKUP(pmPOT)];
	while (temp != NULL) {
		tempsPOT = pmPOT;
		pnodeLog = &(nodeLog[temp->node * units]);
		for (i = 0; i < units; i++) {
			if (*pnodeLog != *tempsPOT)
				goto MICONTINUE;
			pnodeLog++;
			tempsPOT++;
		}
		if (debug)
			(void) printf("mirrorimage %d %d\n",
				move[temp->node], move[newNode]);
		if ((flag == -1) || (move[flag] > move[temp->node]))
			flag = temp->node;
MICONTINUE:	temp = temp->next;
	}
#endif
	return flag;
}

static void memoryError(const char * string) {
	(void) fprintf(stderr, "Not enough memory for %s.\n", string);
	(void) fprintf(stderr, "Try a smaller number of nodes with -n.\n");
	(void) fprintf(stderr, "Exiting\n");
	exit(1);
}

static void putInHash(int index, int node) {
	hashtype **oldHashPt = &(hash[LOOKUP(&nodeLog[node * units])]);
	hashtype **oldChain = &(chain[index]);
	hashtype *htemp = *oldHashPt;
	hashtype *ctemp = *oldChain;

	if (!(htemp = (hashtype *) malloc(sizeof (hashtype)))) {
		memoryError("adding to hash table");
	}
	htemp->node = node;
	if (*oldHashPt == NULL) {
		*oldHashPt = htemp;
		htemp->next = NULL;
		/* nc++; */
	} else {
		htemp->next = *oldHashPt;
		*oldHashPt = htemp;
		c++;
	}
	if (!(ctemp = (hashtype *) malloc(sizeof (hashtype)))) {
		memoryError("adding to chain");
	}
	ctemp->node = node;
	ctemp->next = *oldChain;
	*oldChain = ctemp;
}

static void
putCurrent(int node) {
	LONG *pnodeLog = &(nodeLog[node * units]);
	LONG *psPOT = shrunkPositionOfTile;
	char *ppOT = positionOfTile;
	int j, shift, n = 0;

	for (;;) {
		shift = 0;
		*psPOT = 0;
		for (j = 0; j < numInUnit; j++) {
			*psPOT |= (((LONG) *ppOT) << shift);
			shift += bits;
			ppOT++;
			n++;
			if (n >= areaTiles) {
				*pnodeLog = *psPOT;
				return;
			}
		}
		*pnodeLog = *psPOT;
		psPOT++;
		pnodeLog++;
	}
}

static void
getCurrent(int node) {
	int stack, position, j, shift, n = 0;
	LONG *pnodeLog = &(nodeLog[node * units]);
	LONG *psPOT = shrunkPositionOfTile;
	char *ppOT = positionOfTile;
	char *ptOP = tileOfPosition;

	for (position = 0; position < height; position++)
		for (stack = 0; stack < MAXSTACKS; stack++) {
			*ptOP = EMPTY;
			ptOP++;
		}
	for (;;) {
		shift = 0;
		*psPOT = *pnodeLog;
		for (j = 0; j < numInUnit; j++) {
			*ppOT = (char) TO_LONG((((LONG) *psPOT) >> shift) &
				bitmask);
			tileOfPosition[(int) *ppOT] = numTile(n);
			shift += bits;
			ppOT++;
			n++;
			if (n >= areaTiles)
				return;
		}
		pnodeLog++;
		psPOT++;
	}
}

static void
setPositionOfTile(LONG *psPOT) {
	int stack, j, shift, n = 0;
	char *ppOT = positionOfTile;

	for (;;) {
		shift = 0;
		for (j = 0; j < numInUnit; j++) {
			*ppOT = (char) TO_LONG((((LONG) *psPOT) >> shift) &
				bitmask);
			if (debug > 1)
				(void) printf("positionOfTile %d %d\n",
					numTile(n), *ppOT);
			shift += bits;
			ppOT++;
			n++;
			if (n >= areaTiles)
				return;
		}
		psPOT++;
	}
}

static void
printShrunkStack(LONG *psPOT) {
	int stack, position, j, shift, n = 0;
	char *ppOT = positionOfTile;
	char *ptOP = tileOfPosition;

	for (position = 0; position < height; position++)
		for (stack = 0; stack < MAXSTACKS; stack++) {
			*ptOP = EMPTY;
			ptOP++;
		}
	ptOP = tileOfPosition;
	for (;;) {
		shift = 0;
		for (j = 0; j < numInUnit; j++) {
			*ppOT = (char) TO_LONG((((LONG) *psPOT) >> shift) &
				bitmask);
			if (debug > 1)
				(void) printf("positionOfTile %d %d\n",
					numTile(n), *ppOT);
			tileOfPosition[(int) *ppOT] = numTile(n);
			shift += bits;
			ppOT++;
			n++;
			if (n >= areaTiles)
				goto BREAKOUT;
		}
		psPOT++;
	}
BREAKOUT:
	for (position = 0; position < height; position++)
		for (stack = 0; stack < MAXSTACKS; stack++) {
			(void) printf("%c", tileChar(*ptOP));
			if (stack == MAXSTACKS - 1)
				(void) printf("\n");
			else
				(void) printf(" ");
			ptOP++;
		}
}

static void
printStacks() {
	int position, stack;
	char *ptOP = tileOfPosition;

	for (position = 0; position < height; position++)
		for (stack = 0; stack < MAXSTACKS; stack++) {
			(void) printf("%c", tileChar(*ptOP));
			if (stack == MAXSTACKS - 1)
				(void) printf("\n");
			else
				(void) printf(" ");
			ptOP++;
		}
}

static void
printOtherStacks() {
	int position;
	char *ptOP = tileOfPosition, *temp;

	for (position = 0; position < height; position++) {
		temp = ptOP; temp++; temp++;
		(void) printf("%c ", tileChar(*temp));
		temp = ptOP; temp++;
		(void) printf("%c ", tileChar(*temp));
		(void) printf("%c\n", tileChar(*ptOP));
		ptOP++; ptOP++; ptOP++;
	}
}

#ifdef DEBUG
static void
printKinds() {
	int tile, kind;
	char *ppOT = positionOfTile;

	for (tile = 0; tile < tiles; tile++)
		for (kind = 0; kind < KINDS; kind++) {
			(void) printf("%d", *ppOT);
			if (kind == KINDS - 1)
				(void) printf("\n");
			else
				(void) printf(" ");
			ppOT++;
		}
}

static void
printOtherKinds() {
	int tile, kind, value;
	char *ppOT = positionOfTile;

	for (tile = 0; tile < tiles; tile++)
		for (kind = 0; kind < KINDS; kind++) {
			value = (*ppOT / MAXSTACKS) * MAXSTACKS + MAXSTACKS -
				1 - *ppOT % MAXSTACKS;
			(void) printf("%d", value);
			if (kind == KINDS - 1)
				(void) printf("\n");
			else
				(void) printf(" ");
			ppOT++;
		}
}
#endif

static void
resetTiles() {
	int tile, unit;

	if (nodeLog)
		(void) free((void *) nodeLog);
	if (!(nodeLog = (LONG *) calloc(size * units,
			sizeof (LONG)))) {
		memoryError("nodeLog");
	}

	if (!(hash = (hashtype **) calloc(hashSize,
			sizeof (hashtype *)))) {
		memoryError("hash table");
	}

	if (positionOfTile)
		(void) free((void *) positionOfTile);
	if (!(positionOfTile = (char *) calloc(areaTiles,
			sizeof (char)))) {
		memoryError("positionOfTile");
	}

	if (tileOfPosition)
		(void) free((void *) tileOfPosition);
	if (!(tileOfPosition = (char *) calloc(areaPositions,
			sizeof (char)))) {
		memoryError("tileOfPosition");
	}

	if (shrunkPositionOfTile)
		(void) free((void *) shrunkPositionOfTile);
	if (!(shrunkPositionOfTile = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("shrunkPositionOfTile");
	}

	if (colorPositionOfTile)
		(void) free((void *) colorPositionOfTile);
	if (!(colorPositionOfTile = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("colorPositionOfTile");
	}

	if (mirrorPositionOfTile)
		(void) free((void *) mirrorPositionOfTile);
	if (!(mirrorPositionOfTile = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("mirrorPositionOfTile");
	}

	if (mirrorColorPositionOfTile)
		(void) free((void *) mirrorColorPositionOfTile);
	if (!(mirrorColorPositionOfTile = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("mirrorColorPositionOfTile");
	}

	if (startTiles)
		(void) free((void *) startTiles);
	if (!(startTiles = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("startTiles");
	}
	if (finishTiles)
		(void) free((void *) finishTiles);
	if (!(finishTiles = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("finishTiles");
	}
	if (startMidTiles)
		(void) free((void *) startMidTiles);
	if (!(startMidTiles = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("startMidTiles");
	}
	if (finishMidTiles)
		(void) free((void *) finishMidTiles);
	if (!(finishMidTiles = (LONG *) calloc(units,
			sizeof (LONG)))) {
		memoryError("finishMidTiles");
	}

	if (move)
		(void) free((void *) move);
	if (!(move = (int *) calloc(size, sizeof (int)))) {
		memoryError("move");
	}
	if (backPoss)
		(void) free((void *) backPoss);
	if (!(backPoss = (char *) calloc(size, sizeof (char)))) {
		memoryError("backPoss");
	}
	if (backPath)
		(void) free((void *) backPath);
	if (!(backPath = (int *) calloc(size, sizeof (int)))) {
		memoryError("backPath");
	}
	for (tile = 0; tile < tiles; tile++) {
		setStart((2 * tile), ((tile + blanks) *  MAXSTACKS));
		setStart((2 * tile + 1), ((tile + blanks) * MAXSTACKS + 2));
	}
	for (tile = 0; tile < tiles; tile++) {
		setFinish((2 * tile), ((tile + blanks) * MAXSTACKS + 2));
		setFinish((2 * tile + 1), ((tile + blanks) *  MAXSTACKS));
	}
	if (tiles > 1) {
		for (tile = 0; tile < ((tiles < 4) ? (tiles - 2) : (tiles - 3)); tile++) {
			setStartMid((2 * tile), (2 * tile * MAXSTACKS + 1));
			setStartMid((2 * tile + 1), ((2 * tile + 1) * MAXSTACKS + 1));
		}
	}
	if (tiles > 1 && tiles < 4) {
		tile = tiles - 2;
		setStartMid((2 * tile + 1), (2 * tile * MAXSTACKS + 1));
		setStartMid((2 * tile), ((2 * tile + 1) * MAXSTACKS + 1));
	}
	if (tiles > 3) {
		tile = tiles - 3;
		setStartMid((2 * tile + 1), (2 * tile * MAXSTACKS + 1));
		setStartMid((2 * tile), ((2 * tile + 1) * MAXSTACKS + 1));
		tile = tiles - 2;
		setStartMid((2 * tile), (2 * tile * MAXSTACKS + 1));
		setStartMid((2 * tile + 1), ((2 * tile + 1) * MAXSTACKS + 1));
	}
	if (tiles == 1) {
		setStartMid(0, 0);
		setStartMid(1, 2);
	} else {
		tile = tiles - 1;
		setStartMid((2 * tile), ((2 * tile - 1) * MAXSTACKS));
		setStartMid((2 * tile + 1), ((2 * tile - 1) * MAXSTACKS + 2));
	}
	if (debug) {
		(void) printf("startMid\n");
		printShrunkStack(startMidTiles);
		(void) printf("\n");
	}
	if (tiles > 1) {
		for (tile = 0; tile < ((tiles < 4) ? (tiles - 1) : (tiles - 3)); tile++) {
			setFinishMid((2 * tile), (2 * tile * MAXSTACKS + 1));
			setFinishMid((2 * tile + 1), ((2 * tile + 1) * MAXSTACKS + 1));
		}
	}
	if (tiles > 3) {
		for (tile = tiles - 3; tile < tiles - 1; tile++) {
			setFinishMid((2 * tile + 1), (2 * tile * MAXSTACKS + 1));
			setFinishMid((2 * tile), ((2 * tile + 1) * MAXSTACKS + 1));
		}
	}
	if (tiles == 1) {
		setFinishMid(0, 2);
		setFinishMid(1, 0);
	} else {
		tile = tiles - 1;
		setFinishMid((2 * tile + 1), ((2 * tile - 1) * MAXSTACKS));
		setFinishMid((2 * tile), ((2 * tile - 1) * MAXSTACKS + 2));
	}
	if (debug) {
		(void) printf("finishMid\n");
		printShrunkStack(finishMidTiles);
		(void) printf("\n");
	}
	for (unit = 0; unit < units; unit++) {
		if (start <= 0)
			shrunkPositionOfTile[unit] = startTiles[unit];
		else if (start == 1)
			shrunkPositionOfTile[unit] = startMidTiles[unit];
		else if (start == 2)
			shrunkPositionOfTile[unit] = finishMidTiles[unit];
		else
			shrunkPositionOfTile[unit] = finishTiles[unit];
	}
	setPositionOfTile(shrunkPositionOfTile);
	putCurrent(currNode);
	if (debug) {
		for (unit = 0; unit < units; unit++)
			(void) printf("start %x\n", startTiles[unit]);
		for (unit = 0; unit < units; unit++)
			(void) printf("finish %x\n", finishTiles[unit]);
		for (unit = 0; unit < units; unit++)
			(void) printf("startMid %x\n", startMidTiles[unit]);
		for (unit = 0; unit < units; unit++)
			(void) printf("finishMid %x\n", finishMidTiles[unit]);
	}
	getCurrent(currNode);
	move[0] = 0;
	backPoss[0] = -1; /* needed for first pass */
	loop = 0;
	hashNewIndex = 0;
	c = 0;
	/* nc = 0;
	fc = 0;
	fnc = 0; */
}

static int
checkSolved() {
	int unit;
	LONG *psPOT = shrunkPositionOfTile;
	LONG *pfT = NULL;

	if (finish <= 0)
		pfT = startTiles;
	else if (finish == 1)
		pfT = startMidTiles;
	else if (finish == 2)
		pfT = finishMidTiles;
	else
		pfT = finishTiles;
	for (unit = 0; unit < units; unit++) {
		if (*pfT != *psPOT)
			return FALSE;
		pfT++;
		psPOT++;
	}
	return TRUE;
}

static int
moveATile(int fromStack, int toStack) {
	int fromPosition = -1;
	int position;
	char *ptOPF = &(tileOfPosition[fromStack]);
	LONG *psPOT;

	for (position = 0; position < height; position++) {
		if (*ptOPF != EMPTY) {
			fromPosition = position;
			break;
		}
		ptOPF += MEMSTACKS;
	}
	if (fromPosition < 0) {
		return FALSE;
	} else {
		/* request move */
		int toPosition;
		char *ptOPT = &(tileOfPosition[toStack]);

		/* Do not have to check middle stack as we
		 * allow one to go through it. */
		if (*ptOPT != EMPTY)
			return FALSE;
		{
			/* topOfStack and Algorithme rule */
			toPosition = tiles + blanks - 1;

			position = 1;
			ptOPT += MEMSTACKS;
			for (; position <= toPosition; position++) {
				if (*ptOPT != EMPTY) {
					toPosition = (position - 1);
					/* Algorithme rule */
					if (tileNumber(*ptOPF) > tileNumber(*ptOPT)) {
						return FALSE;
					}
					break;
				}
				ptOPT += MEMSTACKS;
			}
			ptOPT -= MEMSTACKS;
		}
		/* slide tile */
		{
			int tile = tileNum(*ptOPF);
			int offset = bits * (tile % numInUnit);
			int index = tile / numInUnit;

			*ptOPT = *ptOPF;
			position = toStack + toPosition * MAXSTACKS;
			positionOfTile[tile] = position;
			*ptOPF = EMPTY;
			psPOT = &(shrunkPositionOfTile[index]);
			*psPOT &= ~((LONG) bitmask << offset);
			*psPOT |= ((LONG) position << offset);
		}
	}
	return TRUE;
}

/* FURTHER */
static void
printLastData(int node) {
	if (backPath[node] > node || (node != 0 && backPath[node] == node)) {
		return;
	}
	getCurrent(node);
	if (node != 0) {
		(void) printf("%3d: %d, %d\n", move[node],
			fromPoss((int) backPoss[node]),
			toPoss[(int) backPoss[node]]);
	}
	printStacks();
}

/* FURTHER */
static void
printLastHash(int index) {
	hashtype **oldChain = &(chain[index]);
	hashtype *ctemp = *oldChain;
	while (ctemp != NULL) {
		printLastData(ctemp->node);
		ctemp = ctemp->next;
	}
}

static void
printData(int node) {
	if (backPath[node] > node || (node != 0 && backPath[node] == node)) {
		return;
	}
	if (node != 0) {
		printData(backPath[node]);
	}
	getCurrent(node);
	if (node != 0) {
		(void) printf("%3d: %d, %d\n", move[node],
			fromPoss((int) backPoss[node]),
			toPoss[(int) backPoss[node]]);
	}
	printStacks();
}

static void
printOtherData(int node, int otherMove) {
	getCurrent(node);
	printOtherStacks();
	if (backPath[node] > node || (node != 0 && backPath[node] == node)) {
		return;
	}
	if (node != 0) {
		(void) printf("%3d: %d, %d\n", otherMove,
			2 - toPoss[(int) backPoss[node]],
			2 - fromPoss((int) backPoss[node]));
		printOtherData(backPath[node], otherMove + 1);
	}
}

static void
printAllData(int node, int otherNode) {
	printData(node);
	(void) printf("~~~~~~~~~\n");
	printOtherData(otherNode, move[node] + 1);
}

static int
getnum(char * message) {
	char buff[80];

	(void) printf("%s: ", message);
	return((int) strtol(gets(buff), 0, 0));
}

static void
readData(int level, int nodes) {
	int t;

#ifdef ARCH128
	arch = 128;
#elif defined( ARCH64 )
	arch = 64;
#else
	arch = 8 * sizeof( bitmask );
	(void) printf("Assuming %d bit architecture.\n", arch);
#endif
	tiles = level;
	while (tiles < 1)
		tiles = getnum((char *) "Enter: # of tile pairs >");
	t = tiles - 1;
	blanks = tiles - 2;
	if (blanks < 0)
		blanks = 0;
	height = tiles + blanks;
	areaTiles = KINDS * tiles;
	areaPositions = MAXSTACKS * height;
	EMPTY = 2 * tiles;
	if (t > ARRAYTWEAKS) {
		(void) printf("Be careful, you need lots of memory.\n");
		inputSize = getnum((char *)
			"Enter: # of nodes to keep in memory >");
		hashSize = getnum((char *)
			"Enter: size of hash >");
	} else {
		if (tiles < NOFULLSOLUTION && nodes == 0)
			inputSize = minSizes[mode][t];
		else {
			/* 3 moves must be kept in memory to progress */
			LONG minMem = 3 * moveSizes[mode][t];

			if (nodes == 0)
				inputSize = minMem;
			else
				inputSize = nodes;
			if (inputSize < minMem) {
				(void) printf("Nodes %lld < %lld, changing\n",
					inputSize, minMem);
				inputSize = minMem;
			}
			if (inputSize < minSizes[mode][t])
				(void) printf(
					"Full path may not be printed\n");
		}
		hashSize = hashSizes[t];
	}
	size = inputSize;
	bits = MINBITS();
	units = MINUNITS();
	numInUnit = MAXNUMINUNIT();
	if (debug)
		(void) printf("size %d, bits %d, units %d, numInUnit %d\n",
			size, bits, units, numInUnit);
	bitmask = (1 << bits) - 1;
	freqmask = (int) (~((1 << FREQBITS) - 1));
	if (units > 2)
		(void) printf("Key for hash may not be efficient.\n");
	if (debug > 1) /* hashless check */
		unitMem = (sizeof (char)) * units;
	resetTiles();
}

/* hashless, debug level 2 check */
static int
isNewNode() {
	int unit, node;
	LONG *psPOT, *pnodeLog;
	LONG *pnodeLogFirst = &(nodeLog[newNode * units]);

	for (node = newNode; node >= newBottom; node--) {
		/* check if already done */
		pnodeLog = pnodeLogFirst;
		psPOT = shrunkPositionOfTile;
		for (unit = 0; unit < units; unit++) {
			if (*pnodeLog != *psPOT)
				goto CONTINUE; /* I hate gotos */
			pnodeLog++;
			psPOT++;
		}
		return FALSE;
CONTINUE:	pnodeLogFirst -= unitMem;
	}
	return TRUE;
}

static void
searchBreadthFirst() {
    int localMoveSize = 0, moveSize = 0;
    int otherNode;
    char m;
    int foundSolution = 0;

    for (;;) {
	for (m = 0; m < (char) POSSIBILITIES; m++) {
		if (toPoss[(int) backPoss[currNode]] == fromPoss((int) m))
			continue; /* move could be shorter, so ignore */
		if (moveATile(fromPoss((int) m), toPoss[(int) m])) {
			if (!(inHash())) {
				if (debug > 1 && !isNewNode()) {
					/* hashless check */
					(void) printf("got node:%d %ld\n",
						newNode + 1,
						*shrunkPositionOfTile);
				}
				newNode++;
				if (newNode >= size) {
					newNode = 0;
				}
				if ((newNode & freqmask) == newNode) {
					if (newNode == 0)
						(void) printf("loop:%d; ",
							loop + 1);
					(void) printf("node:%d; move:%d; collisions:%lld\n",
						newNode, move[currNode] + 1, c);
				}
				backPoss[newNode] = m;
				backPath[newNode] = currNode;
				move[newNode] = move[currNode] + 1;
				localMoveSize++;
				putCurrent(newNode); /* save node */
				putInHash(hashNewIndex, newNode);
			    if (mode != FURTHER) {
				setMirrorImage(shrunkPositionOfTile,
					mirrorPositionOfTile);
				if (((mode == MEET_AT_MIDDLE) &&
					((otherNode = checkInHash(
					colorPositionOfTile,
					mirrorPositionOfTile)) != -1) &&
					(!foundSolution || (foundSolution >
					move[newNode] + move[otherNode]))) ||
					((mode != MEET_AT_MIDDLE) &&
					(checkSolved())))
				{
				    if (mode == MEET_AT_MIDDLE) {
					(void) printf("loops: %d, size: %d, nodes: %d; moves: %d, otherMoves: %d\n",
		loop, size, newNode + 1, move[newNode], move[otherNode]);
					printAllData(newNode, otherNode);
					(void) printf("loops: %d, size: %d, nodes: %d; moves: %d, otherMoves: %d\n",
		loop, size, newNode + 1, move[newNode], move[otherNode]);
				    } else {
					(void) printf("loops: %d, size: %d, nodes: %d; moves: %d\n",
		loop, size, newNode + 1, move[newNode]);
					printData(newNode);
				    }
					(void) printf("total nodes: %d * %d + %d\n",
						loop, size, newNode + 1);
					(void) printf("total nodes when found solution: %lld\n",
						(long long) loop * size + newNode + 1);
					(void) printf("moveSize: %d\n",
						moveSize);
					(void) printf("collisions: %lld\n",
						c);
					foundSolution = move[newNode] +
						((mode == MEET_AT_MIDDLE) ?
						move[otherNode] : 0);
					(void) printf("total moves: %d\n",
						foundSolution);
					if (mode == STRAIGHT)
						exit(0);
					/* Do not exit as may find a better solution */
					/* Done when number of moves for mode progresses */
				}
			    }
			}
			getCurrent(currNode); /* backout */
		}
	}
	if (newNode != currNode) {
		if (move[currNode] < move[(currNode + 1) % size]) {
			if (mode != FURTHER && foundSolution) {
				(void) printf("total nodes when complete: %lld\n",
					(long long) loop * size + newNode + 1);
				exit(0);
			}
			/* if (move[newNode] == 12)
				printLastHash(hashNewIndex); */
			if (debug > 1) { /* hashless check */
				newBottom = move[(currNode + 1) % size];
			}
			if (localMoveSize > moveSize) {
				moveSize = localMoveSize;
			}
			localMoveSize = 0;
			hashNewIndex = (hashNewIndex + 1) % CHAINS;
			freeHash(hashNewIndex);
		}
		currNode++;
		if (currNode >= size) {
			currNode = 0;
			loop++;
		}
		getCurrent(currNode);
	} else {
		if (mode == FURTHER)
			(void) printf("Farthest point\n");
		else if (!foundSolution)
			(void) printf("No solution\n");
		(void) printf("loops: %d, size: %d, nodes: %d; moves: %d\n",
			loop, size, newNode + 1, move[newNode]);
		if (mode == FURTHER)
			printData(newNode);
		if (size < minSizes[mode][tiles - 1])
			(void) printf("total nodes: %d\n",
				loop * size + newNode + 1);
		else
			(void) printf("total nodes: %d * %d + %d\n",
				loop, size, newNode + 1);
		(void) printf("moveSize: %d\n", moveSize);
		(void) printf("collisions: %lld\n", c);
		(void) printf("total moves: %d\n", move[newNode]);
		if (mode == FURTHER) {
			printLastHash((hashNewIndex + CHAINS - 1) % CHAINS);
		}
		exit(0);
	}
    }
}

static void
usage()
{
	(void) fprintf(stderr,
		"usage: \"algorithme [-d debugLevel] [-m mode] [-n nodesNumber] [-s start] [-f finish] numTilePairs\"\n");
	(void) fprintf(stderr, "\tmode:\t0 meet in middle (default)\n");
	(void) fprintf(stderr, "\t\t1 match finish point\n");
	(void) fprintf(stderr, "\t\t2 find furthest point\n");
	(void) fprintf(stderr, "\tstart:\tposition to start 0-3, default 0\n");
	(void) fprintf(stderr, "\tfinish:\tposition to finish 0-3, default 3\n");
	exit(1);
}

int
main(int argc, char **argv) {
	int level = 0;
	int nodes = 0;

	argc--;
	argv++;
	while (argc) {
		if (argv[0][0] == '-') {
			switch (argv[0][1]) {
			case '\0':
				usage();
				break;
			case 'd':
				if (argc > 1) {
					argc--;
					argv++;
					debug = atoi(argv[0]);
					break;
				} else {
					usage();
				}
			case 'f':
				if (argc > 1) {
					argc--;
					argv++;
					finish = atoi(argv[0]);
					break;
				} else {
					usage();
				}
			case 'm':
				if (argc > 1) {
					argc--;
					argv++;
					mode = atoi(argv[0]);
					if (mode < 0 || mode >= MODES)
						mode = MEET_AT_MIDDLE;
					break;
				} else {
					usage();
				}
			case 'n':
				if (argc > 1) {
					argc--;
					argv++;
					nodes = atoi(argv[0]);
					break;
				} else {
					usage();
				}
			case 's':
				if (argc > 1) {
					argc--;
					argv++;
					start = atoi(argv[0]);
					break;
				} else {
					usage();
				}
			default:
				usage();
			}
		} else
			level = atoi(argv[0]);
		argc--;
		argv++;
	}
	readData(level, nodes);

#if !defined( __MSDOS__ ) && !defined( __CYGWIN__ )
	/* This gets rid of the unwanted buffering in UNIX */
	(void) setlinebuf(stdout);
#endif

	if (!checkSolved()) {
		searchBreadthFirst();
	} else
		printData(0);
#ifdef VMS
	return 1;
#else
	return 0;
#endif
}