/* TestADJ.c */
/* Testing subprogram ADJ, also receive determinant of calling matrix MAT */

#include <stdio.h>

int C[14];
int D[14];

void ADJ(int,int*,int*);
void ALTPERM(int*,int);
int DET(int,int*);

void main(void)
{
 int
  MAT[] =
/*   {3,3,
   6,1,-5, -2,-5,4, -3,3,-1,
     0,0, 0,0, 0,0};
*/

/*   {3,3,
    1,2,3, 1,3,4, 1,4,3,
     0,0, 0,0, 0,0};
*/

/*     {4,4,
   5,0,0,2, 1,1,0,2, 0,0,2,1, 1,0,0,1,
     0,0, 0,0, 0,0};
*/

       {4,4,
    1,2,3,4,  2,0,5,6,  4,8,0,2,  1,4,2,6,   
      0,0, 0,0, 0,0};

 int I,J,N,*P,ADJMAT[200];

 ADJ(MAT[0],&MAT[2],ADJMAT);

 printf("(%d,%d,\n",ADJMAT[0],ADJMAT[1]);
 N=ADJMAT[0]; P = &ADJMAT[2];
 for (I=0; I<N; I++) for (J=0; J<N; J++) {printf("%d,",*P); P++;}
 puts("");
 for (I=1; I<6; I++) {printf("%d,",*P); P++;}  printf("%d)",*P);
 getchar();

}

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

void ADJ(int N, int *PP, int *QQ ) /* N = number of rows, PP -> begin data
    QQ -> ADJMAT in calling program to contain the adjoint of MAT
      in complete array form */
{
 int I,II,J,JJ,N1,*P,*Q,*PPP;
 int SubMAT[200],MATT[200];  /* only data in these arrays;
      MATT to contain the transpose of data in MAT; */
 Q=MATT;  PPP=PP;

 for (J=0; J<N; J++)
  {
   P=PP;
   for (I=0; I<N; I++)
    {
     *Q = *P;
     P = P + N; Q++;
    }
   PP++;
  }         /* PP no longer points to data in MAT */

 *QQ = N; QQ++; *QQ = N;
   QQ++;  /* QQ -> part of array in calling program to receive adjkoint data */  
 
 PP=MATT; N1 = N-1;

 for (II=0; II<N; II++)
 for (JJ=0; JJ<N; JJ++)
  {
   Q=SubMAT;
   for (I=0; I<II; I++)
    {
     P = PP + I*N;
     for (J=0; J<JJ; J++) {*Q = *P; P++; Q++;}
     P++;
     for (J=JJ+1; J<N; J++) {*Q = *P; P++; Q++;}
    }
   for (I=II+1; I<N; I++)
    {
     P = PP + I*N;
     for (J=0; J<JJ; J++) {*Q = *P; P++; Q++;}
     P++;
     for (J=JJ+1; J<N; J++) {*Q = *P; P++; Q++;}
    }


   if ((II+JJ)&1) *QQ = -DET(N1,SubMAT); else *QQ = DET(N1,SubMAT);
   QQ++;

 } /* end II,JJ */

 *QQ = 1;   /* flag to enliven determinant */
 QQ++; *QQ = DET(N,PPP); QQ++;
 for (I=1; I<5; I++) {*QQ = 0; QQ++;}
 


}

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

int DET(int N,int *Q)
/* receive number N of rows in a square matrix, Q -> begin data;
    return determinant value of matrix */
{
 int I,J,K,*P,NFACT=1,SUM=0,PROD,
  PERM[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13};

 for (I=0; I<14; I++) {C[I]=0; D[I]=1;} 
 for (I=2; I<=N; I++) NFACT *= I;

 for (K=1; K<=NFACT; K++)
  {
   PROD=1; ALTPERM(PERM,N);
   for (I=1; I<=N; I++)
    {
     J=PERM[I]-1;
     P = Q + (I-1)*N;
     P = P + J;
     PROD = *P*PROD;
    }
   if (K&1) SUM = SUM - PROD;      /* K is odd */
        else SUM = SUM + PROD;      /* K is even */
  }
 return SUM;
}

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

void ALTPERM(int *A,int N)
    /* Receive integer N, 2<=N<14, and a pointer A to an array of length 14;
       N! repeated calls to ALTPERM produce all permutations of length N;
       The first permutation called is odd; after that the called permutations
        alternate between even,odd; the last permutation is in normal
        increasing order 1,2,...,N and is even */
{
 int I,J,K,COPY;

   K=0; I=N;
  NEW:      /* construct a new PERMutation from the old one */
     C[I] = J = C[I] + D[I];
     if (J==I) {D[I]=-1; goto LAB;}
     if (J) goto FIN;
     D[I]=1; K++;
  LAB:
     if (I>2) {I--; goto NEW;}
     J=1;
  FIN:
     J+=K; COPY=A[J]; A[J]=A[J+1]; A[J+1]=COPY;
 return;
}