#include <iostream.h>
#include <stdio.h>
#include <math.h>

const int size = 56;
const double delta = 0.000001;

void  matrix_inverse(double *matrix, int size, double *inv_matrix);

void main()
{
	//構(gòu)造五對角矩陣A
	double *A = new double[size*size];
	double *inv_A = new double[size*size];
	double a1 = 1.8;
	double a2 = -0.5;
	double a3 = 0.1;
	double *data = new double[size];//data等于五對角矩陣的第一行
	data[0] = a1;
	data[1] = a2;
	data[size-1] = a2;
	data[2] = a3;
	data[size-2] = a3;
	for (int i = 3; i < size-2; i++)
	{
		data[i] = 0;
	}
    //五對角矩陣：上一行的最后一個元素作為下一行的第一個元素，其余元素順序向后移動一位
	for (int h = 0; h < size; h++)
	{
		for (int w = 0; w < size; w++)
		{
			int pos = (w-h+size)%size;
			A[h*size+w] = data[pos]; 
		}
	}
	delete []data;

	
	cout<<"原矩陣A"<<'\n';
	for (int m = 0; m < size; m++)
	{
		for (int n = 0; n < size; n++)
		{
			cout<<A[m*size+n]<<'\t';
		}
		cout<<'\n';
	}



	matrix_inverse(A, size, inv_A);
	
	cout<<"inv_A"<<'\n';
	for (m = 0; m < size; m++)
	{
		for (int n = 0; n < size; n++)
		{
			cout<<inv_A[m*size+n]<<'\t';
		}
		cout<<'\n';
	}
	delete []A;
	delete []inv_A;
	
}


double  matrix_value (double *matrix, int size)
{
	int i,j,m,n,s,t,k=1;
	double f=1,c,x,sn;
	double *matrix1 = new double[size*size];
	for (i = 0; i < size; i++)
	{
		for (j = 0; j < size; j++)
		{
			matrix1[i * size +j] = matrix[i * size +j];
		}
	}
	for (i=0,j=0;i<size&&j<size;i++,j++)
	{
		if (matrix1[i*size+j]==0)
		{
			for (m=i;fabs(matrix1[m*size+j])<delta && m < size;m++)
			if (m==size)
			{
				sn=0;
				//cout<<sn<<'\n';
			}
			else
				for (n=j;n<size;n++)
				{
					c=matrix1[i*size+n];
					matrix1[i*size+n]=matrix1[m*size+n];
					matrix1[m*size+n]=c;
				}
				k*=(-1);
		}
		for (s=size-1;s>i;s--)
		{
			x=matrix1[s*size+j];
			for (t=j;t<size;t++)
				matrix1[s*size+t]-=matrix1[i*size+t]*(x/matrix1[i*size+j]);
		}
	}
	for (i=0;i<size;i++)
	{
		f*=matrix1[i*size+i];
	} 
	sn=k*f;
	return sn;
	delete []matrix1;
}


void matrix_cofactor (double *matrix, int size, double *cof_matrix)
{
	//double *matrix         原矩陣
	//int size              方陣的維數(shù)
	//double *cof_matrix     代數(shù)余子式
	double *temp = new double[(size-1)*(size-1)];
	int z,j,k,l,m,g;
	
	for (z = 0; z < size; z++)
	{
		l = z;
		for (j = 0;j<size;j++)
		{
			m = j;
			for (k = 0;k<size-1; k++)
			{
				for (g = 0;g<size-1;g++)
				{
					if (g >= m && k < l)
					{
						temp[k*(size-1)+g] = matrix[k*size+g+1];
					}
					else if (k>=l && g<m)
					{
						temp[k*(size-1)+g] = matrix[(k+1)*size+g];
					}
					else if (k>=l && g>=m)
					{
						temp[k*(size-1)+g] = matrix[(k+1)*size+g+1];
					}
					else
						temp[k*(size-1)+g] = matrix[k*size+g];
				}//g
			}//k			
			cof_matrix[z*size+j] = matrix_value (temp, size-1);			
		}//j
	}//z	

	for (z = 0; z<size; z++)
	{
		for (j = 0; j<size; j++)
		{			
		    if ((z+j)%2!=0 && cof_matrix[z*size+j]!=0)
			{
				cof_matrix[z*size+j] = (-1)*cof_matrix[z*size+j];
			}			
		}
	}



	/*利用性質(zhì)：循環(huán)矩陣的逆仍然是循環(huán)矩陣	
	    //只求第一行
    	l = 0;
		for (j = 0;j<size;j++)
		{
			m = j;
			for (k = 0;k<size-1; k++)
			{
				for (g = 0;g<size-1;g++)
				{
					if (g >= m)//右側(cè)
					{
						temp[k*(size-1)+g] = matrix[(k+1)*size+g+1];
					}
					else//左側(cè)
					{
						temp[k*(size-1)+g] = matrix[(k+1)*size+g];
					}
				}//g
			}//k
			
			cof_matrix[j] = matrix_value (temp, size-1);
			if (j%2!=0 && cof_matrix[j]!=0)
			{
				cof_matrix[j] = (-1)*cof_matrix[j];
			}
		}//j

		//循環(huán)矩陣：上一行的最后一個元素作為下一行的第一個元素，其余元素順序向后移動一位
	for (int h = 1; h < size; h++)
	{
		for (int w = 0; w < size; w++)
		{
			int pos = (w-h+size)%size;
			cof_matrix[h*size+w] = cof_matrix[pos]; 
		}
	}
*/

	
	
	delete []temp;


}


void  matrix_inverse(double *matrix, int size, double *inv_matrix)
{
	//matrix       原矩陣(大小為size*size)
	//inv_matrix   逆矩陣	
	double value = 0;//行列式的值
	double temp = 0;//用來對矩陣進行轉(zhuǎn)置

	int z,j;

	//value = matrix_value (matrix, size, size);
	value = matrix_value(matrix, size); 
	cout<<"value = "<<value<<'\n';

	//inv_matrix 中存放原矩陣各元素對應(yīng)的“余子式”
	matrix_cofactor (matrix, size, inv_matrix);
	
	/*
	for (z = 0; z<size; z++)
	{
		for (j = 0; j<size; j++)
		{			
		    if ((z+j)%2!=0 && inv_matrix[z*size+j]!=0)
			{
				inv_matrix[z*size+j] = (-1)*inv_matrix[z*size+j];
			}			
		}
	}
	*/
	

    cout<<"代數(shù)余子式："<<'\n';
	for (z = 0; z<size; z++)
	{
		for (j = 0; j<size; j++)
		{			
		    cout<<inv_matrix[z*size+j]<<'\t';			
		}
		cout<<'\n';
	}


	//對inv_matrix轉(zhuǎn)置，此時inv_matrix中存放原矩陣的伴隨矩陣
	for (z = 0; z<size; z++)
	{
		for (j = z+1; j<size; j++)
		{
			temp = inv_matrix[z*size+j];
			inv_matrix[z*size+j] = inv_matrix[j*size+z];
			inv_matrix[j*size+z] = temp;
		}
	}


/*
	cout<<"伴隨矩陣："<<'\n';
	for (z = 0; z<size; z++)
	{
		for (j = 0; j<size; j++)
		{			
		    cout<<inv_matrix[z*size+j]<<'\t';			
		}
		cout<<'\n';
	}
*/


	//求逆矩陣，此時inv_matrix中存放原矩陣的逆矩陣
	for (z = 0; z<size; z++)
	{
		for (j = 0; j<size; j++)
		{
			inv_matrix[z*size+j] = inv_matrix[z*size+j]/value;
		}
	}

}