#include<stdio.h>
#include<malloc.h>

double comput_D(double *q, int c);                                               /*求矩陣的行列式，p為地址，c為行*/
double * Turn(double * q ,int row ,int line );                                   /*求矩陣轉置*/
double * multiply(double * M1,int row1,int line1,double * M2,int row2,int line2);/*求矩陣的乘積*/
double Creat_M(double *p,int m,int n, int k);                                    /*求代數余子式*/
double * Qiuni(double * q,int row );                                             /*求矩陣的逆*/

void main()                                                            
{
	int i,j;
	int n,m;                                                                    /*n系統的階次,設定起始階次為2,m為p矩陣的總數據*/
    double sum;                                                                 /*矩陣的行列式值*/
	double * start;                                                             /*起始地址*/   
	double * p,* Y,* E,*A;                                                      /*p為輸入u(0),u(1),u(2)...u(3n-2)和輸出數據y(0),y(1),y(2)....y(3n-1)所構成的2n×2n矩陣 ,
						                                                          Y為[y(n),y(n+1)....y(3n+1)]的轉置, E為殘差,A為參數矩陣*/
	/********************************/
	/*初始化輸入輸出并賦值*/
	/********************************/
	
	double y[41]={1,0.5,-0.435,1.1649,4.4645,6.667,5.7681,1.1378,-7.4876,-16.8768,
		-16.6648,-11.9831,-2.253,5.4578,7.6442,7.4855,-2.1909,-11.5583,-14.1677,-11.7752,
        -5.9268,0.5166,7.1186,10.8402,7.5873,1.172,-4.1787,-6.6867,-2.8530,4.0195,
		6.3653,5.4713,6.4315,4.091,0.3473,0.9958,3.6766,4.9766,0.7394,-10.1271,-15.1661};
 
	double u[41]={-0.3,0.2,0.99,0.16,0.68,-0.87,0.22,-3.1,0.76,-0.75,
		-0.24,1.18,-1.92,1.81,-1.75,-1.49,0.28,-0.51,-0.08,-0.28,
		0.2,0.9,-0.28,-0.63,0.31,-0.47,0.55,1.05,0.25,-1.06,
		1.4,0.34,-0.76,0.84,0.98,0.19,0.19,-2.18,-1.39,1.36,-0.1};


    /**************開辟空間*****************/
	p=(double *)calloc(40*40,sizeof(double));
	Y=(double *)calloc(40*40,sizeof(double));
	E=(double *)calloc(40*40,sizeof(double));
	A=(double *)calloc(40*40,sizeof(double));

    start=p;
	/*開始循環，設本系統的初始階次為2階*/	
    for(n=1;n<=2;n++)
	{
		m=2*n;	
		/*************p矩陣賦值********************/		
		/********先給矩陣賦輸出數據y（t）**********/                
		for(i=0;i<2*n;i++)                 		
		{
			for(j=0;j<n;j++)			
			{
				*(start+i*2*n+j)=y[j+i];
			}
		}    
		/********給矩陣賦值輸入數據u（t）*********/	
		for(i=0;i<2*n;i++)
		{
			for(j=n;j<2*n;j++)
			{
				*(start+i*2*n+j)=u[j-n+i];
			}
		}



		for(i=0;i<2*n;i++)
		{
			for(j=0;j<2*n;j++)
				printf("%f",*(start+i*2*n+j));
			printf("\n");
		}
		sum=comput_D(p,m);  
		/*由遞退方法求你n可知，判斷sum是否為0，如果sum不為0，矩陣p為滿秩，則n即為所要求的階次，否則繼續循環*/	    
		if(sum!=0)
		{
			printf("\nThe rank of this systerm is:%d",n);
			printf("\nThe value of Matrix is :%f:",sum);
		}
	}
    free(p);
	free(Y);
	free(E);
	free(A);
}


/*****************************************/
/*功能：求矩陣 c X c 的行列式*/
/*入口參數：矩陣首地址 q；矩陣行數 c*/
/*返回值：矩陣的行列式值*/
/*****************************************/
double comput_D(double *q,int c)  
{
     int i,j,m;         /*i--row; j--column*/
     int lop=0;
     double result=0;
     double mid=1;
   
     if (c!=1)
     {
         lop=(c==2)?1:c;     /*控制求和循環次數，若為2階，則循環1次，否則為c次*/

         for(m=0;m<lop;m++)
         {
             mid=1;          /*順序求和*/
             for(i=0,j=m;i<c;i++,j++)
                 mid = mid * ( *(q+i*c+j%c) );
             result+=mid;
         }

         for(m=0;m<lop;m++)
         {                       
             mid=1;          /*逆序相減*/
             for(i=0,j=c-1-m+c; i<c; i++,j--)
                 mid=mid * ( *(q+i*c+j%c));
             result-=mid;
         }
	 }
     else result=*q;
     return(result);
}

/*****************************************/
/*功能：求矩陣 row X line的轉置*/
/*入口參數：矩陣首地址 q；矩陣行數 row,矩陣的列數line*/
/*返回值：矩陣的行列式值*/
/*****************************************/
double * Turn(double * q ,int row ,int line )
{ 
	int i,j;
	double * t;

	t=(double *)calloc(2*row*line,sizeof(double));	
	for(i=0;i<row;i++)
        for(j=0;j<line;j++)        
		  {
		     *(t+j*row+i)=*(q+i*line+j);

		  }
	return(t);
    free(t);
}
/*****************************************/
/*功能：求矩陣 A:row1行，line1列和矩陣B：row2行，line2列的乘積*/
/*入口參數：矩陣A首地址M1；矩陣B的首地址M2；A的行數row1,A的列數line1，B的行數row2,B的列數line2*/
/*返回值：矩陣的乘積值*/
/*****************************************/
double * multiply(double * M1,int row1,int line1,double * M2,int row2,int line2)
{
	double * C;
	int i,j,k,num3;
  	
	num3=2 * row1 * line2;
	C = (double *)calloc(num3, sizeof(double));
	for(i=0;i<row1;i++)          
	        for(j=0;j<line2;j++)
		       *(C+i*line2+j)=0;	
	for(i=0;i<row1;i++)
		{
			for(j=0;j<line2;j++)
			{
				for(k=0;k<line1;k++)
					*(C+i*line2+j)+=(*(	M1+i*line1+k))*(*(M2+k*line2+j));
			}
		}   
    return(C);
	free(C);
}
/*****************************************/
/*功能：求任何一個矩陣的逆*/
/*入口參數：矩陣首地址 q；矩陣行數row*/
/*返回值：方陣的逆*/
/*****************************************/
double * Qiuni(double * q,int row )
{
     double *p;                                      /*定義數組首地址指針變量*/
     int num;                                        /*定義矩陣行數row及矩陣元素個數*/
     int i,j;
     double determ;                                  /*定義矩陣的行列式*/
	 
     num=2 * row * row;
     p = (double *)calloc(num, sizeof(double));     /*分配內存單元*/
     determ=comput_D(q,row);                        /*求整個矩陣的行列式*/
     p=q + row * row;
     if (determ != 0)
     {
         for (i=0;i<row; i++)                      /*求逆矩陣*/
             for (j=0; j<row; j++)
                    *(p+j*row+i)=   Creat_M(q,i,j,row)/determ;   
            
         printf("The determinant is %G\n",determ);

         p=q + row * row;
         printf("\nThe inverse matrix is:\n");
         return(p);                               /*返回該矩陣*/
     }
     else
         return(0);
		 free( p );
}

/******************************************************/
/*功能：求k×k矩陣中元素A(mn)的代數余子式*/
/*入口參數：k×k矩陣首地址；元素A的下標m,n; 矩陣行數 k*/
/*返回值： k×k矩陣中元素A(mn)的代數余子式*/
/*****************************************************/
double Creat_M(double *p,int m,int n,int k)
{
     int len;
     int i,j;
     double mid_result=0;
     int quo=1;
     double *p_creat,*p_mid;

     len=(k-1)*(k-1);
     p_creat = (double *)calloc(len, sizeof(double));     
     p_mid=p_creat;
     for(i=0;i<k;i++)
         for(j=0;j<k;j++)
         {
             if (i!=m && j!=n)
                 *p_mid++ =* (p+i*k+j);            
         }
     quo = (m + n) %2==0 ? 1:-1;
     mid_result = (double ) quo * comput_D(p_creat,k-1);
     free(p_creat);
     return(mid_result);
}