// class  ECHELON will generate matrix with row-reduced-echelon form

#include <iostream.h>
#include <math.h>

const int MAX_ROW = 58;
const int MAX_COL = 65;

class ECHELON64

{
 private:
  int msgm[MAX_ROW][MAX_COL];  // initial matrix
  int kk;    // matrix row.
  int nn;    // matrix column.
  int Right[MAX_COL];  // Right is the index of the rightmost nonzero entry.
  int Left[MAX_COL];	   // Left is the index of the leftmost nonzero entry.
  void toMSGM();
  int echelon[MAX_ROW][MAX_COL];

 public:
  ECHELON64(int[][65], int, int );
  void Find_L_R(int[], int, int *, int *);
  int* get_Left();
  bool IsMatch(int * , int * );
  void display_echelon();
  void Modify();// to sort the rows that can be row reduced form
  void set_echelon();
  int* get_echelon();
  void display_MSGM();//for obtain the state complexity profile
  int* statecomplex();
 
};

ECHELON64::ECHELON64(int G_Matrix[][65], int _k, int _n)
{
  kk = _k;    nn = _n; 
  for(int i = 0; i < kk; i++) {
    for(int j = 0; j < nn; j++) {
      msgm[i][j] = G_Matrix[i][j];
    echelon[i][j] = 0;
    }
  }
  for(int ii = 0; ii < kk; ii++)
    {
      Left[ii] = 0; Right[ii] = 0;
    }
}

void ECHELON64::Find_L_R(int G[], int row, int *L, int *R) //to obtain the leftmost index and rightmost index
{
  for(int i = 0; i < nn; ++i)
    {
      if (G[i] == 1) {
	L[row] = i;
	break;
      }
    }

  for(int j = nn - 1; j >= 0; --j)
    {
      if (G[j] == 1)
	{
	  R[row] = j;
	  break;
	}
    }
}

bool ECHELON64::IsMatch(int *L, int *R)
{
  for(int i = 0; i < kk; i++)
    {
      for(int j = i + 1; j < kk; j++)
	{
	  if((L[i] == L[j])  || (R[i] == R[j]) )
	    {
	      return true;
	      break;
	    }
	}
    }
  return false;
}

	
void ECHELON64::toMSGM() // if it is successful to find a pair of Gi and Gj such that satisfy the above the function, then do either Gi=Gi+Gj or Gj=Gi+Gj

{
  int i, j;

  for(i = 0; i < kk; i++)
    Find_L_R(msgm[i], i, Left, Right);
  while(IsMatch(Left, Right))
    
    {
      for ( i = 0; i < kk; i++)
	{ 
	  for (j = i + 1; j < kk; j++)
	    {
	      if ( Left[i] == Left[j] )
		{
		  if ( Right[i] >= Right[j] )
		    {
		      for(int n1 = 0; n1 < nn; n1++)
			msgm[i][n1] = msgm[i][n1] ^ msgm[j][n1];
		      Find_L_R(msgm[i], i, Left, Right);
		    }
		  else
		    {
		      for(int n1 = 0; n1 < nn; n1++)
			msgm[j][n1] = msgm[i][n1] ^ msgm[j][n1];
		      Find_L_R(msgm[j], j, Left, Right);
		    }	
		}	
	      else
		{
		  if(Right[i] == Right[j]){
		    if ( Left[i] >= Left[j] )
		      {	
			for(int n2 = 0; n2 < nn; n2++)
			  msgm[j][n2] = msgm[i][n2] ^ msgm[j][n2];
			Find_L_R(msgm[j], j, Left, Right);
		      }
		    else
		      {
			for(int n2 = 0; n2 < nn; n2++)
			  msgm[i][n2] = msgm[i][n2] ^ msgm[j][n2];
			Find_L_R(msgm[i], i, Left, Right);
		      }
		    
		  }
		}
	    }
	}
    }
}
void ECHELON64::Modify() // sort the MSGM by the order of Left(row), which is the index of the leftmost nonzero entry of Gi
{
  int numk=0, hold[58][65];
  int i, j;
  
  toMSGM();
  
  for(i=0; i < nn; i++)
    {
      for( j=0; j < kk; j++)
	{
	  if(Left[j] == i)
	    {
	      for(int m=0; m < nn; m++)
		hold[numk][m] = msgm[j][m];
	      numk++;
	      break;
	    }
	}
    }
  for(i=0; i < kk; i++)
    for(j=0; j < nn; j++)
      msgm[i][j] = hold[i][j];
}

void ECHELON64:: set_echelon() //modify the sorted Generated Matrix ( as leftmost indices:refer to the result in the end) to the row-reduced-echelon (systematic-like Generation Matrix), then the trellis can be  drawed followed by this row-reduced-echelon Matrix( systematic like G)
{ 
  int i,j;
  Modify();
  for(i=0; i < kk; i++)
    Find_L_R(msgm[i], i, Left, Right);
  
  for(i=0; i < kk; i++)
    {
      for(j=i-1; j >= 0; j--)
	{
	  if(msgm[j][Left[i]] == 1)
	    for ( int n1=0; n1 < nn; n1++)
	      msgm[j][n1] = msgm[i][n1] ^ msgm[j][n1];
	  
	}
    }
  for(i = 0; i < kk; i++)
    for(j = 0; j < nn; j++)
      echelon[i][j] = msgm[i][j];
}
int* ECHELON64::get_echelon()
{
  set_echelon();
  return &echelon[0][0];
}

int* ECHELON64::get_Left()
{  
  for(int i = 0; i < kk; i++)
    {
      for(int j = 0; j < nn; j++)
	{
	  if( echelon[i][j] == 1 )
	    {
	      Left[i] = j;
	      break;
	    }
	}
    }
  return Left;
}

void ECHELON64::display_echelon() 
{   
  set_echelon();
  for(int i = 0; i < kk; i++)
   {
      for(int j = 0; j < nn; j++)
  	cout << echelon[i][j] << " ";
      cout << "\n";
   }
}
void ECHELON64::display_MSGM()//for obtaining state complexity profile

{
  Modify();
  /*  for(int i = 0; i < kk; i++) */
  /*      { */
  /*        for(int j = 0; j < nn; j++) */
  /*  	cout << msgm[i][j] << " "; */
  /*        cout << "\n"; */
  /*      } */
}
int* ECHELON64::statecomplex()
{
  int statecomplex_file[65];
  int i, j, k;
  int past[65], future[65];
  int node_number = 0; 
  
  for( j = 0; j < kk; j++)
    Find_L_R(msgm[j], j, Left, Right);
  
  for ( i = 0; i < nn; i++)
    {
      past[i] = 0;
      future[i] = 0;
    }
  for ( i = 0; i < nn; i++)
    {
      for ( j = 0; j < kk; j++)
	{
	  if ( Right[j] <= i)
	    past[i]++;
	}
      for ( k = 0; k < kk; k++)
	{
	  if (Left[k] >= i+1)
	    future[i]++;
	}
    }
  statecomplex_file[0] = 0;
  for (i = 1; i <= nn; i++)
    {
      statecomplex_file[i] = kk - past[i-1]-future[i-1];
    }
 
  return &statecomplex_file[0];
}