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


void LU(double**A,double*B,int n)
{
	int i,j,k;
	double**L,**U,*Y,*X;
	U=new double*[n];
	L=new double*[n];
	Y=new double[n];
	X=new double[n];
  	for(i=0;i<n;i++)
	{
		L[i]=new double[n];
		U[i]=new double[n];
	}
	for(k=0;k<n;k++)  //LU分解    
	{
		L[k][k]=1;
		for(j=0;j<n;j++)
		{
			if(j>k) L[k][j]=0;
			if(j<k) U[k][j]=0;
		}
		for(j=k;j<n;j++)
		{
			for(i=0;i<k;i++)
			    A[k][j]-=L[k][i]*U[i][j];			
			U[k][j]=A[k][j];
		}
		for(i=k+1;i<n;i++)
		{
			for(j=0;j<k;j++)
				A[i][k]-=L[i][j]*U[j][k];
			L[i][k]=A[i][k]/U[k][k];
		}
	}
//  ---------

	for(k=0;k<n;k++)        // solve  X
	{
		Y[k]=B[k];
		for(j=0;j<k;j++)
		{
			Y[k]-=L[k][j]*Y[j];
		}
	}
	for(k=n-1;k>-1;k--)
	{
		for(i=k+1;i<n;i++)
			Y[k]-=U[k][i]*X[i];
		X[k]=Y[k]/U[k][k];
	}
	for(i=0;i<n;i++)
	{
		B[i]=X[i];
	}
}
void ccout(double*X,int n)
{
	int i;
	cout<<"The Result is:"<<endl;
	cout<<"  y=";
	cout<<X[0];
	for(i=1;i<n;i++)
	{
		
		if(X[i]>=0)
		
			cout<<"+"<<" ";
		else
			cout<<"";

		cout<<X[i]<<"*x^"<<i;
	}
	cout<<endl;
}

void solve(double **A,double*X1,double*Y1,double*Y , int n,int m)
{
	int i,j;
	double t,T,*D;
	D=new double[2*n+1];
	for(j=0;j<2*n+1;j++) //  求   2n+1 個(gè)數(shù)
	{
		t=0;
		for(i=0;i<m;i++)
		{
			t+=pow(X1[i],j);
		}
		D[j]=t;
		cout<<"D["<<j<<"]="<<D[j]<<endl;
	}
	for(j=0;j<n+1;j++)
	{
		T=0;
		for(i=0;i<m;i++)
		{
			T+=pow(X1[i],j)*Y1[i];		
		}
		Y[j]=T;
	}
	for(i=0;i<n+1;i++)
	{
		cout<<Y[i]<<endl;
	}
	for(i=0;i<n+1;i++)
	{
		for(j=i;j<i+n+1;j++)
		{
			if(i==0)
			{
				A[0][j]=D[j];
			}
			if(i!=0)
			{
				A[i][j-i]=D[j];
			}
		}
	}
	for(i=0;i<n+1;i++)
	{
		for(j=0;j<n+1;j++)
		{
			cout<<"A="<<A[i][j]<<"    ";
		}
		cout<<endl;
	}
}
void main()
{
	int i,n,m;
	double**A,*X1,*Y,*Y1;
//	void LU(double**A,double*B,int n);
//	void solve(double **A,double*X1,double*Y1,double*Y,int n,int m);
//	void ccout(double*X,int n);
/*	cout<<"Input The Index :"<<endl;
	cin>>n;
	cout<<"Input The Mode :"<<endl;
	cin>>m;*/
        n=6;
        m=6;
	A=new double*[n+1];
	X1=new double[n+1];
	Y1=new double[m];
	Y=new double[m];
	for(i=0;i<n+1;i++)
	{
		A[i]=new double[n+1];
	}

	X1[0]=1;X1[1]=1.5;X1[2]=2;X1[3]=34;X1[4]=4;X1[5]=5;
	Y1[0]=0.01466671;Y1[1]=0.116352;Y1[2]=0.257403;Y1[3]=0.467174;Y1[4]=0.590668;Y1[5]=0.668793;
	solve(A,X1,Y1,Y,n,m);
	LU(A,Y,n+1);
	ccout(Y,n+1);
}