/*  POTM ENTRY:		FeudalExpress		 			*/
/*  Name:		Vincent Goffin					*/
/*  email:		vjg@mozart.stc.att.com 				*/
/*  compilation:	make FeudalExpress				*/

void srand48();
double drand48();
#define rand_seed(sd) srand48(sd)
#define rand_int(lo, hi) (int)(lo + (hi-lo)*drand48())

#include 
#include 
#include 

#define u_char	unsigned char
typedef u_char	bool_t;
typedef short	dist_t;
typedef short	tour_t;
#define GRIDSZ	(100+50)
#define N	102
#define NN	(N+(N/2))

#define EDGE_INSERTION		1
#define EDGE_RECONNECTION	2

#define TIME_MAX 10		/* time max in minutes of real time */
#define INS_MAX 5		/* max type of insertion moves */
#define SOL_MAX 10		/* keep this many iteration solutions */

bool_t	active_nodes[INS_MAX+1][N];	/* speed up improvement heuristics */

/*
 * Data structures
 */

typedef struct Pt_s
{
	int	x;
	int	y;
} Pt_t;

typedef struct Sol_s
{
	int	len;		/* path length */
	int	sites;		/* different sites visited */
	tour_t	tour[N];
	Pt_t	*steps;
} Sol_t;

Sol_t	Sol[SOL_MAX];
int	nSol = 0;

typedef struct Edge_s
{
	int	h;	/* head */
	int	t;	/* tail */
} Edge_t;

Pt_t Move[] = { { -1,  2 },
		{  1,  2 },
		{  2,  1 },
		{  2, -1 },
		{  1, -2 },
		{ -1, -2 },
		{ -2, -1 },
		{ -2,  1 } };

typedef struct Mv_s
{
	int	d;		/* actual distance covered by this move */
	int	a;		/* 1st principal direction index */
	int	na;		/* number of a steps */
	int	b;		/* 2nd dir index */
	int	nb;		/* number of b's */
	int	mv[8];		/* 8 principal directions of move */
} Mv_t;

typedef struct Op_s
{
	int	type;		/* INSERTION/RECONNECTION */
	int	i;
	int	j;
	int	n;
	} Op_t;

typedef struct Path_s
{
	short	n;
	short	ovlap;
	Pt_t	p;
	struct Path_s	*prev;
} Path_t;

/*
 * Function prototypes
 */

#ifdef __STDC__
#define _ARG_(x) x
#else
#define _ARG_(x) ()
#endif /* __STDC__ */

int readNodes _ARG_((char* file, Pt_t[]));
int Distance _ARG_((Pt_t, Pt_t));
void Construct _ARG_((dist_t[][N], tour_t[]));
void MinTree _ARG_((dist_t[][N], Edge_t[]));
int OddNodes _ARG_((int, Edge_t*, int[]));
int PerfMatch _ARG_((int n_edge, int n_set, int set[], dist_t[][N], Edge_t*));
void Euler _ARG_((int n_edge, Edge_t edges[], tour_t eTour[]));
void Hamilton _ARG_((int n_edge, tour_t eTour[], tour_t hTour[]));

int tourLen _ARG_((tour_t[], dist_t[][N]));
int Improve _ARG_((tour_t[], dist_t[][N], int));
int edgeInsert _ARG_((dist_t[][N], int n, Op_t* pop, int savD, int minD));
int edgeReconnect _ARG_((dist_t[][N], int n, Op_t* pop, int savD, int minD));
void applyOperation _ARG_((tour_t* bTour, tour_t* cTour, Op_t* pop));

int storeSol _ARG_((tour_t kTour[], dist_t[][N], int *is_dup));
void breedSol _ARG_((dist_t[][N]));
int breedTour _ARG_((int, int, tour_t[], dist_t[][N]));

int findSteps _ARG_((dist_t[][N], Sol_t* keep));
void MinOvlapPath _ARG_((Pt_t, Pt_t, Pt_t*, int));
Path_t** MinOvlapRecurse _ARG_((Pt_t*, int, Pt_t*, int, int));

void printSol _ARG_((Sol_t* keep));
int timeUp(double f);

/*
 * Global variables
 */

Pt_t	pt[N];
short	Grid[GRIDSZ][GRIDSZ];
Path_t*	pathMap[GRIDSZ][GRIDSZ];
int	nNode, nIter;
int	randSeed = 1995;
int	bestD = INT_MAX, bestC = -1;


void main(argc, argv)
int	argc;
char	*argv[];
{
	int	i, j, k, D, bestI;
	tour_t	hTour[N];		/* Hamiltonian tour */
	dist_t	dist[N][N];

	/* datafile */
	nNode = readNodes(argv[argc-1], pt);

	/* symmetric distance matrix (0 on diagonal) */
	for (i = 0; i < nNode; ++i)
	{
		for (j = 0; j < i; ++j)
			dist[j][i] = dist[i][j] = Distance(pt[i], pt[j]);
		dist[i][i] = 0;
	}

	rand_seed((long)randSeed);

	/*
	 * Random independent solutions
	 */
	for (nIter = 1; ; ++nIter)
	{
		Construct(dist, hTour);
		D = Improve(hTour, dist, rand_int(1, INS_MAX+1));
		if (nSol < SOL_MAX || D < Sol[nSol-1].len)
		{
			int	is_dup = 0;
			(void) storeSol(hTour, dist, &is_dup);
			if (D < bestD)
				bestD = D;
		}
		if (timeUp(.7))
			break;
	}

	/*
	 * Breed best of random solutions
	 */
	breedSol(dist);

	/*
	 * Minimize overlap in best of breed
	 */
	bestI = 0;
	for (i = 0; i < nSol && Sol[i].len == Sol[0].len; ++i)
	{
		int C = findSteps(dist, Sol + i);
		if (C > bestC)
		{
			bestC = C;
			bestI = i;
		}
		if (C == Sol[bestI].len)
			break;	/* can't do better */
	}

	/*
	 * Print best solution
	 */
	printSol(Sol + bestI);
	exit(0);
}

readNodes(file, pt)
char	*file;
Pt_t	pt[];
{
	FILE	*fp;
	int	n = 0;

	/* origin (not part of the data) */
	pt[0].x = 0, pt[0].y = 0;
	++n;
	/* other points */
	if (!(fp = fopen(file, "r")))
		exit(1);
	while  (fscanf(fp, "%d %d", &pt[n].x, &pt[n].y) == 2)
		++n;
	if (pt[n-1].x < 0)
		--n;
	fclose(fp);
	return n;
}

int Distance(p, q)
/*
 * Returns distance between points p and q;
 */
Pt_t	p, q;
{
	double	slope, eps = 1e-8;
	int	Dx, Dy, D;		/* distance */
	int	csign = 1, m, n;

	/*
	 * For each pair of points there are 2 (of 8) principal
	 * directions to build the shortest path.
	 * They bracket the actual slope.
	 */

	Dx = q.x - p.x;
	Dy = q.y - p.y;
	/* Dx can always be made positive without loss of generality */
	if (Dx < 0)
	{
		Dx = -Dx, Dy = -Dy;
		csign = -1;
	}

	slope = Dx? ((double)Dy)/Dx: ((Dy > 0)? 99: -99);

	if (slope > 2)
	{
		m = -2*Dx + Dy;
		n = 2*Dx + Dy;
		D = m/4 + n/4;
		if (D > 0)
		{
			D += n%4;
			switch (n%4)
			{
			case 1:
			case 3:
				if (m >= 4)
					m -= 4;
				else
					n -= 4;
				break;
			}

		}
		else
		{
			switch (n%4)
			{
			case 1: D = 3; break;
			case 2: D = 2; break;
			case 3: D = 3; break;
			}
		}
	}
	else if (2 - eps < slope && slope < 2 + eps)
	{
		D = Dx;
	}
	else if (slope > .5)
	{
		m = -Dx + 2*Dy;
		n = 2*Dx - Dy;
		D = m/3 + n/3;
		if (D > 0)
		{
			if (n%3)
				D += 2;
			if (n%3 == 2)
			{
				if (m >= 3)
					m -= 3;
				else
					n -= 3;
			}
		}
		else
		{
			if (n%3 == 1)
			{
				if ((p.x == 0 && p.y == 0)
				    || (q.x == 0 && q.y == 0))
					D = 4;
				else
					D = 2;
			}
			else if (n%3 == 2)
				D = 4;
		}
	}
	else if (.5 - eps < slope && slope < .5 + eps)
	{
		D = Dx/2;
	}
	else if (slope > -.5)
	{
		m = Dx + 2*Dy;
		n = Dx - 2*Dy;
		D = m/4 + n/4;
		if (D > 0)
		{
			D += n%4;
			switch (n%4)
			{
			case 1:
			case 3:
				if (m >= 4)
					m -= 4;
				else
					n -= 4;
				break;
			}

		}
		else
		{
			switch (n%4)
			{
			case 1: D = 3; break;
			case 2: D = 2; break;
			case 3: D = 3; break;
			}
		}
	}
	else if (-.5 - eps < slope && slope < -.5 + eps)
	{
		D = Dx/2;
	}
	else if (slope > -2)
	{
		m = 2*Dx + Dy;
		n = -(Dx + 2*Dy);
		D = m/3 + n/3;
		if (D > 0)
		{
			if (n%3)
				D += 2;
			if (n%3 == 2)
			{
				if (m >= 3)
					m -= 3;
				else
					n -= 3;
			}
		}
		else
		{
			switch (n%3)
			{
			case 1: D = 2; break;
			case 2: D = 4; break;
			}
		}
	}
	else if (-2 - eps < slope && slope < -2 + eps)
	{
		D = Dx;
	}
	else			/* slope < -2 */
	{
		m = 2*Dx - Dy;
		n = -(2*Dx + Dy);
		D = m/4 + n/4;
		if (D > 0)
		{
			D += n%4;
			switch (n%4)
			{
			case 1:
			case 3:
				if (m >= 4)
					m -= 4;
				else
					n -= 4;
				break;
			}
		}
		else
		{
			switch (n%4)
			{
			case 1: D = 3; break;
			case 2: D = 2; break;
			case 3: D = 3; break;
			}
		}
	}
	return D;
}

void Construct(dist, hTour)
/*
 * Construct a reasonable Hamiltonian tour
 * using a simplified Christofides heuristic.
 */
dist_t	dist[][N];
tour_t	hTour[];
{
		tour_t	eTour[2*N];
		int	n_odd;
		int	oddNodes[N];
	static	int	n_edge;
	static	Edge_t	edges[2*N];

	/* select a new minimum spanning tree every 10 iterations */
	if (nIter%10 == 1)
	{
		/* Construct a minimum spanning tree */
		MinTree(dist, edges);
		n_edge = nNode - 1;

		/* Determine the set of odd degree nodes remaining */
		if ((n_odd = OddNodes(n_edge, edges, oddNodes)) > 0)
		{
			/* Determine a perfect matching on the odd nodes */
			n_edge += PerfMatch(n_edge, n_odd, oddNodes, dist, edges);
		}
	}

	/* Find an Eulerian tour */
	Euler(n_edge, edges, eTour);

	/* Form a Hamiltonian tour */
	Hamilton(n_edge, eTour, hTour);
}

void MinTree(dist, edges)
/*
 * Find a (randomly selected) minimum spanning tree.
 */
dist_t	dist[][N];
Edge_t	edges[];
{
	int	i, j, e, nsav;
	Edge_t	esav[N];
	bool_t	used[N];

	memset(used, 0, sizeof(used));
	used[0] = 1;			/* start from node 0 */

	for (e = 0; e < nNode-1; ++e)	/* tree edges */
	{
		/* find an unused node that is closest to any used node */
		int	d = INT_MAX;
		for (i = 0; i < nNode; ++i)
			if (used[i])
				for (j = 0; j < nNode; ++j)
				{
					if (!used[j] && dist[i][j] <= d)
					{
						if (dist[i][j] < d)
						{
							d = dist[i][j];
							nsav = 0;
						}
						if (dist[i][j] == d)
						{
							esav[nsav].h = i;
							esav[nsav].t = j;
							nsav++;
						}
					}
				}

		/* include the node in the current path */
		if (d < INT_MAX)
		{
			nsav = rand_int(0, nsav);
			edges[e] = esav[nsav];
			used[esav[nsav].t] = 1;
		}
	}
}

int OddNodes(n_edge, edges, oddNodes)
/* an odd node is a node reached by odd number of edges */
int	n_edge;
Edge_t	*edges;
int	oddNodes[];
{
	int	e, n;
	int	degree[N];
	int	n_odd = 0;

	memset(degree, 0, sizeof(degree));
	for (e = 0; e < n_edge; ++e)
	{
		degree[edges[e].h]++;
		degree[edges[e].t]++;
	}

	for (n = 0; n < nNode; ++n)
		if (degree[n] & 1)
			oddNodes[n_odd++] = n;
	return n_odd;
}

int PerfMatch(n_edge, n_set, set, dist, edges)
/*
 * Perfect match (a paring of the nodes) using a simple algorithm
 */
int	n_edge;
int	n_set, set[];
dist_t	dist[][N];
Edge_t	*edges;
{
	int	i, j, k;
	int	d, nsav, isav[N], jsav[N];
	dist_t	dist2[N][N];
	bool_t	used[N];

	memset(used, 0, sizeof(used));

	/* make distance sub-matrix */
	for (i = n_set; --i; )
	{
		dist_t *ip = dist[set[i]];
		for (j = i; j--; )
			dist2[i][j] = dist2[j][i] = ip[set[j]];
	}

	for (k = 0; k < n_set/2; ++k)
	{
		d = INT_MAX;
		for (i = 1; i < n_set; ++i)
		{
			if (!used[i])
				for (j = 0; j < i; ++j)
					if (!used[j])
					{
						if (dist2[i][j] < d)
						{
							nsav = 0;
							d = dist2[i][j];
						}
						if (dist2[i][j] == d)
						{
							isav[nsav] = i;
							jsav[nsav] = j;
							nsav++;
						}
					}
		}
		nsav = rand_int(0, nsav);
		edges[n_edge + k].h = set[isav[nsav]];
		edges[n_edge + k].t = set[jsav[nsav]];
		used[isav[nsav]] = used[jsav[nsav]] = 1;
	}

	return n_set/2;
}

void Euler(n_edge, edges, eTour)
/*
 * Find an Euler tour in an Eulerian graph.
 */
int	n_edge;
Edge_t	edges[];
tour_t	eTour[];
{
	int	i, j, k, e, last_node;
	int	nodes, ncycle, tcycle;
	int	node_path[NN], cycle_beg[N], cycle_end[N];
	bool_t	edge_used[NN], node_used[N];
	int	esav[N], nsav;

	memset(node_used, 0, sizeof(node_used));
	memset(edge_used, 0, sizeof(edge_used));
	memset(eTour,	  0, n_edge*sizeof(tour_t));
	memset(node_path, 0, sizeof(node_path));

	nodes = ncycle = 0;

	last_node = i = node_path[nodes++] = edges[0].h;
	cycle_beg[ncycle++] = nodes - 1;
	node_used[i] = 1;

	/* find an unused edge originating from node last_node */
next_cycle:
	nsav = 0;
	for (e = n_edge; e--; )
		if (!edge_used[e]
		    && (last_node == edges[e].h || last_node == edges[e].t))
			esav[nsav++] = e;

	if (nsav)
	{
		nsav = rand_int(0, nsav);
		e = esav[nsav];

		if (!node_used[edges[e].h])
			last_node = edges[e].h;
		else
			last_node = edges[e].t;
		node_path[nodes++] = last_node;
		edge_used[e] = 1;
		node_used[last_node] = 1;
		goto next_cycle;
	}

	/* the end of a cycle is the same as the beginning */
	--nodes;

	cycle_end[ncycle-1] = nodes;

	/* start a new cycle:
	   search for the first unused edge connected to a used node */
	for (e = 0; e < n_edge; ++e)
	{
		if (!edge_used[e])
		{
			if (node_used[edges[e].h])
			{
				j = edges[e].h;
				k = edges[e].t;
			}
			else if (node_used[edges[e].t])
			{
				j = edges[e].t;
				k = edges[e].h;
			}
			else
				continue;
			/* this e begins a cycle that will have
			   to be 1) traced, 2) inserted to the tour */
			node_path[nodes++] = j;
			cycle_beg[ncycle++] = nodes-1;
			edge_used[e] = node_used[k] = 1;
			node_path[nodes++] = last_node = k;
			goto next_cycle;
		}
	}

	/* form the Euler tour */
	/* combine the identified cycles into one tour */

	for (tcycle = ncycle, ncycle = 1; ncycle < tcycle; ++ncycle)
	{
		int nl = 0, np = 0, ni;

		ni = cycle_beg[ncycle];
		while ((eTour[nl++] = node_path[np++]) != node_path[ni])
			;

		/* skip 1st one */
		for (i = ni+1; i < cycle_end[ncycle]; ++i)
		{
			eTour[nl++] = node_path[i];
			node_path[i] = -1;
		}
		/* now place 1st one */
		eTour[nl++] = node_path[ni];
		node_path[ni] = -1;

		while (np < nodes)
			if (node_path[np++] != -1)
				eTour[nl++] = node_path[np-1];

		/* re-copy */
		for (i = 0; i < nodes; ++i)
			node_path[i] = eTour[i];
	}

	for (i = 0; i < nodes; ++i)
		eTour[i] = node_path[i];
}

void Hamilton(n_edge, eTour, hTour)
/*
 * Form a Hamiltonian tour from an Eulerian tour.
 * A standard algorithm.
 */
int	n_edge;
tour_t	eTour[], hTour[];
{
	int	i, j, k;
	bool_t	used[N];

	memset(used, 0, sizeof(used));
	for (j = k = 0; j < n_edge; ++j)
	{
		if (!used[eTour[j]])
		{
			hTour[k++] = eTour[j];
			used[eTour[j]] = 1;
		}
	}
}

int tourLen(tour, dist)
tour_t	tour[];
dist_t	dist[][N];
{
	int	i = nNode, D = dist[tour[i-1]][tour[0]];
	while (--i)
		D += dist[tour[i-1]][tour[i]];
	return D;
}

int Improve(bTour, dist, ins_max)
/* Improvement heuristics */
tour_t	bTour[];	/* best Tour */
dist_t	dist[][N];
int	ins_max;
{
	int	i, j, n;
	dist_t	dist2[N][N];
	int	oldD, minD;
	tour_t	cTour[N];	/* current Tour */
	Op_t	op;

	oldD = minD = tourLen(bTour, dist);
	memset(active_nodes, 1, sizeof(active_nodes));
	memcpy(cTour, bTour, sizeof(cTour));

	for (i = nNode; i--; )
	{
		dist_t *pd = dist[bTour[i]];
		for (j = i; j--; )
			dist2[i][j] = dist2[j][i] = pd[bTour[j]];
	}

	while (1)
	{
		int savD = minD;

		for (i = nNode; i--; )
			if (bTour[i] != cTour[i])
			{
				dist_t *pd = dist[bTour[i]];
				for (j = ins_max+1; j--;)
					active_nodes[j][i] = 1;
				for (j = nNode; j--; )
					dist2[i][j] = dist2[j][i] = pd[bTour[j]];
			}
		memcpy(cTour, bTour, N*sizeof(tour_t));

		/* node & edge insertions */
		for (n = 0; n < ins_max; ++n)
			minD = edgeInsert(dist2, n, &op, savD, minD);

		/* edge reconnections */
		minD = edgeReconnect(dist2, n, &op, savD, minD);

		if (minD == savD)
			break;

		applyOperation(bTour, cTour, &op);
	}

	return minD;
}

int edgeInsert(dist, n, pop, savD, minD)
/*
 * Node & edge insertions.
 * Remove edge segment (i,i+n) and insert it before j
 * Allow  0 <= n <= ins_max
 *	n = 0 is a node insertion
 *	n > 0 are edge insertions
 */
dist_t	dist[][N];
int	n;
Op_t	*pop;
int	minD, savD;
{

	int	i, j;
	int	cst, cst2, cst3;
	bool_t	*active = active_nodes[n],
		activeim1 = active[nNode-1];

	for (i = 0; i < nNode-n; ++i)
	{
		int	im1 = i? (i-1): (nNode-1),
			ipe = (i + n + 1)%nNode;
		if (!activeim1)
		{
			for (j = 0; j < n+1; ++j)
				if (active[i+j])
					break;
			if (j == n+1 && !active[ipe])
				continue;
		}
		activeim1 = active[i];
		active[i] = 0;
		cst = dist[im1][ipe] - dist[im1][i] - dist[i+n][ipe];
		/*
		 * if removing the edge didn't help
		 * then re-inserting it elsewhere is not going to help
		 */
		if (cst >= dist[i][i+n])
			continue;
		for (j = 0; j < nNode; ++j)
		{
			int	jm1 = j? (j-1): (nNode-1);
			if ((i <= j && j <= i+n) || j == ipe)
				continue;
			cst3 = cst - dist[jm1][j];
			if (cst3 >= 0)
				continue;
			cst2 = cst3 + dist[jm1][i] + dist[i+n][j];
			if (cst2 < 0)
			{
				active[i] = 1;
				if (savD + cst2 < minD)
				{
					pop->type = EDGE_INSERTION;
					pop->i = i;
					pop->j = j;
					pop->n = n;
					minD = savD + cst2;
				}
			}
			if (n > 0)
			{
				/* flip edge before reconnecting */
				cst2 = cst3 + dist[jm1][i+n] + dist[i][j];
				if (cst2 < 0)
				{
					active[i] = 1;
					if (cst2 < 0 && savD + cst2 <= minD)
					{
						pop->type = EDGE_INSERTION;
						pop->i = i;
						pop->j = j;
						pop->n = -n;	/* flip flag */
						minD = savD + cst2;
					}
				}
			}
		}
	}
	return minD;
}

int edgeReconnect(dist, n, pop, savD, minD)
/*
 * Edge reconnections.
 * Remove edges (i,i+1) and (j,j+1) and reconnect them as (i,j) and (i+1,j+1).
 * There's only 1 way to do this w/o breaking the path into separate cycles
 */
dist_t	dist[][N];
int	n;
Op_t	*pop;
int	minD, savD;
{
	int	i, j, cst, cst2;
	bool_t	*active = active_nodes[n];

	for (i = 0; i < nNode-2; ++i)
	{
		if (!active[i] && !active[i+1])
			continue;
		active[i] = 0;
		cst = dist[i][i+1];
	for (j = i+2; j < nNode; ++j)
	{
		int	jp1 = (j+1)%nNode;

		cst2 = dist[i][j] + dist[i+1][jp1] - dist[j][jp1] - cst;
		if (cst2 < 0)
		{
			active[i] = 1;
			if (savD + cst2 < minD)
			{
				pop->type = EDGE_RECONNECTION;
				pop->i = i;
				pop->j = j;
				pop->n = 1;
				minD = savD + cst2;
			}
		}
	}
	}
	return minD;
}

void applyOperation(bTour, cTour, pop)
/*
 * carry out the best reconnection or insertion
 */
tour_t	*bTour, *cTour;
Op_t	*pop;
{
	memcpy(bTour, cTour, N*sizeof(tour_t));
	if (pop->type == EDGE_INSERTION)
	{
		int	i = pop->i;
		int	j = pop->j;
		int	n = pop->n;
		int	y;
		/* remove (i,i+n) and insert before j */
		if (i < j)
		{
			if (n >= 0)	/* reinsert as (i,i+n) */
			{
				for (y = i; y < j-n-1; ++y)
					bTour[y] = cTour[y+n+1];
				for (y = j-n-1; y < j; ++y)
					bTour[y] = cTour[y+i+n-j+1];
			}
			else	/* reinsert as (i+n,i) */
			{
				n = -n;
				for (y = i; y < j-n-1; ++y)
					bTour[y] = cTour[y+n+1];
				for (y = j-n-1; y < j; ++y)
					bTour[y] = cTour[i+j-1-y];
			}
		}
		else	/* i > j */
		{
			if (n >= 0)	/* reinsert as (i,i+n) */
			{
				for (y = j; y < j+n+1; ++y)
					bTour[y] = cTour[y+i-j];
			}
			else		/* reinsert as (i+n,i) */
			{
				n = -n;
				for (y = j; y < j+n+1; ++y)
					bTour[y] = cTour[i+j+n-y];
			}
			/*
			for (y = j+n+1; y < i+n+1; ++y)
				bTour[y] = cTour[y-n-1];
			*/
			memcpy(bTour+j+n+1, cTour+j, (i-j)*sizeof(tour_t));
		}
	}
	else	/* EDGE_RECONNECTION */
	{
		int	i = pop->i;
		int	j = pop->j;
		int	y, z;
		/* reconnect (i,i+1) to (j,j+1) as (i,j) and (i+1,j+1) */
		for (y = i+1, z = j; z >= i+1; --z, ++y)
			bTour[y] = cTour[z];
	}
}

void breedSol(dist)
/*
 * "Breed" the initial solutions to see if a better solution emerges
 */
dist_t	dist[][N];
{
	int	i, j, stable = 0;
	int	D, minD, timedOut = 0;
	tour_t	bTour[N];

	while (!stable && !timedOut)
	{
		stable = 1;
		for (i = 0; i < 2 && !timedOut; ++i)
		{
			minD = tourLen(Sol[i].tour, dist);
			for (j = i+1; j < nSol; ++j)
			{
				D = breedTour(i, j, bTour, dist);
				if (D < bestD)
				{
					int	k, is_dup = 0;
					bestD = D;
					k = storeSol(bTour, dist, &is_dup);
					if (!is_dup)
					{
						stable = 0;
						i = k-1;
						break;
					}
				}
				if (timeUp(.95))
				{
					timedOut = 1;
					break;
				}
			}
		}
	}
}

int breedTour(iT, jT, Tour3, dist)
int	iT, jT;
tour_t	Tour3[];
dist_t	dist[][N];
{
	int	i, j, k;
	int	step = nNode/8;
	int	D, minD;
	tour_t	dTour[2*N], tTour[N];
	bool_t	used[N];

	minD = tourLen(Sol[iT].tour, dist);
	memcpy(dTour, Sol[iT].tour, nNode*sizeof(tour_t));
	memcpy(dTour + nNode, Sol[jT].tour, nNode*sizeof(tour_t));
	for (i = 0; i < nNode; i += step)
	{
		memset(used, 0, sizeof(used));
		for (j = i, k = 0; k < nNode; j = (j+1)%(2*nNode))
		{
			if (!used[dTour[j]])
			{
				tTour[k++] = dTour[j];
				used[dTour[j]] = 1;
			}
		}
		D = Improve(tTour, dist, INS_MAX);
		if (D < minD)
		{
			minD = D;
			memcpy(Tour3, tTour, nNode*sizeof(tour_t));
		}
	}
	return minD;
}

int storeSol(kTour, dist, is_dup)
tour_t	kTour[];
dist_t	dist[][N];
int	*is_dup;
{
	int	i, j, n, n_dup;
	tour_t	tTour[N];
	int	D = tourLen(kTour, dist);

	*is_dup = 0;
	if (nSol < SOL_MAX)
	{
		Sol[nSol].len = D;
		memcpy(Sol[nSol].tour, kTour, sizeof(tTour));
		Sol[nSol].sites = 0;
		Sol[nSol].steps = 0;
		++nSol;
	}
	else
	{
		/* replace longest */
		int	isav, lmax = 0;
		for (i = 0; i < nSol; ++i)
			if (Sol[i].len > lmax)
			{
				lmax = Sol[i].len;
				isav = i;
			}
		if (D < lmax)
		{
			Sol[isav].len = D;
			memcpy(Sol[isav].tour, kTour, sizeof(tTour));
		}
	}

	/* sort Sol[] by tour length (ascending order) */

	n_dup = 0;
	for (i = 0; i < nSol-1; ++i)
	for (j = i+1; j < nSol; ++j)
		if (Sol[j].len < Sol[i].len)
		{
			n = Sol[i].len;
			Sol[i].len = Sol[j].len;
			Sol[j].len = n;
			memcpy(tTour, Sol[i].tour, sizeof(tTour));
			memcpy(Sol[i].tour, Sol[j].tour, sizeof(tTour));
			memcpy(Sol[j].tour, tTour, sizeof(tTour));
		}
		else if (Sol[j].len == Sol[i].len && Sol[j].len != INT_MAX)
		{
			if (memcmp(Sol[j].tour, Sol[i].tour, nNode*sizeof(tour_t)) == 0)
			{
				Sol[j].len = INT_MAX;
				n_dup++;
			}
		}

	/* all duplicates should have bubbled to the top */
	if (n_dup)
		*is_dup = 1;
	nSol -= n_dup;

	for (i = 0; i < nSol; ++i)
		if (memcmp(Sol[i].tour, kTour, nNode*sizeof(tour_t)) == 0)
			break;
	return i;
}

int findSteps(dist, best)
/*
 * Determine the individual steps of the tour.
 * Store result in best.tour[]
 */
dist_t	dist[][N];
Sol_t	*best;
{
	int	i, j;
	int	pass = 0, old_sites;
	int	D = best->len, cumDist;
	tour_t	*bTour = best->tour;
	Pt_t	*steps;

	steps = (Pt_t *)calloc(D+1, sizeof(Pt_t));
	memset(Grid, 0, sizeof(Grid));

	best->sites = 0;
	do
	{
		old_sites = best->sites;
		cumDist = 0;
		for (i = 0; i < nNode; ++i)
		{
			Pt_t	*step = steps + cumDist;
			int	a = bTour[i];
			int	b = bTour[(i+1)%nNode];
			int	d = dist[a][b];

			if (pass)
				for (j = 0; j < d+1; ++j)
					Grid[step[j].x][step[j].y] -= 1;

			/* visit as many sites as possible */
			MinOvlapPath(pt[a], pt[b], step, d+1);

			for (j = 0; j < d+1; ++j)
				Grid[step[j].x][step[j].y] += 1;

			cumDist += d;
		}

		/* count how many visited sites */
		best->sites = 0;
		for (i = 0; i < GRIDSZ; ++i)
		for (j = 0; j < GRIDSZ; ++j)
			if (Grid[i][j])
				best->sites++;

		best->steps = steps;
		pass++;
	} while (best->sites > old_sites && best->sites < cumDist);

	return best->sites;
}

void printSol(best)
/*
 * print a solution
 */
Sol_t *best;
{
	int	i, j;
	int	D = best->len;
	Pt_t	*steps = best->steps;

	for (i = 0; i < D; i++)
		if (steps[i].x == 0 && steps[i].y == 0)
			break;
	for (j = i; j < D+i; ++j)
		printf("%d %d\n", steps[j%D].x, steps[j%D].y);
	printf("0 0\n");
}

Path_t *new_path(n, ovlap, p, prev)
/* path allocator/constructor */
int	n, ovlap;
Pt_t	p;
Path_t	*prev;
{
	static	int n_paths = 0, i_paths = 0;
	static	Path_t *p_paths[1024];
	Path_t	*path;

	if (n == 0)	/* re-initialization */
	{
		int	i;
		for (i = 0; i < i_paths; ++i)
			free((char*)p_paths[i]);
		i_paths = n_paths = 0;
		return 0;
	}
	if (n_paths == 1024)
		n_paths = 0;
	if (n_paths == 0)
		p_paths[i_paths++] = (Path_t*)malloc(1024*sizeof(Path_t));

	path = p_paths[i_paths-1] + n_paths++;
	path->n = n;
	path->ovlap = ovlap;
	path->p = p;
	path->prev = prev;
	return path;
}

void MinOvlapPath(pa, pb, step, d)
/*
 * Find best path between points p and q.
 */
Pt_t	pa, pb, *step;
int	d;
{
	int	i, n;
	Path_t	**paths, *p, *q, *q2;

	if (pa.x == pb.x && pa.y == pb.y)
		return;

	memset(pathMap, 0, sizeof(pathMap));
	pathMap[pa.x][pa.y] = new_path(1, 0, pa, 0);
	pathMap[pb.x][pb.y] = new_path(-1, 0, pb, 0);

	paths = MinOvlapRecurse(&pa, 1, &pb, 1, 1);
	for (p = paths[0], q = paths[1]; q; q = q2)
	{
		q2 = q->prev;
		q->prev = p;
		p = q;
	}

	if (p->p.x == pa.x && p->p.y == pa.y)
	{
		for (; p; p = p->prev)
			*step++ = p->p;
	}
	else if (p->p.x == pb.x && p->p.y == pb.y)
	{
		step += d;
		for (; p; p = p->prev)
			*--step = p->p;
	}
	else
		exit(2);
	new_path(0, 0, 0, 0);
}

Path_t** MinOvlapRecurse(pSites, np, qSites, nq, n)
/*
 *	return 1 best path
 */
Pt_t	*pSites, *qSites;
int	np, nq, n;
{
		int	i, j, npNew;
		int	contact = 0;
		Pt_t	p, *pSitesNew;
	static	Path_t	*resPaths[2];
		int	resOvlap = INT_MAX;

	/* breadth first search from both ends for a/the meeting point/s */

	pSitesNew = (Pt_t*)calloc(8*np, sizeof(Pt_t));
	npNew = 0;

	for (i = 0; i < np; ++i)
	{
		Pt_t	iP;
		Path_t *iPath, *jPath;

		iP = pSites[i];
		iPath = pathMap[iP.x][iP.y];

		for (j = 0; j < 8; ++j)
		{
			p.x = iP.x + Move[j].x;
			p.y = iP.y + Move[j].y;

			/* path step must be on the allowed board */
			if (p.x<0 || p.y<0 || p.x>=GRIDSZ || p.y>=GRIDSZ)
				continue;

			jPath = pathMap[p.x][p.y];
			if (!contact && (!jPath || jPath->n == n))
			{
				int ovlap = iPath->ovlap;
				/* site is new or was visited during
				   this loop */
				if (Grid[p.x][p.y] != 0)
					ovlap++;
				if (!jPath)
				{
					pSitesNew[npNew++] = p;
					pathMap[p.x][p.y] =
						new_path(n, ovlap, p, iPath);
				}
				else if (jPath->ovlap > ovlap)
				{
					jPath->ovlap = ovlap;
					jPath->prev = iPath;
				}
				/* else: keep current jPath */
			}
			else if (jPath && jPath->n*n < 0)
			{
				/* we met! */
				contact = 1;
				if (jPath->ovlap + iPath->ovlap < resOvlap)
				{
					resOvlap = jPath->ovlap + iPath->ovlap;
					resPaths[0] = iPath;
					resPaths[1] = jPath;
				}
				pSitesNew = 0;
			}
			/* else: site previously visited, ignore it */
		}
	}

	if (!contact) /* switch args */
		MinOvlapRecurse(qSites, nq, pSitesNew, npNew, n<0? 1-n:-(1+n));
	free((char*)pSitesNew);
	return resPaths;
}

int timeUp(f)
double f;
{
	struct tms tbuff;
	times(&tbuff);
	return (tbuff.tms_utime + tbuff.tms_stime)/(CLK_TCK*60) > TIME_MAX*f;
}