//the pivot gauss-elimination method
//basic conceptions:only diag-elements are non_zero!
#include <iostream.h>
#include <iomanip.h>
#include <math.h>

const int N=4;
const double eps=1e-10;

void main()
{
	double a[N][N+1],x[N];

	cout<<"input coefficient matrix!"<<endl;
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
			cin>>a[i][j];	//input coefficient matrix
	cout<<"input constant of right side!"<<endl;
	for(i=0;i<N;i++)
		cin>>a[i][N];		//input constant of right side

	cout<<"the steps of elimination: "<<endl;
	for(int k=0;k<N;k++)
	{
		//select the maximum-column element
		double max=a[k][k];
		int mLine=k;
		for(int i=k+1;i<N;i++)
			if(fabs(a[i][k])>=max)
			{
				max=a[i][k];
				mLine=i;
			}
		//maximum-column element exception
		if(fabs(max)<eps)
			return;
		if(mLine!=k)
		{
			for(int j=k;j<N+1;j++)	//k <==> mLine
			{
				double temp=a[k][j];
				a[k][j]=a[mLine][j];
				a[mLine][j]=temp;
			}
		}
		//消元
		for(i=0;i<N;i++)	//jordan elimination(i starts from 0)
		{
			if(i==k)
				continue;
			double c=a[i][k]/a[k][k];	//消元系數(shù)
			for(int j=k;j<=N;j++)			//第i 行消元
				a[i][j]-=a[k][j]*c;
		}

		//消元后輸出
		for(i=0;i<N;i++)
		{
			for(int j=0;j<N+1;j++)
				cout<<setw(8)<<a[i][j]<<"  ";
			cout<<endl;
		}
		cout<<endl;
	}

	//只有對角線上的元素非零，故直接求解
	for(i=0;i<N;i++)
		x[i]=a[i][N]/a[i][i];

	//輸出解

	cout<<"The Key to Equations: "<<endl;
	for(k=0;k<N;k++)
		cout<<"x["<<k<<"]= "<<x[k]<<endl;
	cout<<endl;
}