/* TestCofactor.c */

#include <stdio.h>
#define INF 0XEFFFFF
#define NN 300

int Adjoint(int*,int*);
void DisplayMat(int*);
int Det(int*);

void main(void)
{
 int Adj[NN],
  Mat[] = {303,  1,2,3,   0,4,5,  1,0,6};
   /*{404, 1,2,3,4,  5,6,7,8,  9,10,11,12,  13,14,15,16,INF}; */


 if (!Adjoint(Mat,Adj)) goto FIN;

 puts("For the matrix:");
 DisplayMat(Mat);
 printf("its determinant value is %d\n",Det(Mat));
 puts("and its adjoint is:");
 DisplayMat(Adj);



FIN:
 getchar();
}

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

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];
}

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

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)          /* A -> linear matrix file */
   /*  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<2 || N>NN)
  {printf("Input %d  not between limits  2 and %d\n",N,NN); return 0;}
 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;
}