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

/***************************************************
A data structure for the entries of the path matrix:
***************************************************/

#define INFINITY 0x7ffffff /* A VERY large number of edge crossings */

typedef struct
	{
	int pathLength; /* Length of path, measured in number of crossed edges */
	int nextPolygon; /* Index of the next polygon to go to along path */
	} MatrixEntry;

/* "Constructor" functions for the MatrixEntry structure: */

MatrixEntry *trivialEntry(MatrixEntry *entry,int polygonIndex)
	{
	entry->pathLength=0;
	entry->nextPolygon=polygonIndex;
	return entry;
	}

MatrixEntry *invalidEntry(MatrixEntry *entry)
	{
	entry->pathLength=INFINITY;
	entry->nextPolygon=-1;
	return entry;
	}

MatrixEntry *setEntry(MatrixEntry *entry,int sPathLength,int sNextPolygon)
	{
	entry->pathLength=sPathLength;
	entry->nextPolygon=sNextPolygon;
	return entry;
	}

MatrixEntry *copyEntry(MatrixEntry *dest,const MatrixEntry *source)
	{
	dest->pathLength=source->pathLength;
	dest->nextPolygon=source->nextPolygon;
	return dest;
	}

/* Access functions for the MatrixEntry structure: */

int isValid(MatrixEntry* entry)
	{
	return entry->pathLength!=INFINITY;
	}

/* Algebraic operations on the MatrixEntry structure: */

MatrixEntry *multEntries(MatrixEntry *result,const MatrixEntry *e1,const MatrixEntry *e2)
	{
  /* In this function we have to fudge a bit to handle */
  /* both undefined and trivial paths correctly: */
  if(e1->pathLength==INFINITY||e2->pathLength==INFINITY)
    {
    /* If one entry is undefined, the result is undefined, too: */
		setEntry(result,INFINITY,-1);
    }
  else if(e1->pathLength==0)
    {
    /* Entry e1 is trivial, so take only e2: */
    copyEntry(result,e2);
    }
  else if(e2->pathLength==0)
    {
    /* Entry e2 is trivial, so take only e1: */
    copyEntry(result,e1);
    }
  else
    {
    /* Concatenate e1 and e2 by adding their lengths and going */
    /* in direction of entry e1 first: */
		setEntry(result,e1->pathLength+e2->pathLength,e1->nextPolygon);
    }
  return result;
	}

MatrixEntry *addEntries(MatrixEntry *result,const MatrixEntry *e1,const MatrixEntry *e2)
  {
	/* Return the shorter path of the two: */
  if(e1->pathLength<=e2->pathLength)
		copyEntry(result,e1);
  else
		copyEntry(result,e2);
	return result;
  }


/*******************************
Define a data type for matrices:
*******************************/

/* We only consider a special case here, to make life easier :) : */

#define NUMPOLYS 11
typedef MatrixEntry Matrix[NUMPOLYS][NUMPOLYS];

/* "Constructors" of the Matrix type: */

void copyMatrix(Matrix dest,const Matrix source)
	{
	memcpy(dest,source,sizeof(Matrix));
	}

/********************************************************************
The main program creates a map and calculates its transitive closure:
********************************************************************/

/* The map given here is the one from the web page (the indices are entered
   1-based to make things more consistent; this is changed on-the-fly,
	 as the algorithm itself uses 0-based indices): */

#define NUMPAIRS 11
int polyPairs[NUMPAIRS][2]={{1,2},{2,3},{3,7},{3,10},{4,5},{4,11},{6,7},{6,11},
                            {7,8},{8,9},{10,11}};

/* Two matrices to hold intermediate results: */

Matrix matrix1;
Matrix matrix2;

void initMatrix(Matrix matrix)
	{
	int i,j;
	
	/* Initialize the matrix to all trivial/invalid entries
	   (that is, to the identity matrix, algebraically speaking): */
	for(i=0;i<NUMPOLYS;++i)
		for(j=0;j<NUMPOLYS;++j)
			{
			if(i==j) /* This is a trivial entry! */
				trivialEntry(&matrix[i][j],i);
			else /* This is an invalid entry! */
				invalidEntry(&matrix[i][j]);
			}
	
	/* Enter all polygon pairs into the matrix (and convert indices to 0-based): */
	for(i=0;i<NUMPAIRS;++i)
		{
		setEntry(&matrix[polyPairs[i][0]-1][polyPairs[i][1]-1],1,polyPairs[i][1]-1);
		setEntry(&matrix[polyPairs[i][1]-1][polyPairs[i][0]-1],1,polyPairs[i][0]-1);
		}
	}

void multiplyMatrices(Matrix result,const Matrix mat1,const Matrix mat2)
	{
	int i,j,k;
	MatrixEntry mult,sum;
	
	/* Use the standard matrix multiplication algorithm: */
	for(i=0;i<NUMPOLYS;++i)
		for(j=0;j<NUMPOLYS;++j)
			{
			invalidEntry(&sum); /* "sum = 0" */
			for(k=0;k<NUMPOLYS;++k)
				addEntries(&sum,&sum,multEntries(&mult,&mat1[i][k],&mat2[k][j])); /* "sum += mat1[i][k] * mat2[k][j]" */
			copyEntry(&result[i][j],&sum); /* "result[i][j] = sum" */
			}
	}

void printMatrix(Matrix mat)
	{
	int i,j;
	
	/* Print all matrix entries, converting to 1-based for output: */
	for(i=0;i<NUMPOLYS;++i)
		{
		for(j=0;j<NUMPOLYS;++j)
			{
			if(mat[i][j].pathLength<INFINITY)
				printf(" %2d/%2d ",mat[i][j].pathLength,mat[i][j].nextPolygon+1);
			else
				printf("       ");
			}
		printf("\n\n");
		}
	}

void printPath(Matrix mat,int from,int to)
	{
	int currentPoly;
	
	/* Convert query pair to 0-based: */
	--from;
	--to;
	
	/* Print all polygons visited along the path: */
	if(isValid(&mat[from][to])) /* Is there a path connecting those two? */
		{
		printf("Path of length %d: ",mat[from][to].pathLength);
		currentPoly=from;
		while(currentPoly!=to)
			{
			printf("%2d, ",currentPoly+1);
			currentPoly=mat[currentPoly][to].nextPolygon;
			}
		printf("%2d\n",currentPoly+1);
		}
	else
		printf("There is no valid path from %2d to %2d!\n",from+1,to+1);
	}

int main(void)
	{
	int i,from,to;
	
	/* Initialize (and print) the original connectivity matrix: */
	initMatrix(matrix1);
	printMatrix(matrix1);
	getchar();
	
	/* Calculate the transitive closure by repeated squaring: */
	for(i=1;i<NUMPOLYS;i*=2)
		{
		multiplyMatrices(matrix2,matrix1,matrix1);
		copyMatrix(matrix1,matrix2);
		printMatrix(matrix1);
		getchar();
		}
	
	/* Answer the user's path queries: */
	while(1)
		{
		printf("Start polygon index: ");
		scanf("%d",&from);
		if(from<=0)
			break;
		printf("End polygon index: ");
		scanf("%d",&to);
		printPath(matrix1,from,to);
		}
	return 0;
	}

