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

void col_expunction(int n,double** a,double* b)    //列主消元法函數
{
	double* x=new double[n];
	int* q=new int [n];
	for(int i=0;i<n;i++)
		q[i]=i;

	if(n==1)x[0]=b[0]/a[0][0];                     //當n=1時的處理

	for(int k=0;k<n-1;k++){                       //選最大系數列
		double max=fabs(a[q[k]][k]); int p=q[k];
		for(int j=k+1;j<n;j++){
			if(max<fabs(a[q[j]][k])){
				max=a[q[j]][k];
				p=q[j];
			}
		}

		for(i=k+1;i<n;i++){
			if(q[i]==p){
				q[i]=q[k];
				q[k]=p;
				break;
			}
		}

		for(i=k+1;i<n;i++){                          //方程求解過程
			double l=(a[q[i]][k])/a[p][k];
			for(j=k+1;j<n;j++)
				a[q[i]][j]=a[q[i]][j]-l*(a[p][j]);
			b[q[i]]=b[q[i]]-l*b[p];
		}

		x[n-1]=b[q[n-1]]/a[q[n-1]][n-1];
		for(i=n-2;i>=0;i--){
			double y=0;
			for(j=i+1;j<n;j++)
				y+=a[q[i]][j]*x[j];
			x[i]=(b[q[i]]-y)/a[q[i]][i];
		}
	}

	for(i=0;i<n;i++)        
		b[i]=x[i];
}


void main()
{
	char ch;
	do{
		int n;
		cout<<"Input n:";
		cin>>n;

		double** a=new double* [n];
		for(int i=0;i<n;i++)
			a[i]=new double[n];
		double* b=new double [n];
			
		cout<<endl<<"自動生成Hilbert矩陣?(y=yes||anykey=no):";
		cin>>ch;

		if(ch!='y'&&ch!='Y'){
			cout<<"\nInput a[i][j](1<=i,j<=n):\n";
			for(i=0;i<n;i++)
				for(int j=0;j<n;j++)
					cin>>a[i][j];
		
			cout<<"\nInput b[i](1<=i<=n):\n";
			for(i=0;i<n;i++)
				cin>>b[i];
		}
		else{
		
			for(i=0;i<n;i++)
				for(int j=0;j<n;j++)
					a[i][j]=(double)1/(i+j+1);

			for(i=0;i<n;i++){
				b[i]=0;
				for(int j=0;j<n;j++)
					b[i]+=(double)(j+1)/(i+j+1);
			}
		}

		col_expunction(n,a,b);

		cout<<endl<<"方程組的解為:"<<endl;
		for(i=0;i<n;i++)
			cout<<"x["<<i<<"]="<<b[i]<<"\t";
		cout<<endl;

		cout<<endl<<"Continue?(y=yes||anykey=no):";
		cin>>ch;
		cout<<endl;

	}while(ch=='y'||ch=='Y');

}