/****************************************************************************
 * This program does linear algebra with modular arithmetic.	
 * It was written by a former CMSC 443 student. 			  
 *  			  
 *  			 Course : CMSC 443
 *			  
 * 	     	    
 *****************************************************************************
   This project was designed to perform some matrix operations for the
 modulo systems. It reads  and saves data from data files.
 However, it was designed to use the command line to
 read and save data.
 this program can perform  matrix multiplication, calculation of inverse, and
 solving linear systems.

 The input file should be structured as:
   size of matrix on first line
   followed by n by n + 1 matrix where n is size of the matrix
      and the last column represents the right hand side of a
      linear system that is to be solved. If the matrix is to
      be inverted and there is no system to be solved then create
      a dummy column.
   the last line should contain the modulus.

 Compile this program with 'gcc' or 'cc'.
 At the command line use 'a.out nameofinputfile nameofoutputfile'
  (without quotes) 
 *****************************************************************************/

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

#define  size  50

void  SaveResult( int [size][size], int );
FILE  *InFile,  *OutFile;

int main( int argc, char *argu[] )
{
   int   M[size][size];		   /* main matrix */
   int   Product_M[size][size];    /* store the product of two matries */
   int   Inv_M[size][size];        /* store the inverse of main matrix */
   int   flag=1;
   int   M_size,		   /* size of main matrix */
         modulus;		   /* modulus system */
   int   selection;

   void  Initial_Matrix( int [size][size], int, FILE * );
   int   Determinant( int [size][size], int, int );
   void  Inverse( int [size][size], int [size][size], int, int );
   void  Multiplication( int [size][size], int [size][size],
				 int, int, FILE * );
   void  Solve_Equation( int [size][size], int, int );
   
   if (argc != 3)
   {
      printf("\n*****  Usage : %s infile outfile  *****\n", argu[0]);
      printf("You forgot to indicate the input and output files.\n");
      exit(1);
   }
   else if ((InFile = fopen(argu[1], "r")) == NULL)
   {
      printf("Couldn't open %s file.\n", argu[1]);
      exit(1);
   }
   else if ((OutFile = fopen(argu[2], "w")) == NULL)
   {
      printf("Couldn't open %s file.\n", argu[2]);
      exit(1);
   }

   fscanf( InFile, "%d", &M_size );
   if( M_size < 0 )
      printf("Error input for matrix with size 0\n");
   else
   {
      Initial_Matrix( M, M_size, InFile );    /* read in data for main matrix */
      fscanf( InFile, "%d", &modulus );       /* read in modulus limit */
      printf("\tModulus = %d\n", modulus);
      while(flag)
      {
         printf("\n\t***********************************\n");
         printf("\t*                                 *\n");
         printf("\t*       (1) Multiplication        *\n");
         printf("\t*       (2) Inverse matrix        *\n");
         printf("\t*       (3) Solve equcations      *\n");
         printf("\t*       (4) Exit                  *\n");
         printf("\t*                                 *\n");     
         printf("\t***********************************\n");
         printf("\n\n--> Please select one operation < ");
         while( (selection = getc( stdin )) < '1' || selection > '4' ) 
         {
            if( selection == '\n' );
            else
            {
               printf("\n** Your selection was not accepted.\n");
               printf("--> Please enter your selection again between 1 to 4 < ");
	    }
         }

         switch(selection)
         {
             case '1' :   /* Multiplication */
                   Multiplication( M, Product_M, M_size, modulus, InFile );
                   break;
             case '2' :   /* Inverse matrix */
              /* if determinant is 0, it is a singlar matrix, can't do inverse */
                   if( Determinant( M, M_size, modulus ) )
                       Inverse( M, Inv_M, M_size, modulus );
                   else
                       printf("\n\tSorry!! The input matrix is a singlar matrix\n");
                   break;
             case '3' :   /* Solve equcations */
                   Solve_Equation( M, M_size, modulus );
                   break;
             case '4' :   /* Exit */
                   printf("\n\t********** Thank you! Bye! **********\n");
                   exit(1);
             default:
                   break;
         }
      }
   }
   return 0;
}

/****************************  INITIAL_MATRIX  *****************************
   This function is used to read in matrix data from input file. The main
 matrix is called M. I passed the size of matrix, main matrix, and file
 pointer as its arguments.
 ***************************************************************************/
void  Initial_Matrix( int M[size][size], int M_size, FILE *InFile )
{
   int  i, j, loop=0;

   for( i=0 ; i<M_size ; i++ )                  
      for( j=0 ; j<M_size+1 ; j++ )
         fscanf(InFile, "%d", &M[i][j]);   /* read in matrix element */

   printf("\n***** The input %d*%d matrix is *****\n", M_size, M_size);
   for( i=0 ; i<M_size ; i++ )
   {
      printf("\t(");
      for( j=0 ; j<M_size+1 ; j++ )
      {
         if( j == M_size )
             printf(" | ");
         printf(" %3d ", M[i][j]);   /* display man matrix */
      }
      printf(")\n");
   }
} 

/****************************  MULTIPLICATION  *****************************
    This function tries to do imatrix mul lication.  It takes main matrix,
 product matrix, file pointer, matrix size, and modulus as its parameters.
 It will read in the multiplicand matrix from input file, and it processes
 the multiplication by row timing column.
 ***************************************************************************/
void  Multiplication ( int M[size][size], int Product_M[size][size],
		          int M_size, int modulus, FILE *InFile )
{
   int  i, j, k;
   int  Mult[size][size];   /* store the multiplicand matrix */
   int  temp;               /* store temporary sum of products */
   int  option;             /* option to save the final product */

   printf("\n\t*****************************************\n");
   printf("\t*            Multiplication             *\n");
   printf("\t*****************************************\n");

   for( i=0 ; i<M_size ; i++ )
      for( j=0 ; j<M_size ; j++ )
	 fscanf( InFile, "%d", &Mult[i][j]);  /* read in multiplicand */

   printf("\n*****  The multiplicand matrix is  *****\n");
   for( i=0 ; i<M_size ; i++ )
   {
      printf("\t(");
      for( j=0 ; j<M_size ; j++ )
         printf(" %3d ", Mult[i][j]);   /* display multiplicand */
      printf(")\n");
   }

   for( i=0 ; i<M_size ; i++ )
      for( j=0 ; j<M_size ; j++ )
      {
	 temp = 0;
         for( k=0 ; k<M_size ; k++ )
	    temp += M[i][k] * Mult[k][j];   /* sum each sub_products */
         Product_M[i][j] = temp % modulus;  /* put result into Product matrix*/
      }

   printf("\n*****  The product matrix is  *****\n");
   for( i=0 ; i<M_size ; i++ )
   {
      printf("\t(");
      for( j=0 ; j<M_size ; j++ )
         printf(" %3d ", Product_M[i][j]);   /* display product matrix */
      printf(")\n");
   }

   printf("Do you want to save the product matrix? ('y' or 'n') < ");
   while( 1 ) 
   {
       option=getc(stdin);
       if( option == '\n' );
       else if( option == 'y' || option == 'n' )    break;
       else 
            printf("\n*** Please enter 'y' to save otherwise enter 'n' ***");
   }

   if( option == 'y' )
   {
       fprintf(OutFile, "\n*****  The product matrix is  *****\n");
       SaveResult( Product_M, M_size );
   }
}

/*****************************  DETERMINANT  *******************************
   This function calculates the determinant of the input matrix. It takes the
 main matrix , matrix size, and the modulus as auguments, and return the
 absolute value of determinant.
 ***************************************************************************/
int  Determinant( int M[size][size], int M_size, int modulus )
{
   int  i, j, p=0;
   int  temp,        /* temporary sub_product */
        positive=0,  /* sum of positive sub_products */
 	negative=0;  /* sum of negative sub_products */
      
   for( i=0 ; i<M_size ; i++, p++ )  /* calculate positive part */
   {
      temp = 1;
      for( j=0 ; j<M_size ; j++, p++)
	  temp *= M[j][ p%M_size ];   /* calculate sub_product */
      positive += temp;               /* sum all positive part of sub_product */
   }

   for( i=0 ; i<M_size ; i++, p++ )  /* calculate negative part */
   {
      temp = 1;
      for( j=(M_size-1) ; j>=0 ; j--, p++)
	  temp *= M[j][ p%M_size ];  /* calculate sub_product */
      negative += temp;              /* sum all negative part of sub_product */
   }

   return  abs( positive - negative );  /* return absolute value of the
                                           determine */
}

/*********************************  GCD  ***********************************
   This function return the GCD of two numbers. It implements the Euclid's
 algorithm by recursion.
 ***************************************************************************/
int GCD(int x, int y)
{
   int  P=0, z=0, G1=x, G2=y;

   if(x < y)
   {
      G1 = y;  
      G2 = x;
   }
   while(G2)
   {
      z = G1 % G2;
      G1 = G2;
      G2 = z;
      GCD(G1, G2);
   }
   return G1;
}

/*****************************  Find_INVERSE  ******************************
   This function can find the inverse of the input number in mod system. It
 implements the Euclid's algorithm. It returns the inverse number of input
 data in modulus system.
 ***************************************************************************/
int  Find_Inverse( int x , int modulus )
{
   int  G1=modulus, G2=x, G3;
   int  U1=1, U2=0, U3;
   int  V1=0, V2=1, V3;
   int  quotient, inver;

   while( G2 )
   {
	quotient = G1 / G2;
        G3 = G1 - G2 * quotient;
   	U3 = U1 - quotient * U2;
   	V3 = V1 - quotient * V2;
	G1 = G2;
 	G2 = G3;
	U1 = U2;
 	U2 = U3;
	V1 = V2;
	V2 = V3;
   }

   inver = V1;
   if( inver >= 0 )     return inver;
   else                 return (inver += modulus);
} 

/******************************  INVERSE  **********************************
   This function finds the inverse of the input matrix. It takes the main
 matrix, modulus, and size as auguments, and the result store in inverse
 matrix.
 ***************************************************************************/
void Inverse(int M[size][size], int Inv_M[size][size], int M_size, int modulus)
{
   int i, j, k, inver,
       p, a, b;
   int option;               /* option to save the inversed matrix */
   int temp_M[size][2*size]; /* temporary matrix to operate before save */
   int Find_Inverse(int, int);
   int GCD( int, int );

   printf("\n\t*****************************************\n");
   printf("\t*            Inverse Matrix             *\n");
   printf("\t*****************************************\n");

   for( i=0 ; i<M_size ; i++ )
      for( j=0 ; j<2*M_size ; j++ )
      {
         if( j < M_size )
	    temp_M[i][j] = M[i][j]; /* assign element to the temporary matrix */
	 else                       /* make temp_M to be (size * 2size) */
            if( (j - M_size) == i )    temp_M[i][j] = 1;
	    else                       temp_M[i][j] = 0;
      }

   for( i=0 ; i<M_size ; i++ )      /* process inverse */
   {
      for( j=i ; j<2*M_size ; j++ )
      {
	  if( i == j )
          {
             if( temp_M[i][j] < 0 )
               for( b=i ; b<2*M_size ; b++ )
  		    temp_M[i][b] = -temp_M[i][b];

      /* if GCD of the element and the modulus is 1, it will have inverse
         else there is no inverse */
             if( GCD( temp_M[i][j], modulus) == 1 )
             {
	        inver = Find_Inverse( temp_M[i][j], modulus );
	        temp_M[i][j] = 1;
	     }
             else
             {
                printf("temp_M=%d, modulus=%d\n", temp_M[i][j], modulus);
                printf("\tCAN NOT BE INVERSE IN MODULUS %d\n", modulus);
        	return;
	     }
          }
          else
	      temp_M[i][j] = (inver * temp_M[i][j]) % modulus;
      }

      for( k=0 ; k<M_size ; k++ )
      {
         a = 0;
         if( k != i )
  	 {
            p = temp_M[k][i];
            for( j=i ; j<2*M_size ; j++ )
 	       temp_M[k][j] = (temp_M[k][j] - p*temp_M[i][j]) % modulus;
	 }
      }
   }

   printf("\n***** Inverse of your input matrix is *****\n");
   for( i=0 ; i<M_size ; i++ )
   {
      printf("\t(");
      for( j=M_size ; j<2*M_size ; j++ )
      {
         if( temp_M[i][j] < 0 )
	     Inv_M[i][j-M_size] = temp_M[i][j] + modulus;
 	 else   
	     Inv_M[i][j-M_size] = temp_M[i][j];
         printf(" %3d ", Inv_M[i][j-M_size]);
      }
      printf(")\n");
   }

   printf("Do you want to save the inversed matrix? ('y' or 'n') < ");
   while( 1 ) 
   {
       option=getc(stdin);
       if( option == '\n' );
       else if( option == 'y' || option == 'n' )    break;
       else 
            printf("\n*** Please enter 'y' to save otherwise enter 'n' ***");
   }

   if( option == 'y' )
   {
       fprintf(OutFile, "\n*****  The Inverse matrix is  *****\n");
       SaveResult( Inv_M, M_size );
   }
}

/*****************************  SOLVE_EQUATION  ****************************
   This function implemets the Gaussian Elimination and the inverse function
 to solve equatons.
 ***************************************************************************/
void  Solve_Equation( int M[size][size], int M_size, int modulus )
{
   int  i, j, k, a, b, p;
   int  temp_M[size][size+1]; /* temporary processing matrix */
   int  inverse;
   char sign;

   printf("\n\t*****************************************\n");
   printf("\t*            Solve Equations             *\n");
   printf("\t*****************************************\n");

   for( i=0 ; i<M_size ; i++ )
      for( j=0 ; j<M_size+1 ; j++ )
         temp_M[i][j] = M[i][j];  /* assign element to temporary matrix */

   for( i=0 ; i<M_size ; i++ )
   {
      for( j=i ; j<M_size+1 ; j++ )
      {
	  if( i == j )
          {
             if( temp_M[i][j] < 0 )
               for( b=i ; b<M_size+1 ; b++ )
  		    temp_M[i][b] = -temp_M[i][b];

             if( GCD( temp_M[i][j], modulus) == 1 )
             {
	        inverse = Find_Inverse( temp_M[i][j], modulus );
	        temp_M[i][j] = 1;
	     }
             else
             {
                printf("\tCAN NOT BE INVERSE IN MODULUS %d\n", modulus);
        	exit(1);
	     }
          }
          else
	     temp_M[i][j] = (inverse * temp_M[i][j]) % modulus;
      }

      for( k=0 ; k<M_size ; k++ )
      {
         if( k != i )
  	 {
            p = temp_M[k][i];
            for( j=i ; j<M_size+1 ; j++ )
	       temp_M[k][j]=((-p)*temp_M[i][j]+temp_M[k][j])%modulus;
	 }
      }
   }

   printf("\n***** Solution of your input matrix is *****\n");
   for( i=0 ; i<M_size ; i++ )
   {
      printf("\tX%d = ", i+1);
      for( j=M_size ; j<M_size+1 ; j++ )
      {
         if( temp_M[i][j] < 0 )
	     temp_M[i][j] = temp_M[i][j] + modulus;
 	 else   
	     temp_M[i][j] = temp_M[i][j];
         printf(" %3d ", temp_M[i][j]);
      }
      printf("\n");
   }
}

/*******************************  SAVERESULT  ******************************
   This function just write the result matrix into output file. 
 ***************************************************************************/
void  SaveResult( int temp[size][size], int M_size )
{
   int  i, j;

   for( i=0 ; i<M_size ; i++ )
   {
      fprintf(OutFile, "\t(");
      for( j=0 ; j<M_size ; j++ )
         fprintf(OutFile, " %3d ", temp[i][j]);
      fprintf(OutFile, ")\n");
   }
}