#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//------------------------------
//Data[i][j];
//i - Stroka  (Row)
//j - Stolbec (Col)
//------------------------------
//Minimum Value On Diagonal Elements
#define MINVAL 0.0000000000000000000000000000001f
//------------------------------
//Boolean Type
typedef enum
{
	FALSE=0,
	TRUE
} BOOL;
//------------------------------
//Matrix Type
struct QMATRIX
{
	float **Data;
	int n;
};
//------------------------------
//Prototypes
void FreeQMatrix(QMATRIX *m);
BOOL LoadQMatrix(char *FileName,QMATRIX *m);
BOOL PrintQMatrix(QMATRIX *m);
BOOL PrintQResult(QMATRIX *m);
BOOL QSwapRows(QMATRIX *m,int i,int j);
BOOL QIterateRows(QMATRIX *m,int From);
BOOL QGaussDown(QMATRIX *m,int from);
BOOL QGaussUp(QMATRIX *m);
BOOL QGauss(QMATRIX *m);
//------------------------------
//Free Matrix
void FreeQMatrix(QMATRIX *m)
{
	int i;

	if(m==NULL)
		return;
	m->n=0;
	if(m->Data==NULL)
		return;
	for(i=0;i<m->n;i++)
		if(m->Data[i]!=NULL)
			free(m->Data[i]);
	free(m->Data);
}
//------------------------------
//Load Matrix From File
BOOL LoadQMatrix(char *FileName,QMATRIX *m)
{
	FILE *f;
	int i,j;

	//Free Last Matrix
	FreeQMatrix(m);

	//Open File
	f=fopen(FileName,"rt");
	if(f==NULL)
		return FALSE;

	//Scan Size
	if(fscanf(f,"%i",&(m->n))!=1)
	{
		m->n=0;
		fclose(f);
		return FALSE;
	}

	//Allocate Memory
	m->Data=(float**)malloc(sizeof(float*)*m->n);
	if(m->Data==NULL)
	{
		m->n=0;
		fclose(f);
		return FALSE;
	}

	for(i=0;i<m->n;i++)
	{
		m->Data[i]=(float*)malloc(sizeof(float)*(m->n+1));
		if(m->Data[i]==NULL)
		{
			//If Error - Kill Allocated Memory
			for(j=i-1;j>=0;j--)
				free(m->Data[i]);
			free(m->Data);
			return FALSE;
		}
	}

	//Load Data From File
	for(i=0;i<m->n;i++)
	for(j=0;j<=m->n;j++)
		if(fscanf(f,"%f",&(m->Data[i][j]))!=1)
		{
			FreeQMatrix(m);
			fclose(f);
			return FALSE;
		}

	//Close File And Exit
	fclose(f);
	return TRUE;
}
//------------------------------
//Printing Matrix To Console
BOOL PrintQMatrix(QMATRIX *m)
{
	int i,j;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	for(i=0;i<m->n;i++)
	{
		printf("%3.3f*X%i",m->Data[i][0],1);
		for(j=1;j<m->n;j++)
			printf("%+3.3f*X%i",m->Data[i][j],j+1);
		printf(" =%3.3f\n",m->Data[i][m->n]);
	}
	return TRUE;
}
//------------------------------
//Print Result To Console
BOOL PrintQResult(QMATRIX *m)
{
	int i;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;

	for(i=0;i<m->n;i++)
		printf("X%i=%3.3f ",i+1,m->Data[i][m->n]);
	printf("\n");

	return TRUE;
}
//------------------------------
//Swap Rows In QMatrix
BOOL QSwapRows(QMATRIX *m,int i,int j)
{
	float *Tmp;

	if((m==NULL)||(i<0)||(j<0)||(i>=m->n)||(j>=m->n)||(m->Data==NULL))
		return FALSE;
	//Swap Data
	Tmp=m->Data[i];
	m->Data[i]=m->Data[j];
	m->Data[j]=Tmp;

	return TRUE;
}
//------------------------------
//Iterate Rows (Recursion):
BOOL QIterateRows(QMATRIX *m,int From)
{
	int i;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	//if End Calculating - Exit
	if(From>=m->n)
		return QGaussDown(m,From);
	//if Can Calculating - Calculate And Exit 
	if(QGaussDown(m,From))
		return TRUE;
	//Else Iterate Rows
	for(i=From+1;i<m->n;i++)
	{
		QSwapRows(m,From,i);
		//And Try Calculate
		if(QGaussDown(m,From))
			return TRUE;
		//QRIterateRows(m,From+1);
		QSwapRows(m,From,i);
	}
	return FALSE;
}
//------------------------------
//Gsuss Pass 1 (recursion):
BOOL QGaussDown(QMATRIX *m,int from)
{
	int i,j;
	float *k;
	
	//If End Calculate - Exit
	if(from>=m->n)
		return TRUE;
	//If Zero On Main Diagonal - Return False
	if(fabs(m->Data[from][from])<=MINVAL)
		return FALSE;
	//Alloc Mem For Coefficients
	k=(float*)malloc(sizeof(float)*m->n);
	if(k==NULL)
		return FALSE;
	//ReCalculate Matrix To New Zero Column
	for(j=from+1;j<m->n;j++)
	{
		k[j-from-1]=m->Data[j][from]/m->Data[from][from];
		for(i=from;i<=m->n;i++)
			m->Data[j][i]=m->Data[j][i]-m->Data[from][i]*k[j-from-1];
	}
	//Try Calculate Next Columns
	if(QIterateRows(m,from+1))
	{
		free(k);
		return TRUE;
	}
	else
	{
		//If Cannot Calculate - Return Back
		for(j=from+1;j<m->n;j++)
		{
			for(i=from;i<=m->n;i++)
				m->Data[j][i]=m->Data[j][i]+m->Data[from][i]*k[j-from-1];
		}
		//Free Coefficients
		free(k);
	}
	//Return False
	return FALSE;
}
//------------------------------
//Gauss Pass 2: Final
BOOL QGaussUp(QMATRIX *m)
{
	int i,j,q;
	float k;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	
	for(j=(m->n-1);j>=0;j--)
	{
		for(q=j-1;q>=0;q--)
		{
			k=m->Data[q][j]/m->Data[j][j];
			for(i=j;i<=m->n;i++)
				m->Data[q][i]=m->Data[q][i]-m->Data[j][i]*k;
		}
		m->Data[j][m->n]=m->Data[j][m->n]/m->Data[j][j];
		m->Data[j][j]=1.0f;
	}
	return TRUE;
}
//------------------------------
//Gausses Method
BOOL QGauss(QMATRIX *m)
{
	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	if(!QIterateRows(m,0))
		return FALSE;
	QGaussUp(m);
	return TRUE;
}
//------------------------------
//Entry Point
void main(void)
{
	QMATRIX m={NULL,0};

	if(!LoadQMatrix("d:\\gin.txt",&m))
	{
		printf("Cannot Load Matrix!\n");
		return;
	}

	PrintQMatrix(&m);
	printf("\n");

	if(QGauss(&m))
		PrintQResult(&m);
	else
		printf("Cannot Evaluate Matrix\n");
	
	FreeQMatrix(&m);
}
//------------------------------