/* 
   Date: 1-27-2004 
   Author: Yong Zhang
   This program is based on graph6node.c.
   This program is used to find a weighted acyclic graph with
   9 vertices that satisfies the following property:
   The distances from any two vertices are all distinct and belong
   to the set {1,...,36}.
   This program can be modified to solve the similar problem for
   a bigger input.

   This is the final version of the program. 

*/

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

#define FALSE 0
#define TRUE 1
#define ERROR -2

#define DIM 9 /* the dimension of the matrix */

#define min(x,y) ((x) < (y) ? (x) : (y))
#define max(x,y) ((x) > (y) ? (x) : (y))

int compute_path(int A[DIM][DIM], int B[DIM][DIM]);
int multiply_matrix(int A[DIM][DIM], int B[DIM][DIM]);
int multiply_column(int col_a[DIM], int col_b[DIM]);
int isvalid_matrix(int d[DIM][DIM]);
int isequal_matrix(int A[DIM][DIM], int B[DIM][DIM]);
int nextvalue_matrix(int d[DIM][DIM]);
void assign_value(int d[DIM][DIM], int val, int location);
void initiate_matrix(int d[DIM][DIM]);
void copy_matrix(int A[DIM][DIM], int B[DIM][DIM]);
void print_matrix(int d[DIM][DIM]);



/*
 *********************************************************
 *                                                       *
 *   In addition to the above DIM Macro, the following   *
 *   needs to be modified when scale up the program.     *
 *                                                       *
 *********************************************************
*/



/* there are 36 paths in a graph of nine nodes */
int positions=36;

/* Here is how we match the position and the ijth entry
   of the matrix.

   The upper triangular of the original matrix:
   (0,0) (0,1) (0,2) (0,3) (0,4) (0,5) (0,6) (0,7) (0,8)
         (1,1) (1,2) (1,3) (1,4) (1,5) (1,6) (1,7) (1,8)
               (2,2) (2,3) (2,4) (2,5) (2,6) (2,7) (2,8)
                     (3,3) (3,4) (3,5) (3,6) (3,7) (3,8)
                           (4,4) (4,5) (4,6) (4,7) (4,8)
                                 (5,5) (5,6) (5,7) (5,8)
                                       (6,6) (6,7) (6,8)
                                             (7,7) (7,8)
                                                   (8,8)  

  and the corresponding position:

  # 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 28
                # 30 31 32
                   # 33 34
                      # 35
                         #

*/

/* the pos array serves as a lookup table which maps
   {0,...,35} to ji'th entry of the matrix */ 
int pos[36]={10,20,30,40,50,60,70,80,21,
	     31,41,51,61,71,81,32,42,52,
	     62,72,82,43,53,63,73,83,54,
	     64,74,84,65,75,85,76,86,87};


/* the num array was used by isvalid_matrix() */
int num[37];

/* m1 is our input matrix */
int m1[DIM][DIM]={
  {0,1,-1,-1,-1,-1,-1,-1,-1},
  {1,0,-1,-1,-1,-1,-1,-1,-1},
  {-1,-1,0,-1,-1,-1,-1,-1,-1},
  {-1,-1,-1,0,-1,-1,-1,-1,-1},
  {-1,-1,-1,-1,0,-1,-1,-1,-1},
  {-1,-1,-1,-1,-1,0,-1,-1,-1},
  {-1,-1,-1,-1,-1,-1,0,-1,-1},
  {-1,-1,-1,-1,-1,-1,-1,0,-1},
  {-1,-1,-1,-1,-1,-1,-1,-1,0}
  }; 


/* in every loop we need to save the intermediate matrix */
int m2[DIM][DIM],m3[DIM][DIM],m4[DIM][DIM];
int m5[DIM][DIM], m6[DIM][DIM],m7[DIM][DIM];
int m8[DIM][DIM], result[DIM][DIM];

/* positions of eight edges*/
int p1,p2,p3,p4,p5,p6,p7,p8;
/* weight of eight edges */
int v1,v2,v3,v4,v5,v6,v7,v8;


/* for debug use */
int testval;
int test[DIM][DIM]={
  {0,1,-1,-1,-1,-1,-1,-1,-1},
  {1,0,2,-1,-1,-1,-1,-1,-1},
  {-1,2,0,3,4,-1,-1,-1,-1},
  {-1,-1,3,0,-1,-1,-1,-1,-1},
  {-1,-1,4,-1,0,5,6,-1,-1},
  {-1,-1,-1,-1,5,0,-1,-1,8},
  {-1,-1,-1,-1,6,-1,0,7,-1},
  {-1,-1,-1,-1,-1,-1,7,0,-1},
  {-1,-1,-1,-1,-1,8,-1,-1,0}
};


int test2[DIM][DIM]={
  {0,1,3,5,7,9,11,13,15},
  {1,0,17,19,21,23,25,27,29},
  {3,17,0,31,33,35,2,4,6},
  {5,19,31,0,8,10,12,14,16},
  {7,21,33,8,0,18,20,22,24},
  {9,23,35,10,18,0,26,28,30},
  {11,25,2,12,20,26,0,32,34},
  {13,27,4,14,22,28,32,0,36},
  {15,29,6,16,24,30,34,36,0}
};


int main(int argc, char *argv[])
{
  int i;



  system("date");

  /* for debug use, this is to test if compute_path() and isvalid() works 
  print_matrix(test);
  compute_path(result, test);
  print_matrix(result);
  testval=isvalid_matrix(result);
  printf("%d\n", testval);
  testval=isvalid_matrix(test2);
  printf("%d\n", testval);
  return; */

  /* initiate the num array */
  for (i=0;i<37;i++)
    num[i]=0;


  /* assign the first edge and weight, in fact we don't 
     need to use them, it is included in the input matrix,
     we state them here for the sake of clearance */
  p1=0;
  v1=1;


  /* In fact, p2 only needs to be 1 or 15, v2 is 2 */ 
  p2=1;
  v2=2;

  /* make a new copy of input matrix m1 */
  copy_matrix(m2,m1);

  /* assign v2 to given location p2 */
  assign_value(m2,v2,p2);

  /* when p2=1, v3 has to be 4 */
  v3=4;

  for (p3=1;p3<positions;p3++)
    {
      /* p3 only needs to be 2,9,15 or 21 */
      if (p3!=2 && p3!=9 && p3!=15 && p3!=21)
	continue;

      /* make a new copy of m2 */
      copy_matrix(m3,m2);

      /* assignment the third value */
      assign_value(m3,v3,p3);
      
      /* calculate the next value */
      v4=nextvalue_matrix(m3);

      /* if there is a contradiction, go to next p3 */
      if (v4==ERROR) 
	continue;

      for (p4=1;p4<positions;p4++)
	{
	  /* make sure p3 is a distinct location */
	  if (p4==p2 || p4==p3) 
	    continue;
	  
	  /* make a new copy of m3 */
	  copy_matrix(m4,m3);
	  
	  /* assign the fourth value */
	  assign_value(m4,v4,p4);
	  
	  /* calculate the next value */
	  v5=nextvalue_matrix(m4);
	  
	  /* if there is a contradiction, go to next p4 */
	  if (v5==ERROR) 
	    continue;

	  for (p5=1;p5<positions;p5++)
	    {
	      /* make sure p5 is a distinct location */
	      if (p5==p2 || p5==p3 || p5==p4) 
		continue;
	      
	      /* make a new copy of m3 */
	      copy_matrix(m5,m4);
	      
	      /* assign the fifth value */
	      assign_value(m5,v5,p5);
	      
	      /* calculate the next value */
	      v6=nextvalue_matrix(m5);
	      
	      /* if there is a contradiction, go to next p5 */
	      if (v6==ERROR) 
		continue;
	      
	      for (p6=1;p6<positions;p6++)
		{
		  /* make sure p6 is a distinct location */
		  if (p6==p2 || p6==p3 || p6==p4 || p6==p5) 
		    continue;
		  
		  /* make a new copy of m5 */
		  copy_matrix(m6,m5);
		  
		  /* assign the sixth value */
		  assign_value(m6,v6,p6);
		  
		  /* calculate the next value */
		  v7=nextvalue_matrix(m6);
		  
		  /* if there is a contradiction, go to next p6 */
		  if (v7==ERROR) 
		    continue;
		  
		  for (p7=1;p7<positions;p7++)
		    {
		      /* make sure p7 is a distinct location */
		      if (p7==p2 || p7==p3 || p7==p4 || p7==p5 || p7==p6) 
			continue;
		      
		      /* make a new copy of m6 */
		      copy_matrix(m7,m6);
		      
		      /* assign the seventh value */
		      assign_value(m7,v7,p7);
		      
		      /* calculate the next value */
		      v8=nextvalue_matrix(m7);
		      
		      /* if there is a contradiction, go to next p7 */
		      if (v8==ERROR) 
			continue;
		      
		      for (p8=1;p8<positions;p8++)
			{
			  /* make sure p8 is a distinct location */
			  if (p8==p2 || p8==p3 || p8==p4 || p8==p5 || p8==p6 || p8==p7) 
			    continue;
			  
			  /* make a new copy of m7 */
			  copy_matrix(m8,m7);
			  
			  /* assign the last value */
			  assign_value(m8,v8,p8);
			  
			  /* compute all paths */
			  if (compute_path(result,m8)==ERROR)
			    continue;
			  
			  /* check if it is the desired matrix */
			  if (isvalid_matrix(result)) 
			    {
			      printf("found!\n");
			      print_matrix(m8);
			      print_matrix(result);
			      return;
			    }
			}
		    }
		}
	    }
	}
    }


  system("date");	  
  printf("half\n");

  /* the other value for p2 and v2 */
  p2=15;
  v2=2;

  /* make a new copy of input matrix m1 */
  copy_matrix(m2,m1);

  /* assign v2 to given location p2 */
  assign_value(m2,v2,p2);

  /* when p2=15, v3 has to be 3 */
  v3=3;

  for (p3=1;p3<positions;p3++)
    {
      /* p3 only needs to be 1,3,16,26 */
      if (p3!=1 && p3!=3 && p3!=16 && p3!=26)
	continue;

      /* make a new copy of m2 */
      copy_matrix(m3,m2);

      /* assignment the third value */
      assign_value(m3,v3,p3);
      
      /* calculate the next value */
      v4=nextvalue_matrix(m3);

      /* if there is a contradiction, go to next p3 */
      if (v4==ERROR) 
	continue;

      for (p4=1;p4<positions;p4++)
	{
	  /* make sure p3 is a distinct location */
	  if (p4==p2 || p4==p3) 
	    continue;
	  
	  /* make a new copy of m3 */
	  copy_matrix(m4,m3);
	  
	  /* assign the fourth value */
	  assign_value(m4,v4,p4);
	  
	  /* calculate the next value */
	  v5=nextvalue_matrix(m4);
	  
	  /* if there is a contradiction, go to next p4 */
	  if (v5==ERROR) 
	    continue;

	  for (p5=1;p5<positions;p5++)
	    {
	      /* make sure p5 is a distinct location */
	      if (p5==p2 || p5==p3 || p5==p4) 
		continue;
	      
	      /* make a new copy of m3 */
	      copy_matrix(m5,m4);
	      
	      /* assign the fifth value */
	      assign_value(m5,v5,p5);
	      
	      /* calculate the next value */
	      v6=nextvalue_matrix(m5);
	      
	      /* if there is a contradiction, go to next p5 */
	      if (v6==ERROR) 
		continue;
	      
	      for (p6=1;p6<positions;p6++)
		{
		  /* make sure p6 is a distinct location */
		  if (p6==p2 || p6==p3 || p6==p4 || p6==p5) 
		    continue;
		  
		  /* make a new copy of m5 */
		  copy_matrix(m6,m5);
		  
		  /* assign the sixth value */
		  assign_value(m6,v6,p6);
		  
		  /* calculate the next value */
		  v7=nextvalue_matrix(m6);
		  
		  /* if there is a contradiction, go to next p6 */
		  if (v7==ERROR) 
		    continue;
		  
		  for (p7=1;p7<positions;p7++)
		    {
		      /* make sure p7 is a distinct location */
		      if (p7==p2 || p7==p3 || p7==p4 || p7==p5 || p7==p6) 
			continue;
		      
		      /* make a new copy of m6 */
		      copy_matrix(m7,m6);
		      
		      /* assign the seventh value */
		      assign_value(m7,v7,p7);
		      
		      /* calculate the next value */
		      v8=nextvalue_matrix(m7);
		      
		      /* if there is a contradiction, go to next p7 */
		      if (v8==ERROR) 
			continue;
		      
		      for (p8=1;p8<positions;p8++)
			{
			  /* make sure p8 is a distinct location */
			  if (p8==p2 || p8==p3 || p8==p4 || p8==p5 || p8==p6 || p8==p7) 
			    continue;
			  
			  /* make a new copy of m7 */
			  copy_matrix(m8,m7);
			  
			  /* assign the last value */
			  assign_value(m8,v8,p8);
			  
			  /* compute all paths */
			  if (compute_path(result,m8)==ERROR)
			    continue;
			  
			  /* check if it is the desired matrix */
			  if (isvalid_matrix(result)) 
			    {
			      printf("found!\n");
			      print_matrix(m8);
			      print_matrix(result);
			      return;
			    }
			}
		    }
		}
	    }
	}
    }


  system("date");
  printf("done\n");

   
}



/*
 *********************************************************
 *                                                       *
 *   Routines that scale up.                             *
 *                                                       *
 *********************************************************
*/


/*
  Given the adjancency matrix B as input, this routine will 
  output A the weight of all paths between any two vertices
  using all-pair shortest path algoritm. The return value 
  is used to propagate the error up. If it returns ERROR, 
  there is a cycle in the graph.
*/
int compute_path(int A[DIM][DIM], int B[DIM][DIM])
{
  int i,val;
  copy_matrix(A, B);

  /* multiply matrix DIM-1 times */
  for (i=1;i<DIM;i++)
    {
      val=multiply_matrix(A,B);
      /* if error, then propogate the error up. */
      if (val==ERROR)
	return ERROR;
      else if(val==FALSE)
	break; 
    }

  return TRUE;
}


/*
  Do the matrix multiplication and store the result in A.
  Please note that the multiplication is not the usual
  matrix multiplication. The return value is used to propagate 
  the error up. If it returns ERROR, there is a cycle in the
  graph. If it returns FALSE, no work is done, i.e., A remains
  the same. Otherwise, it returns TRUE.
*/
int multiply_matrix(int A[DIM][DIM], int B[DIM][DIM])
{
  /* declare a tmp matrix to store the result */
  int i,j,tmp[DIM][DIM];

  /* set flag to 0, this flag is to indicate if we did some
   useful work. */
  int flag=FALSE;

  copy_matrix(tmp,A);

  /* calculate the result */
  for (i=0; i<DIM-1; i++)
    for (j=i+1; j<DIM; j++)
      {
	/* since the graph is acyclic, the path between any 
	   two vertices is unique */
	if (tmp[i][j]!=-1) 
	  continue;

	tmp[i][j]=multiply_column(A[i],B[j]);

	/* if error, then propogate the error up, if we did
	 some useful work, set flag to TRUE. */
	if (tmp[i][j]==ERROR) 
	  return ERROR;
	else if(tmp[i][j]!=-1)
	  flag=TRUE;

	tmp[j][i]=tmp[i][j];
      }

  /* copy the result back to A */
  copy_matrix(A,tmp);
  return flag;
}


/*
  This routine is used by multiply_matrix() to calculate
  the multiplication of two rows. (In fact one row and one
  column, but here the matrix is symmetric.)
*/
int multiply_column(int row_a[DIM], int row_b[DIM])
{
  /* initialize result to be some impossible value */
  int i,result=-1;

  /* The multiplication we define is to add the two rows
     component-wise and take the minimum. */
  for (i=0;i<DIM;i++){
    /* If either one is -1, i.e., the two vertices are not 
       connected, then skip. */
    if (row_a[i]<=0 || row_b[i]<=0) continue;

    if (result==-1)
      result=row_a[i]+row_b[i];
    /* if there is a cycle in the graph, return some impossible
       value to indicate error. */
    else
      return ERROR;
  }
  return result;
}


/*
  The routine is used to check if the given matrix is
  what we want, i.e., the upper triangular contains 
  integers from 1 to DIM choose 2.
*/
int isvalid_matrix(int d[DIM][DIM])
{
  int i,j,length=positions+1;

  /* initiate the num array first */
  for (i=0;i<length;i++)
    num[i]=0;

  for (i=0; i<DIM-1; i++)
    for (j=i+1; j<DIM; j++)
      {
	/* if there is no path between i and j, or the weight is bigger
	   desired value, fail */
	if (d[i][j]==-1 || d[i][j]>positions) return FALSE;
	/* if there are two path of the same weight, fail */
	if (num[d[i][j]]!=0) return FALSE;
	/* otherwise put the weight in the array */
	num[d[i][j]]=d[i][j];
      }  
  return TRUE;
}


/*
  This routine is used to check if the two given matrices 
  are the same. 
*/
int isequal_matrix(int A[DIM][DIM], int B[DIM][DIM])
{
  int i,j;

  /* since all matrices are symmetric, we will only
     check the upper triangular part. */
  for (i=0; i<DIM-1; i++)
    for (j=i+1; j<DIM; j++)
      {
	if (A[i][j]!=B[i][j]) return FALSE;
      }
  return TRUE;
}


/*
  This routine is used to next value we should assign
  to the given matrix, which is the smallest integer that
  can not be represented by the given graph. 
*/
int nextvalue_matrix(int d[DIM][DIM])
{
  int i,j;
  int tmp[DIM][DIM];

  /* calculate the weight of all path by multiplying the matrices */
  if (compute_path(tmp,d)==ERROR) 
    return ERROR;

  /* find the next value, which is the smallest integer that can not
     be found in the tmp matrix */
  /* initiate the num array first */
  for (i=0;i<positions+1;i++)
    num[i]=0;

  /* put the positive values in tmp matrix in the num[] array */
  for (i=0; i<DIM-1; i++)
    for (j=i+1; j<DIM; j++)
      {
	if (tmp[i][j]<0)
	  continue;
	/* if we found a large path, fail */
	else if (tmp[i][j]>positions) 
	  return ERROR;
	/* if we find two paths of the same weight, fail */
	else if (num[tmp[i][j]]!=0)
	  return ERROR;
	else
	  num[tmp[i][j]]=tmp[i][j];
      }

  /* return the smallest integer that is not in the tmp matrix */
  for (i=2;i<positions+1;i++)
    {
      if (num[i]==0)
	return i;
    }
  return 0;
}


/*
  This routine is used to assign given value to the 
  given location (and its transpose location) of the 
  matrix. 
*/
void assign_value(int d[DIM][DIM], int val, int location)
{
  int i,j;

  /* decode the location */
  i=pos[location]%10;
  j=(pos[location]-i)/10;

  /* assign value */
  d[i][j]=val;
  d[j][i]=val;
}


/*
  Printout the content of the given matrix. Mainly for 
  purpose of debug.
*/
void print_matrix(int d[DIM][DIM])
{ 
  int i,j;
  for (i=0;i<DIM;i++){
    for (j=0;j<DIM;j++)
      printf("%d   ", d[i][j]);
    printf("\n");
  }
  printf("\n");
}


/*
  This routine is used to initiate the given matrix to 0.
*/
void initiate_matrix(int d[DIM][DIM])
{
  int i,j;
  for (i=0;i<DIM;i++)
    for (j=0;j<DIM;j++)
      d[i][j]=0;
}


/*
  This routine is used to copy matrix B into matrix A.
*/
void copy_matrix(int A[DIM][DIM], int B[DIM][DIM])
{
  int i,j;
  for (i=0;i<DIM;i++)
    for (j=0;j<DIM;j++)
      A[i][j]=B[i][j];
}