/* Determine whether a digraph has cycles.

   Suppose that the vertices of a digraph are numbered 1,2,...,N.
   Let A be the connection matrix of the graph, where

	A(i,j) = 1, if there is a directed edge from vertex j to vertex i
	       = 0, if not.

   DEFINITION. Write j->i if there is a directed path in the digraph leading
   from vertex j to vertex i.

   LEMMA. j->i if and only if A^p(i,j)>0 for some integer power p.

   PROOF. A^p(i,j)>0 if and only if there is a p-step path i->...->j. []

   DEFINITION. An cycle in the digraph is an ordered list of vertices
   C={i1,...,iK} such that j->i for every pair i,j in C.

   DEFINITION. Vertex i is a pure source if there is no vertex j with j->i. 

   DEFINITION. Vertex j is a pure sink if there is no vertex i with j->i. 


   LEMMA.  Vertex i is a pure source if and only if row i of A is zero.
   Vertex j is a pure sink if and only if column j of A is zero.

   PROOF. If row i of A is zero, then row i of A^p will be zero for all p>0.
   Likewise, if column i of A is zero, then column i of A^p will be zero for
   all p>0. []


   THEOREM.  A has a cycle {i1,...,iK} if and only if A/(Id - c*A) has a
   strictly positive submatrix consisting of the rows and columns labeled
   i1,i2,...,iK, where 0<c< 0.5/||A||.

   PROOF. Put B = A/(Id - c*A) = A + A*cA +A*cA*cA +... and note that this
   series of matrices with nonnegative coefficients converges.

   (==>) If there is a cycle {i1,...,iK}, then for each pair (i,j) of
   vertices in the cycle, there will be some power p=p(i,j) for which the
   matrix A^p will have a strictly positive coefficient A^p(i,j).  But then
   B(i,j)>0, since it is the sum of nonnegative numbers.

   (<==) Suppose B(i,j)>0 for all pairs (i,j) taken from {i1,...,iK}.  Then
   for each such pair (i,j) there must be some finite power p=p(i,j) such
   that A^p(i,j)>0. []

   DEFINITION.  The eventual connection matrix C of a digraph A has the
   coefficients C(i,j)=1, if j->i, with C(i,j)-0 otherwise.

   LEMMA.  Put B=A/(Id-cA) for 0<c<0.5/||A||.  Then C(i,j)=1 if and only if
   B(i,j)>0. []

   DEFINITION.  Let r=r(C) be the vector of row sums of the {0,1}-matrix C.
   Then r(C) = C*(1,...,1)'.  Let r0=r0(C) be the nondecreasing rearrangement
   of r(C).

   DEFINITION.  Let s=s(C) be the vector of column sums of the matrix C.
   Then s(c) = (1,...,1)*C.  Let s0=s0(C) be the nondecreasing rearrangement
   of s(C).


   DEFINITION.  A trivial cycle is a single vertex i for which A(i,i)=1.

   THEOREM. A contains no cycles if its eventual connection matrix C
   satisfies r0(C) < (1,2,...,N).

 */

#include <stdio.h>
/* #define N 42 */
#define N 38
typedef int Matrix[N][N];

void
multiply( Matrix out, Matrix inleft, Matrix inright )
{
  int i,j,k;
  for(i=0; i<N; i++)
    for(j=0; j<N; j++){
      out[i][j]=0;
      for(k=0; k<N; k++) out[i][j] += inleft[i][k]*inright[k][j];
    }
  return;
}

void
addto( Matrix out, Matrix in )
{
  int i,j;
  for(i=0; i<N; i++) 
    for(j=0; j<N; j++) out[i][j] += in[i][j];
  return;
}

void
nonzero( Matrix out, Matrix in )
{
  int i,j;
  for(i=0; i<N; i++)
    for(j=0; j<N; j++) out[i][j]= (in[i][j]!=0);
  return;
}

void
assign( Matrix out, Matrix in )
{
  int i,j;
  for(i=0; i<N; i++)
    for(j=0; j<N; j++) out[i][j] = in[i][j];
  return;
}

void
read(Matrix out)
{
  int i,j;
  for(i=0; i<N; i++)
    for(j=0; j<N; j++) scanf("%d",&out[i][j]);
  return;
}

void
write(char *msg,  Matrix out)
{
  int i,j;
  char label[11]="123456789 ";
  if(msg){
    puts(msg);
    puts("+-123456789 123456789 123456789 123456789 12");
  }
  for(i=0; i<N; i++){
    if(msg) {
      putchar(label[i%10]);
      putchar(' ');
    }
    for(j=0; j<N; j++) printf("%1d",out[i][j]);
    putchar('\n');
  }
  return;
}

void
dump(Matrix out)
{
  int i,j;
  for(i=0; i<N; i++){
    for(j=0; j<N; j++) printf("%1d ",out[i][j]);
    putchar('\n');
  }
  return;
}

int
main(void)
{
  int n;
  Matrix a, b, c, id;

  read(a);
  write("ORIGINAL MATRIX:",a);
  assign(b,a);

  for(n=0; n<N; n++) {
    multiply(c,b,a);
    addto(c,a);
    nonzero(b,c);
  }
  write("EVENTUAL CONNECTION MATRIX:",b);
}
