#include<stdio.h>
#include<math.h>
#define MAX_N 20
#define EPS 0.001
#define Y 2
void main()
{ 
	int i,j,k,n,I0,r,t=0;
    float a[MAX_N][MAX_N],d[MAX_N][MAX_N],x[MAX_N],p,W,w,c,s[MAX_N],m[MAX_N],z[MAX_N];
	printf("\n請輸入系數矩陣的階數n:");
	scanf("%d",&n);
	printf("請輸入增廣矩陣:\n");
	for(i=1;i<=n;i++)
		for(j=1;j<=n+1;j++)
			scanf("%f",&a[i][j]);
    for(i=1;i<=n;i++)
		for(j=1;j<=n+1;j++)
			d[i][j]=a[i][j];
	//列主元解方程
		for(k=1;k<=n-1;k++)
		{
			p=a[k][k];
			I0=k;
			for(i=k;i<=n;i++)
			{
				if(fabs(a[i][k])>fabs(p))
				{
					p=a[i][k];
					I0=i;
				}
			}
		    if(fabs(p)<=EPS)
				printf("EXI=1");
			if(I0==k)
			{
				for(i=k+1;i<=n;i++)
				{
					a[i][k]=a[i][k]/a[k][k];
				    for(j=k+1;j<=n+1;j++)
						a[i][j]=a[i][j]-a[i][k]*a[k][j];
				}
			}
		    else
			{
				for(j=k;j<=n+1;j++)
				{
					w=a[k][j];
			        a[k][j]=a[I0][j];
				    a[I0][j]=w;
				}
	            for(i=k+1;i<=n;i++)
				{
					a[i][k]=a[i][k]/a[k][k];
				    for(j=k+1;j<=n+1;j++)
						a[i][j]=a[i][j]-a[i][k]*a[k][j];
				}
			}
		}
	if(a[n][n]==0)
		printf("EXI=1");
	else
	{
	  a[n][n+1]=a[n][n+1]/a[n][n];
	  for(k=n-1;k>=1;k--)
	  {
		  W=0;
		  for(j=k+1;j<=n;j++)
			  W=W+a[k][j]*a[j][n+1];
		  a[k][n+1]=a[k][n+1]-W;
		  a[k][n+1]=a[k][n+1]/a[k][k];
	  }
	}
	for(k=1;k<=n;k++)
		x[k]=a[k][n+1];
   //迭代
	while(t!=Y)                 //循環控制
   {
		for(i=1;i<=n;i++)               //Doolittle解方程
	    	for(j=1;j<=n+1;j++)
		    	a[i][j]=d[i][j];
	    for(i=1;i<=n;i++)
		{
			m[i]=0;
			for(k=1;k<=n;k++)
				m[i]=m[i]+a[i][k]*x[k];
		}
		for(i=1;i<=n;i++)
			a[i][n+1]=a[i][n+1]-m[i];
		for(r=1;r<=n-1;r++)
		{
			for(i=r;i<=n;i++)
			{
				s[i]=a[i][r];
			    for(k=1;k<=r-1;k++)
				{
					s[i]=s[i]-a[i][k]*a[k][r];
				}
		        a[i][r]=s[i];
			}
		    p=0;
	    	for(i=r;i<=n;i++)
			{
				if(fabs(s[i])>p)  
					p=s[i];
			        I0=i;
			}
      		if(fabs(p)<=EPS)
				printf("EXI=1");
		    if(I0!=r)
			{
				for(j=1;j<=n+1;j++)
				{
					w=a[r][j];
			        a[r][j]=a[I0][j];
				    a[I0][j]=w;
				}
			}
	    	for(i=r+1;i<=n;i++)
				a[i][r]=a[i][r]/a[r][r];
		    for(j=r+1;j<=n+1;j++)
			{
				c=a[r][j];
			    for(k=1;k<=r-1;k++)
				{		
					c=c-a[r][k]*a[k][j];
				}
		        a[r][j]=c;
			}
		}
	    for(j=n;j<=n+1;j++)
		{
			c=a[n][j];
		    for(k=1;k<=n-1;k++)
			{
				c=c-a[n][k]*a[k][j];
			}
	        a[n][j]=c;
		}
        if(fabs(a[n][n])<=EPS)
			printf("EXI=1");
	    a[n][n+1]=a[n][n+1]/a[n][n];
        for(i=n-1;i>=1;i--)
		{
			for(k=i+1;k<=n;k++)
				a[i][n+1]=a[i][n+1]-a[i][k]*a[k][n+1];
		    a[i][n+1]=a[i][n+1]/a[i][i];
		}
		for(i=1;i<=n;i++)
			z[i]=a[i][n+1];
    	for(i=1;i<=n;i++)
	    	x[i]=x[i]+z[i];
        t=t+1;
	}
	for(i=1;i<=n;i++)
		printf("x[%d]=%f\n",i,x[i]);
}