#include<stdio.h>
#include<malloc.h>
#include<math.h>
#include<iostream.h>

#define MAX_N 20

void main()
{
	cout<<"==========================================="<<endl;
	cout<<"     Courant求解線性方程組AX=b"<<endl;
	cout<<"==========================================="<<endl;
    cout<<endl;
	int n,i,j,k;
	double **a,*b,*x,*y,**l,**u;
	cout<<"輸入階數(shù)n  (0<n<=20):"<<endl;
	cin>>n;
	if(n > MAX_N || n <=0)
	{
		cout<<"n 不合法!"<<endl;
		return;
	}
	a = (double **)malloc(sizeof(double*)*n);
	for(i = 0;i<n;i++)
	{
		a[i] = (double *)malloc(sizeof(double)*n);
	}
	b= (double *)malloc(sizeof(double)*n);
	x = (double *)malloc(sizeof(double)*n);
	y = (double *)malloc(sizeof(double)*n);
	u = (double **)malloc(sizeof(double*)*n);
	for(i = 0;i<n;i++)
	{
		u[i] = (double *)malloc(sizeof(double)*n);
	}
	l = (double **)malloc(sizeof(double*)*n);
	for(i = 0;i<n;i++)
	{
		l[i] = (double *)malloc(sizeof(double)*n);
	}
	cout<<"輸入系數(shù)矩陣a(i,j)(按行輸入)(每輸入一個即回車)"<<endl;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			cin>>a[i][j];

	cout<<"輸入 b(j)(每輸入一個即回車)"<<endl;
	for(i = 0;i<n;i++)
		cin>>b[i];
	for(i = 0;i<n;i++) u[i][j]=1;
	for(k=0;k<n;k++)
	{
		for(i=k;i<n;i++)
		{
			l[i][k] = a[i][k];
			for(j=0; j< k; j++)
				l[i][k] -= (l[i][j] * u[j][k]);
		}

		for(j = k+1; j<n ; j++)
		{	
			u[k][j] = a[k][j];
			for(i = 0; i < k; i++)
				u[k][j] -= (l[k][i]*u[i][j]);
			u[k][j] /= l[k][k];
		}

	}
	for(i = 0; i<n ; i++)
	{
		y[i] = b[i];
		for(j = 0; j < i; j++)
			 y[i] -= (l[i][j] * y[j]);
		y[i] /= l[i][j];
	}

	for(i = n-1; i>=0 ; i--)
	{
		x[i] = y[i];
		for(j = i+1;j<n;j++)
			x[i] -= (u[i][j] * x[j]);
	}

	cout<<"方程組的解為:"<<endl;
	for(i = 0; i < n-1; i++)
	{
		cout<<x[i]<<" , ";
	}
    cout<<x[n-1]<<endl;

}
/*
  完成人：pb03011103  肖璟
*/