#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <Malloc.h>
#include <iostream.h>
/*double f(double x)
{
 return exp(x)+10*x-2;
}
double f(double x,double y)
{
   return y-2.0*x/y;
}
double df(int n,double x)
{
	double h=0.0001;
	if (n==1) return (f(x+h)-f(x-h))/(2*h);
	else return df(n-1,x);
} */
void main()
{
	int i,j,k,d,N;
	double **a,*b,temp,max=0,S=0;
	cout<<"******高斯消去法解線性方程組******"<<endl;
	cout<<"請輸入未知數個數：";
	cin>>N;
	a=(double **)malloc((N+1)*sizeof(double *));
    for(i=0;i<=N;i++)
		*(a+i)=(double *)malloc((N+1)*sizeof(double));
	b=(double *)malloc((N+1)*sizeof(double));
	cout<<"請輸入系數增廣矩陣："<<endl;
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++)
			cin>>a[i][j];
		cin>>b[i];
	}
	for(k=1;k<=N;k++)
	{
		//列主元消去
        if(fabs(a[k][k])<0.000001)
		{
			for(i=k+1;i<=N;i++)
                if(fabs(a[i][k])>max) {max=a[i][k];d=i;}
			for(j=k;j<=N;j++)
			{temp=a[k][j];a[k][j]=a[d][j];a[d][j]=temp;}
			max=0;
		}
		//化為上三角矩陣
		for(j=k+1;j<=N;j++)
			a[k][j]=1.0*a[k][j]/a[k][k];
		b[k]=1.0*b[k]/a[k][k];
		for(i=k+1;i<=N;i++)
		{
			for(j=k+1;j<=N;j++)
				a[i][j]=a[i][j]-a[i][k]*a[k][j];
		    b[i]=b[i]-a[i][k]*b[k];
		}
	}
	//求根
	for(i=N-1;i>=1;i--)
	{
		for(j=i+1;j<=N;j++)
            S+=a[i][j]*b[j];
		b[i]=b[i]-S;
		S=0;
	}
	for(i=1;i<=N;i++)
		cout<<"x"<<i<<"="<<b[i]<<endl;
	free(b);
	for(i=0;i<=N;i++)
		free(*(a+i));
}