/* MatInverse.c */

#include <stdio.h>

#define DataFile "MatrixDataFile2.txt"
#define NN 1001

    /* maximum number entries = NN-1 */
#define INF 0XEFFFFF

void KBToMat(int*);
int FileToMat(int*);
void DisplayMat(int*);
void MatAtoFile(int*);
int Det(int*);
int Adjoint(int*,int*);
int MatProd(int*,int*,int*);

void main(void)
{
 int   K,Mat[NN],Adj[NN];   char CH;

 puts("Where is the data for the matrix?");
 puts("If it will be typed in on the (K)eyboard, type K");
 puts("If it is already in the GraphData (F)ile, type F");
INX:
 CH=(char)getchar(); getchar();
 if ((CH=='K') || (CH=='k'))
  {puts("Input square matrix, 2x2 or larger"); KBToMat(Mat); goto NXT;}
 if ((CH=='F') || (CH=='f')) {if (FileToMat(Mat)) goto NXT; else goto FIN;}
 puts("Error: enter either K or F");  goto INX;

NXT:
 puts("\nThe given matrix is:");
 DisplayMat(Mat);
 puts("");
 if (Mat[0]/100 != Mat[0]%100)
  {puts("ERROR: matrix is not square; there is no inverse matrix"); goto FIN;}
 K=Det(Mat);
 if (!K) {puts("Matrix is singular; there is no inverse matrix"); goto FIN;}

 Adjoint(Mat,Adj);
 puts("For the given matrix");
 printf("  its determinant value is %d\n",K);
 puts("  and its adjoint matrix is:");
 DisplayMat(Adj);
 puts("");
 if (K==1) {puts("The adjoint matrix is also the inverse matrix"); goto FIN;}
 printf("Divide the determinant value %d into the adjoint matrix\n",K);
 puts("  to get the inverse matrix");
 
FIN:
 getchar();
}

/*****************************************************************************/

int FileToMat(int *P)
{
 FILE *F;   char TITLF[]=DataFile;
 char CH,*S,GrStr[NN];  int I,K,M,N=0,CNT,Sign,   *PP;

  if ((F=fopen(TITLF,"r"))==NULL)
    {printf("Cannot find file %s\n",TITLF); return 0;}

     /************* list input-datas already in file ***********/

CH=(char)fgetc(F);     PP=P;

LP2:
 while ((!feof(F)) && (CH != '[')) CH=(char)fgetc(F);     /* find next [ */

 K=0;
 if (!feof(F))
  {
   N++; printf("%d  ",N);
   while (CH != ']') {putchar(CH); CH=(char)fgetc(F);}
   puts("]\n");
   goto LP2;
  }

 INX:
  puts("Input number before the selected data-input (zero to abort):");
  scanf("%d",&K); getchar();
  if (!K) {fclose(F); return 0;}
  if ((K<0) || (K>N)) {printf("Number must be between 1 and %d\n",N); goto INX;}

  rewind(F);

  CH = (char)fgetc(F);
  for (I=1; I<K; I++)
   {
    while (CH != ']') CH = (char)fgetc(F);
   CH = (char)fgetc(F);
  }

  while (CH != '[') CH = (char)fgetc(F);
  S=GrStr;    *S=CH;         /* CH =  [ */
  S++;

  CH = (char)fgetc(F); if ((CH=='m') || (CH=='M')) {*S='M'; S++;}
  else
   {
    puts("Input-array in file not marked 'M' (after '[') for Matrix data");
   return 0;
  }

    /**** Copy character array in file to GrStr ********/

  CH = (char)fgetc(F);
  while (CH!=']') {*S=CH; S++; CH = (char)fgetc(F);}
  *S = CH;
  fclose(F);
  S++; *S='\0';
/*
  puts("Input-data from file is:");
  puts(GrStr);   /* getchar();   /* check data file */

           /*********GrStr to A ******************/

  S=GrStr;


 S++;

 /* Get Size of matrix */
 while ((*S<'0') || (*S>'9')) S++;
 K=0; while ((*S>='0') && (*S<='9')) {K = 10*K + *S-48; S++;}
 M=K/100;  N=K%100;


 *P=100*M + N; P++;

 CNT=M*N;       /* CNT = total number of entries */

 while (CNT)
  {
   while ((*S<'0') || (*S>'9')) S++;
   S--; if (*S=='-') Sign=-1; else Sign=1; S++;
   K=0; while ((*S>='0') && (*S<='9')) {K = 10*K + *S-48; S++;}
   if (Sign==1)*P=K; else *P=-K;  P++;
   CNT--;
  }
 *P=INF;
 return *PP;
}

/*****************************************************************************/

void KBToMat(int *P)
{
 int I,J,K,M,N;


INX:
 puts("Enter size of matrix ROWSxCOLUMNS (0x0 termimates input)  ");
 scanf("%dx%d",&M,&N); getchar();
 if ((M>12) || (N>12)) {puts("No size larger than 12"); goto INX;}
 *P = 100*M + N;  P++;

 for (I=1; I<=M; I++)
  {
   printf("Prepare to enter %d numbers in row %d\n",N,I);
   for (J=1; J<=N; J++)
    {
     if (J==1) puts("first number");
     if ((J>1) && (J<N)) puts("next number");
     if (J==N) puts("last number");
     scanf("%d",&K); getchar();
     *P=K; P++;
    }
  }

 *P=INF;
 return;
}

/*************************** Print A as a Matrix *****************************/

void DisplayMat(int *P)        /* P -> linear matrix array */
 {
  int I,J,K,M,N,Max=0,T;    char S1[] = "%1d%2d%3d%4d%5d%6d%7d%8d%9d",S2[3];

  M=*P/100;  N=*P%100;  P++;
  for (I=1; I<=M; I++)
   {
    for (J=1; J<=N; J++)
     {
      if (*P>=0) K=*P; else K=-*P;
      if (K>Max) Max=K;
    }
  }

 T=1; K=1;
 while (T<=Max) {T=T*10; K++;} 
 for (I=0; I<3; I++) S2[I] = S1[3*K+I];  S2[3]='\0';


 for (I=1; I<=M; I++)
   {
    for (J=1; J<=N; J++) {printf(S2,*P); P++;}
    puts("");
  }

 return;
}

/**************************************************************************/

int Det(int *A)
   /*  NN has been defined globally so that NN-1 is the maximum length of
         any permutation (to prevent rollover of counters).
       Receive integer N, 2<=N<NN, and a pointer A to an array of length NN+1;
       The first permutation is odd; after that the odd-even alternate */
{
 int CNT,H,I,J,K,N,PROD,COPY;    int ODD=1,SUM=0;
 int R[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0},
     D[]={1,1,1,1,1,1,1,1,1,1,1,1,1,1},
  PERM[]={0,1,2,3,4,5,6,7,8,9,10,11,12,13};

 N=A[0]%100;

 if (N==1) return A[1];

 CNT=1; for (I=1; I<=N; I++) CNT*=I;   /* CNT = N! */

 while (CNT)
  {              /* form a sum of products of N! terms from matrix A */
   K=0; I=N;
  INDEX:      /* but first, construct a new PERMutation from the old one */
     R[I] = J = R[I] + D[I];
     if (J==I) {D[I]=-1; goto LOOP;}
     if (J) goto TRANSPOSE;
     D[I]=1; K++;
  LOOP:
     if (I>2) {I--; goto INDEX;}
     J=1;
  TRANSPOSE:
     J+=K; COPY=PERM[J]; PERM[J]=PERM[J+1]; PERM[J+1]=COPY;
     /* note: odd and even permutations alternate during the generation */

   /* Using permutation PERM form a product PROD of entries from the matrix A */
     PROD=1;
     for (H=1; H<=N; H++) PROD*= A[(H-1)*N + PERM[H]];   /* row H */
    /* PROD = A[1,PERM[I]]* **** * A[N,PERM[N]] */

/*If permutation even, PROD is added to SUM; if odd then -PROD is added */
     if(ODD) SUM-=PROD; else SUM+=PROD;
     ODD=1-ODD;
     CNT--;
  } /* end of while CNT */   
 return SUM;
}

/***************************************************************************/

int Adjoint(int *Mat, int *Adj)
{
 int I,II,JJ,J,M,N,      *P,*Q,*R,                Buf[200], Cof[NN];

 M = Mat[0]/100;  N = Mat[0]%100;
 if (M != N) {puts("Error: Matrix not square"); return INF;}

      Q=Cof;  *Q=Mat[0];       Buf[0] = (M-1)*100 +N-1;
 for (II=1; II<=M; II++)
  for (JJ=1; JJ<=N; JJ++)
   {
     P=Mat; R=Buf;
     for (I=1; I<=M; I++)
    for (J=1; J<=N; J++)
     {
      P++;
      if ((I != II) && (J != JJ)) {R++; *R = *P;}
     }

    Q++;  if((II+JJ)%2) *Q = -Det(Buf); else *Q=Det(Buf);
   } /* end II,JJ */

/* DisplayMat(Cof); */

  P=Cof; P++;   Adj[0] = M*100 + N;

 for (J=1; J<=N; J++)
  {
   Q=&Adj[J];
   for (I=1; I<=M; I++) {*Q=*P; P++; Q = Q + N;}
  }
 return Adj[0];
}

/********************** MatProd C = Mat A x Mat B *************************/

int MatProd(int *A, int *B, int *C)
{
 int I,J,K,M1,N1,M2,N2,Sum,   *PP, *P, *Q, *R;  

 M1 = A[0]/100; N1 = A[0]%100;  M2 = B[0]/100; N2 = B[0]%100;

 if (N1 != M2)
  {
   puts("ERROR:");
   puts("First matrix not conformable to second matrix for multiplication");
   return 0;
  }

  PP= &A[1];    C[0] = 100*M1 + N2;  R = &C[1];

 for (I=1; I<=M1; I++)
 {
  for (J=1; J<=N2; J++)
   {
    P=PP; Q=&B[J]; Sum=0;
    for (K=1; K<=N1; K++)
     {
      Sum = Sum + *P**Q;  /*  printf("I=%d J=%d *P=%d *Q=%d\n",I,J,*P,*Q); */
      P++; Q = Q + N2;
     }
   *R = Sum;  R++;          /*    printf("I=%d J=%d Sum = %d\n\n",I,J,Sum); */
   } /* end J */

  PP = PP + N1;

 } /* end I */
 *R = INF;

 return C[0];
}

/*****************************************************************************/