#include"stdio.h"
#include"time.h"

double Surplus(double A[],int m,int n) /*求矩陣行列式*/ 
{ 
	
	int i,j,k,p,r; 
	double X,temp=1,temp1=1,s=0,s1=0; 
	
	if(n==2) 
	{for(i=0;i<m;i++) 
	for(j=0;j<n;j++) 
		if((i+j)%2) temp1*=A[i*n+j]; 
		else temp*=A[i*n+j]; 
		X=temp-temp1;} 
	else{ 
		for(k=0;k<n;k++) 
		{for(i=0,j=k;i<m,j<n;i++,j++) 
		temp*=A[i*n+j]; 
		if(m-i) 
		{for(p=m-i,r=m-1;p>0;p--,r--) 
		temp*=A[r*n+p-1];} 
		s+=temp; 
		temp=1; 
		} 
		
		for(k=n-1;k>=0;k--) 
		{for(i=0,j=k;i<m,j>=0;i++,j--) 
		temp1*=A[i*n+j]; 
		if(m-i) 
		{for(p=m-1,r=i;r<m;p--,r++) 
		temp1*=A[r*n+p];} 
		s1+=temp1; 
		temp1=1; 
		} 
		
		X=s-s1;} 
	return X; 
} 
double * MatrixInver(double A[],int m,int n) /*矩陣轉置*/ 
{ 
	int i,j; 
	double *B=NULL; 
	B=(double *)malloc(m*n*sizeof(double)); 
	
	for(i=0;i<n;i++) 
	{
		for(j=0;j<m;j++) 
			B[i*m+j]=A[j*n+i]; 
	}
	return B; 
} 
double * MatrixOpp(double A[],int m,int n) /*逆矩陣*/
{ 
	int i,j,x,y,k; 
	double *SP=NULL,*AB=NULL,*B=NULL,X,*C; 
	SP=(double *)malloc(m*n*sizeof(double)); 
	AB=(double *)malloc(m*n*sizeof(double)); 
	B=(double *)malloc(m*n*sizeof(double)); 
	
	X=Surplus(A,m,n); 
	X=1/X; 
	
	for(i=0;i<m;i++) 
	{ 
		for(j=0;j<n;j++) 
		{
			for(k=0;k<m*n;k++) 
				B[k]=A[k]; 
			{
				for(x=0;x<n;x++) 
					B[i*n+x]=0; 
				for(y=0;y<m;y++) 
					B[m*y+j]=0; 
				B[i*n+j]=1; 
				SP[i*n+j]=Surplus(B,m,n); 
				AB[i*n+j]=X*SP[i*n+j]; 
			} 
		}
	}
	C=MatrixInver(AB,m,n); 
	
	return C; 
}

double * MatrixMul(double A[],double B[],int m,int n,int p) /*矩陣相乘*/ 
// A:m×n  B:n×p  C:m×p  C等于A乘以B
{ 
	int i,j,k; 
	double *C=NULL; 
	C=(double*)malloc(m*p*sizeof(double)); 
	
	for(i=0;i<m;i++) 
		for(j=0;j<p;j++) 
		{   C[i*p+j]=0;
			for(k=0;k<n;k++)
				C[i*p+j]=C[i*p+j]+A[i*n+k]*B[k*p+j]; 
		}

	return C; 
} 

double * MatrixMulNum(double A[],int m,int n,double num) /*矩陣數乘*/ 
// A:m×n, m行n列
{ 
	int i,j; 
	double *C=NULL; 
	C=(double*)malloc(m*n*sizeof(double)); 	
	for(i=0;i<m;i++) 
		for(j=0;j<n;j++) 
		  C[i*n+j]=num*A[i*n+j];
	return C; 
}

double * MatrixAdd(double A[],double B[],int m,int n) /*矩陣相加*/ 
// A:m×n  B:m×n C:m×n  C等于A加B
{ 
	int i,j; 
	double *C=NULL; 
	C=(double*)malloc(m*n*sizeof(double)); 	
	for(i=0;i<m;i++) 
		for(j=0;j<n;j++) 
		   C[i*n+j]=A[i*n+j]+B[i*n+j];		
	return C; 
} 

double * MatrixSub(double A[],double B[],int m,int n) /*矩陣相減*/ 
// A:m×n  B:m×n C:m×n  C等于A減B
{ 
	int i,j; 
	double *C=NULL; 
	C=(double*)malloc(m*n*sizeof(double)); 	
	for(i=0;i<m;i++) 
		for(j=0;j<n;j++) 
		   C[i*n+j]=A[i*n+j]-B[i*n+j];		
	return C; 
} 

double Matrix2Num(double A[])/*將一行一列的矩陣A轉化為數*/
{
	return A[0];
}

double * MatrixRand(int m,int n)/*用(-1,1)之間的隨機數生成填充一個m×n的矩陣*/
{   int i,j;
    double* C=NULL;
	C=(double*)malloc(m*n*sizeof(double));
	srand((unsigned)(time(NULL)));
    for (i=0;i<m;i++) 
	{
		for (j=0;j<n;j++) 
		{
			//rand()函數產生的隨機數是整數，所以需要將之做適當的變換，變到-1～1之間
			C[i*n+j]= (double)(rand())/(32767/2) - 1;
		}
	}
	return C;	
}

double ABS(double a) /*返回a的絕對值*/
{  if(a>=0) return a;
   else return -a;
}

int* Rand2(int a,int b)/*產生b個0到a-1的隨機數*/
{   int i;
    int* C=NULL;C=(int*)malloc(b*sizeof(int));
	srand((unsigned)(time(NULL)));
	for (i=0;i<b;i++) 
		{
			//rand()函數產生的隨機數是整數，所以需要將之做適當的變換，變到0到a-1之間
			C[i]=int(a*rand()/32767);
		}
	return C;
}