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

#define  DataType long double
#define N 12
int GaussEliminant(DataType a[N][N],DataType b[N],DataType c[N]);//利用高斯消元法解多元一次線性方程組

int main(int argc,char *argv[])
{ 
	DataType control[7][4]= // 6個(gè)控制點(diǎn)的(XYZ)坐標(biāo)
	{
		{0,0,0,0},
		{0,1.311,1.751,1.951},
		{0,5.121,5.531,5.941},
		{0,10.511,10.751,10.191},
		{0,15.311,15.551,15.971},
		{0,20.711,20.551,20.981},
		{0,25.311,25.541,25.951},
	};
	DataType photo_control[7][3] = //6個(gè)控制點(diǎn)在圖像內(nèi)的坐標(biāo)(XY)
	{
		{0,0,0},
		{0,1.211,1.531},
		{0,5.141,5.551},
		{0,10.161,10.751},
		{0,15.181,15.591},
		{0,20.121,20.541},
		{0,25.151,25.561},
	};
	DataType b[N];
	DataType L[N];
	for(int ii=0; ii<N; ii++)
	{
		b[ii] = 0;//b 是方程右邊
		L[ii] = 0;//L 是L系數(shù)
	}
	
	DataType a[N][N];//相當(dāng)于方程Al = -B 左側(cè)的A a[0]行與列不參與運(yùn)算
	for (ii=1; ii<=11; ii++)//初始化為0
		for(int jj=1;jj<=11;jj++)
			a[ii][jj] =0;

	//以下是給a[11][11]賦值過程
	for ( ii=1; ii<=11; ii++)//ii 行下標(biāo)
		for( int jj=1;jj<=11;jj++)//jj 列下標(biāo)
		{
			if (ii%2)//奇數(shù)行
			{
				if (jj<=3)//奇數(shù)行的前3項(xiàng)與實(shí)際控制點(diǎn)坐標(biāo)相同
				{
					a[ii][jj] = control[ii/2+1][jj];
				}
				else if (jj==4)//奇數(shù)行的第4項(xiàng)為1
				{
					a[ii][jj] = 1;
				}
				else if (jj>4&&jj<8)//第4到第8項(xiàng)為0
				{
					a[ii][jj] = 0;
				}
				else if (jj>8&&jj<N)//第8到第11項(xiàng)為一乘積
				{
					a[ii][jj] = control[ii/2+1][jj-8]*photo_control[ii/2+1][1];
				}		
			}
			else if (ii%2 == 0)
			{
				if (jj<=4)
				{
					a[ii][jj] = 0;
				}
				else if (jj>4&&jj<8)
				{
					a[ii][jj] = control[ii/2][jj-4];
				}
				else if (jj==8)
				{
					a[ii][jj] = 1;
				}
				else if (jj>8&&jj<N)
				{
						a[ii][jj] = control[ii/2][jj-8]*photo_control[ii/2][2];
				}
			}				
		}//N維數(shù)組A賦值過程結(jié)束
		//給B數(shù)組賦值
		for ( ii=1;ii<N;ii++)
		{
			if (ii%2)
			{
				b[ii] = photo_control[(ii+1)/2][1];
			}
			else if(ii%2==0)
			{
				b[ii] = photo_control[ii/2][2];
			}
		}
		char szBuff[50];//臨時(shí)數(shù)組用于格式化輸出DOUBLE
		cout<<"方程左邊矩陣A為：";
		for (ii=1;ii<N;ii++)
		{
			cout<<endl;
			for (int jj=1;jj<N;jj++)
			{
				sprintf(szBuff,"%f",a[ii][jj]);
				cout<<szBuff<<" ";
			}
		}
		cout<<endl<<"方程右邊矩陣B為："<<endl;
		for (ii=0;ii<N;ii++)
		{
			cout<<b[ii]<<endl;
		}
		GaussEliminant(a,b,L);
		return 0;
}

int GaussEliminant(DataType a[N][N],DataType b[N],DataType c[N])
{
	DataType d=0,t=0,sum=0;
	int k=1,l=0,i=0,j=0,m,n;
	n=N-1;//n代表方程的組數(shù) 
	while(k<=n)
	{
		//*****選主元素*****
		d=a[k][k];
		l=k;
		for(i=k+1;i<=n;i++)
		{
			if(fabs(a[i][k])>fabs(d))
			{
				d=a[i][k];//現(xiàn)在d中存放的是：a[k]列中最大的系數(shù) 就是對(duì)角的那一列
				l=i;// l 中存放：最大數(shù)的行標(biāo)
			}
		}
		if(l!=k)//如果對(duì)角元素在這一列中不是最大的
		{
			for(j=k;j<=n;j++)//則將對(duì)角元素所在式 與 最大元素所在式進(jìn)行 系數(shù) 交換 //兩式對(duì)調(diào)
			{					//目的是：使最大元素在對(duì)角線上
				t=a[l][j];
				a[l][j]=a[k][j];
				a[k][j]=t;
			}
			//下面三步是交換方程 等號(hào) 右面的系數(shù)
			t=b[k];
			b[k]=b[l];
			b[l]=t;
		}
		//*****消元*****
		for(j=k+1;j<=n;j++)
		{
			a[k][j]=a[k][j]/a[k][k];
		}
		b[k]=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];
			}
			j=1;
			b[i]=b[i]-a[i][k]*b[k];
		}
		k++;
	}
	//*****回代*****
	for(i=n-1;i>=1;i--)
	{
		sum=0;
		for(j=i+1;j<=n;j++)
		{
			sum=sum+a[i][j]*b[j];
		}
		b[i]=b[i]-sum;
	}
	cout<<endl<<endl<<"最后結(jié)果:"<<endl;
	for(int loop=1;loop<=n;loop++)
	{
		c[loop] = b[loop];
		cout<<"x"<<loop<<"="<<c[loop]<<endl;
	}
	return 0;
}