/* this program does linear algebra over mod systems */
/* it was written by a former CMSC 443 student       */
/* it could use a few more comments                  */
/* to compile this program use gcc, you do not need  */
/* to link the math library                          */
/*                                                   */
#include <stdio.h>
#include <stdlib.h>

int **cal(n)
{
int i,j;
int **c;
c=calloc(n,sizeof(int));
for (i=0;i<n;i++)
   c[i]=calloc(n,sizeof(int));
return c;
}

int *modinv(int m[],int r)
{
int i,j;
for (i=i;i<r;i++)
 m[i]=0;
for (i=1;i<r;i++){
 if (m[i]==0){
   for (j=1;j<r;j++){
     if ((i*j)%r==1){
       m[i]=j;
       m[j]=i;
       break;
     }
   }
 }
}
return m;
} 

int **mult(int **a,int **b,int r,int n)
{
int i,j,k,q=0;
int **c;
c=cal(n);
for(i=0;i<n;i++){
 for (j=0;j<n;j++){
  for (k=0;k<n;k++)
     q=q+((a[i][k])*(b[k][j]))%r;
  c[i][j]=q%r; 
  q=0;
 } 
}
return c;
} 

int **invs(int **a,int m[],int r,int n)
{
int i,j,k,q,cnt,min,minj;
int **z;
z=cal(n);
for(i=0;i<n;i++){
 for(j=0;j<n;j++){
  if (i==j)
    z[i][j]=1;
  else
    z[i][j]=0;
 }
}
for(i=0;i<n;i++){
  do{
   cnt=0;
   min=r;
   minj=-1;
   for(j=i;j<n;j++){  
    if((a[j][i]<min)&&(a[j][i]!=0)){
      min=a[j][i];
      minj=j;
    }
   }
   if(minj==-1){
     printf("%s\n","Sorry! Matrix is singular.");
     z[0][0]=-1;
     return z;
   }
   for(j=i;j<n;j++){
    if(j!=minj){
      q=(a[j][i])/min;
      for(k=0;k<n;k++){
       a[j][k]=(a[j][k]-q*a[minj][k])%r;
       z[j][k]=(z[j][k]-q*z[minj][k])%r;
       if(a[j][k]<0)
         a[j][k]=a[j][k]+r;
       if(z[j][k]<0)
         z[j][k]=z[j][k]+r;
      } 
      if(a[j][i]==0)
        cnt=cnt+1;
    }
   }
  }while(cnt<(n-i-1));
   if(minj!=i){
     for(k=0;k<n;k++){
       q=a[minj][k];
       a[minj][k]=a[i][k];
       a[i][k]=q;  
       q=z[minj][k];
       z[minj][k]=z[i][k];
       z[i][k]=q;
     }
   }
}  
for(i=0;i<n;i++){
  if (m[a[i][i]]==0){
    printf("%s\n","Sorry! Matrix is singular!");
    z[0][0]=-1;
    return z;  
  }
  else{ 
    q=m[a[i][i]]; 
    for(j=0;j<n;j++){ 
      a[i][j]=(a[i][j]*q%r);
      z[i][j]=(z[i][j]*q%r); 
    }
  }
  for(j=0;j<i;j++){  
   q=a[j][i];
   for(k=0;k<n;k++){
     a[j][k]=(a[j][k]-q*a[i][k])%r; 
     if (a[j][k]<0)
       a[j][k]=a[j][k]+r;
     z[j][k]=(z[j][k]-q*z[i][k])%r;
     if (z[j][k]<0)
       z[j][k]=z[j][k]+r; 
   }
  }
}
return z;
}
   
int **solve(int **a,int **b,int m[],int r,int n)
{
int **c,**d;
c=cal(n);
d=cal(n);
d=invs(a,m,r,n);
c=mult(d,b,r,n);
return c;
}

main()
{
int i,j,n,r,k,s;
int **a,**b,**c,*m;
char resp;
printf("%s\n","This program is designed to perform one of three functions:");
printf("%s\n","They are all matrix operations in the system");
printf("%s\n","integers modulo r (r is of your choice).");  
printf("%s\n","1. To multiply two nxn matrices (n is of your choice).");
printf("%s\n","2. To invert an nxn matrix ");
printf("%s\n","    (or alert the user that the matrix is singular(non-invertible)).");
printf("%s\n","3. To solve a matrix equation Ax=b for some nxn matrix A and");
printf("%s\n","    nx1 matrix b (provided the solution is unique, i.e. A is invertible).");
printf("%s\n","All nxn matrices must be entered row by row (reading across).");
printf("%s\n","Integers in a mod r system must be represented by their unique");
printf("%s\n","representatives in the set {0,1,..,r-1}");
printf("%s\n","(for instance, in mod 12 do not write 3 as 15,and 12 is entered as 0).");
printf("%s\n","This is precisely how values will be returned to you."); 
label:printf("%s\n","Perform function 1, 2, or 3?");
label2:scanf("%d",&s);
if (s<1||s>3){
  printf("%s\n","Come Again?");
  goto label2;
}  
printf("%s\n","Input matrix size");
scanf("%d",&n);
printf("%s\n","Input modulus");
scanf("%d",&r);
m=calloc(r,sizeof(int));
m=modinv(m,r);
a=cal(n);
b=cal(n);
c=cal(n);
if (s==1){
  printf("%s\n","Input the elements of matrix 1 (by rows)");
  for (i=0;i<n;i++){
   for (j=0;j<n;j++){
     L1:scanf("%d",&k);
     if ((k<0)||(k>=r)){
       printf("%s\n","Out of mod range. Ignored.");
       goto L1;
     }  
     a[i][j]=k;
   }
  } 
  printf("%s\n","Do the same for matrix 2");
  for (i=0;i<n;i++){
   for (j=0;j<n;j++){
     L2:scanf("%d",&k);
     if ((k<0)||(k>=r)){
        printf("%s\n","Out of mod range. Ignored.");
        goto L2;
     }
     b[i][j]=k; 
   }
  }
  c=mult(a,b,r,n);
  printf("%s\n","Product matrix is:");
  for(i=0;i<n;i++){
   for(j=0;j<n;j++)
    printf("%d%s",c[i][j]," ");  
   printf("%s\n"," ");
  }
}

if (s==2){
  printf("%s\n","Input matrix elements (by rows)");
  for (i=0;i<n;i++){
   for (j=0;j<n;j++){
     L3:scanf("%d",&k);
     if ((k<0)||(k>=r)){
       printf("%s\n","Out of mod range. Ignored."); 
       goto L3;
     }
     a[i][j]=k;  
   }
  }
  c=invs(a,m,r,n);
  if(c[0][0]!=-1){
  printf("%s\n","Inverse matrix is:");
    for(i=0;i<n;i++){
     for(j=0;j<n;j++)
      printf("%d%s",c[i][j]," ");
     printf("%s\n"," ");
    } 
  }
}

if (s==3){
  printf("%s\n","Input matrix A's elements (by rows)");
  for (i=0;i<n;i++){
   for (j=0;j<n;j++){
     L4:scanf("%d",&k);
     if ((k<0)||(k>=r)){
       printf("%s\n","Out of mod range. Ignored.");
       goto L4;
     }
     a[i][j]=k;
   }
  }	
  printf("%s\n","Input matrix b's elements");
  for (i=0;i<n;i++){
     L5:scanf("%d",&k);
     if ((k<0)||(k>=r)){
       printf("%s\n","Out of mod range. Ignored.");
       goto L5;
     }
     b[i][0]=k;
  }
  c=solve(a,b,m,r,n);
  printf("%s\n","Solution vector is:");
  for(i=0;i<n;i++){
  printf("%s%d%s%d\n","x",i+1," = ",c[i][0]);  
  }
}
printf("%s\n","Run Again? (type y or n)");
scanf("%s",&resp);
if (resp=='y')
  goto label; 
}