#define M 20
#include "stdio.h"
/******************************************************************************************/
void inpmat(double sm[M][M],int n)  /*輸入函數*/
{
int i,j;
double tmpi;
printf("Input elements of matrix:\n");
for(i=0;i<n;i++)
  for(j=i;j<n;j++)
    {if(i==j)
      {printf("a%d%d=",i+1,i+1);
       scanf("%lf",&tmpi);
       sm[i][i]=tmpi;
       }
     else
      {printf("a%d%d=a%d%d=",i+1,j+1,j+1,i+1);
       scanf("%lf",&tmpi);
       sm[i][j]=tmpi,sm[j][i]=tmpi;
       }
     }
/*for(i=0;i<n;i++)
  for(j=0;j<n;j++)
  {printf("a%d%d=",i+1,j+1);
    scanf("%lf",&tmpi);
    sm[i][j]=tmpi;
   }
*/
}
/******************************************************************************************/
void ouptsm(double sm[M][M],int n)  /*輸出原矩陣*/
{
int i,j;
printf("The original matrix is:\n");
for(i=0;i<n;i++)
  {for(j=0;j<n;j++)
    {printf("a%d%d=%.4E\t",i+1,j+1, sm[i][j]);
//     if(n>5)delay(500);
     }
   printf("\n");
  }
}
/******************************************************************************************/
double valdet(double mtx[M][M],int n)  /*求行列式值*/
{
int i,j,cgi,cgj,r,chgop;
double temp,detv;
chgop=0;
if(mtx[0][0]==0)
  {chgop=chgop+1;
   for(cgj=1;cgj<n;cgj++)
    {if(mtx[0][cgj]!=0)
      {for(cgi=0;cgi<n;cgi++)
         temp=mtx[cgi][0],mtx[cgi][0]=mtx[cgi][cgj],mtx[cgi][cgj]=temp;
      }
    }
  }
for(r=0;r<n-1;r++)
  {for(i=r+1;i<n;i++)
    {for(j=n;j>=r;j--)
       mtx[i][j]=mtx[i][j]-mtx[i][r]/mtx[r][r]*mtx[r][j];
     if(mtx[r+1][r+1]==0)
       {chgop=chgop+1;
        for(cgj=r+1;cgj<n;cgj++)
         {if(mtx[r+1][cgj]!=0)
	   {for(cgi=0;cgi<n;cgi++)
	      temp=mtx[cgi][r+1],mtx[cgi][r+1]=mtx[cgi][cgj],mtx[cgi][cgj]=temp;
            }
	 }
       }
     }
   }
if(chgop%2!=0)mtx[n-1][n-1]=mtx[n-1][n-1]*-1;
detv=1;
for(r=0;r<n;r++)
  detv=detv*mtx[r][r];
if(detv==0)detv=0;
return(detv);
}
/******************************************************************************************/
void juzeni(double tm[M][M],double dm[M][M],int n)  /*求出代數余子式和逆矩陣*/
{
int fij,fi,fj,ti,tj,i,j;
double tp[M][M],rank;
ti=0,tj=0;
for(i=0;i<n;i++)
  for(j=0;j<n;j++)
    tp[i][j]=tm[i][j];
rank=valdet(tp,n);
printf("The outcome of elementary transform of determinant is:\n");
for(i=0;i<n;i++)
  {for(j=0;j<n;j++)
     printf("a%d%d=%.4E\t",i+1,j+1,tp[i][j]);
   printf("\n");
  }
printf("The value of the determinant of the matrix is:\t%.4E\n",rank);
for(fi=0;fi<n;fi++)
  for(fj=0;fj<n;fj++)
    {for(i=0;i<n;i++)
     if(i!=fi)
      {for(j=0;j<n;j++)
        {if(j!=fj)
          {tp[ti][tj]=tm[i][j];tj++;}
         }
       tj=0;
       ti++;
       }
     ti=0;
     fij=fi+fj;
     dm[fj][fi]=valdet(tp,n-1);
     if(fij%2!=0)
       dm[fj][fi]=dm[fj][fi]*(-1)/rank;
     else
       dm[fj][fi]=dm[fj][fi]/rank;
     }
}
/******************************************************************************************/
void pdzdjz(double tm[M][M],int n)  /*正定矩陣判定*/
{
int i,j,rank;
double tp[M][M],hlsz;
printf("\tPrimary child determinants calculation\n");
for(rank=1;rank<=n;rank++)
  {for(i=0;i<rank;i++)
     for(j=0;j<rank;j++)
       tp[i][j]=tm[i][j];
   hlsz=valdet(tp,rank);
   if(hlsz<=0)
     {for(i=0;i<rank;i++)
       {for(j=0;j<rank;j++)
         printf("a%d%d=%.4E\t",i+1,j+1,tm[i][j]);
        printf("\n");
//        if(n>5)delay(500);
       }
      printf("The value of %d rank of determinant is:\t%.4E\n",rank,hlsz);
      goto labelout;
     }
   else
      printf("The value of %d rank of determinant is:\t%.4E\n",rank,hlsz);
  }
labelout:
if(rank>n)
  printf("\tThis is a definite angle shot matrix.\n");
else
  printf("The value of this rank of primary child determinant is not over 0.\n\tSo this matrix is not a definite angle shot matrix.\n");
}
/******************************************************************************************/
void jzxchn(double sm[M][M],double dm[M][M],int n)  /*原逆矩陣相乘檢驗結果*/
{
int i,j,k;
double tm[M][M],acc;
for(i=0;i<n;i++)
  for(j=0;j<n;j++)
    {tm[i][j]=0;
     for(k=0;k<n;k++)
      {acc=sm[i][k]*dm[k][j];
       tm[i][j]=tm[i][j]+acc;
      }
    }
printf("The product of original matrix and negative matrix is:\n");
for(i=0;i<n;i++)
  {for(j=0;j<n;j++)
     printf("a%d%d=%.4E\t",i+1,j+1,tm[i][j]);
   printf("\n");
  }
}
/******************************************************************************************/
main()  /*主函數*/
{
int i,j,n;
double sm[M][M],dm[M][M],tm[M][M],detv;
printf("**********************************BEGIN**********************************\n");
n=0;
while(n<=0||n>20)
 {printf("Input rank of matrix(0<n<20):");
  scanf("%d",&n);
 }
inpmat(sm,n);/*輸入矩陣*/
ouptsm(sm,n); /*輸出原矩陣*/
printf("Press any key to continue...");
getchar();
getchar();
for(i=0;i<n;i++)
  for(j=0;j<n;j++)
    tm[i][j]=sm[i][j];
pdzdjz(tm,n); /*正定矩陣判定*/
printf("Press any key to continue...");
getchar();
juzeni(tm,dm,n);/*求出逆矩陣*/
printf("The negative matrix is:\n");
for(i=0;i<n;i++)
 {for(j=0;j<n;j++)
   {printf("a%d%d=%.4E\t",i+1,j+1,dm[i][j]);
//    if(n>5)delay(500);
    }
  printf("\n");
}
printf("Press any key to continue...");
getchar();
jzxchn(sm,dm,n);/*逆矩陣與原矩陣相乘判斷結果*/
printf("-----------------------------------END-----------------------------------\n Press any key to exit.");
getchar();
return (0);
}