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

void mt(float *a,float *b,int row,int column)//矩陣求轉(zhuǎn)置函數(shù)
{
	int i,j;//a是轉(zhuǎn)置前矩陣 b是轉(zhuǎn)置后的矩陣 row為矩陣的行數(shù) column是矩陣的列數(shù)
    float temp;
    for(i=0;i<row;i++)
		for(j=0;j<column;j++)
		{
			temp=*(a+i*column+j);
			*(b+j*row+i)=temp;
		}
}

void subb(float *a,float *b,int n,float *c)//矩陣相減函數(shù)
{                                          //c中存放差矩陣
   int i;
   for(i=0;i<n;i++)
   {
	   *(c+i)=*(a+i)-*(b+i);
   }
}

void madd(float *a,float *b,int n,float *c)//矩陣相加函數(shù)
{                                          //c中存放和矩陣
   int i;
   for(i=0;i<n;i++)
   {
	   *(c+i)=*(a+i)+*(b+i);//fabs是浮點(diǎn)數(shù)絕對值函數(shù)
   }
}

void kmul(float *a,int n,float c,float *b)//矩陣數(shù)乘函數(shù)
{
   int i;
   for(i=0;i<n;i++)
   {
	   *(b+i)=(*(a+i))*c;
   }
}

void bmul(float *a,float *b,int m,int n,int k,float *result)//矩陣相乘矩陣
{ 
	int i,j,l,u;              //矩陣result=矩陣a*矩陣b
    for (i=0; i<=m-1; i++)    //m是矩陣a的行數(shù)，n是矩陣a的列數(shù) k是矩陣b的列數(shù)
		for (j=0; j<=k-1; j++)
		{
			u=i*k+j;
			result[u]=0.0;
			for (l=0; l<=n-1; l++)
				result[u]=result[u]+a[i*n+l]*b[l*k+j];
		}
    return;
}

int inv(float *a,int n) //矩陣求逆函數(shù) n是矩陣a的階數(shù)
{
	int *is,*js,i,j,k,l,u,v;
    float d,p;
    
    is=new int [n];
    js=new int [n];
  
    for (k=0; k<=n-1; k++)
	{
		d=0.0;
        for (i=k; i<=n-1; i++)
			for (j=k; j<=n-1; j++)
			{
				l=i*n+j;
				p=fabs(a[l]);
				if (p>d) { d=p; is[k]=i; js[k]=j; }
			}
        if (d+1.0==1.0)
		{
			delete is;
			delete js;
			cout<<"error matrix can not inv"<<endl;
            return(0);
		}
        if (is[k]!=k)
			for (j=0; j<=n-1; j++)
			{
				u=k*n+j; v=is[k]*n+j;
				p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (js[k]!=k)
			for (i=0; i<=n-1; i++)
            {
				u=i*n+k; v=i*n+js[k];
				p=a[u]; a[u]=a[v]; a[v]=p;
            }
        l=k*n+k;
        a[l]=1.0/a[l];
        for (j=0; j<=n-1; j++)
		  if (j!=k)
		  { u=k*n+j; a[u]=a[u]*a[l]; }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                {
				  u=i*n+j;
                  a[u]=a[u]-a[i*n+k]*a[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
		  { u=i*n+k; a[u]=-a[u]*a[l];}
	}
    for (k=n-1; k>=0; k--)
	{
		if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            {
			  u=k*n+j; v=js[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            {
			  u=i*n+k; v=i*n+is[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
	}
    delete is; delete js;
    return(1);
}


void main()
{
	float x1[4][2]={{0,0},{2,0},{2,2},{0,2}};
	float x2[4][2]={{4,4},{6,4},{6,6},{4,6}};
	float x[2];
	float xr[2][2];
	float *temp;
	float p1=0.5;
	float p2=0.5;
	float m1[2]={0,0};
	float m2[2]={0,0};
	float m1r[2][2],m2r[2][2];
	float m[2];
	float c1[2][2]={{0,0},{0,0}};
	float c2[2][2]={{0,0},{0,0}};
	float mc1[2],mc2[2];
	float result[2];
	float con1,con2;
	float N1=4;
	float N2=4;
	int i,j,k=0;

	for(i=0;i<N1;i++)
	{
		temp=&x1[i][0];
		x[0]=*temp;
		x[1]=*(temp+1);
		madd(&m1[0],&x[0],2,&m1[0]);
		bmul(&x[0],&x[0],2,1,2,&xr[0][0]);
		madd(&c1[0][0],&xr[0][0],4,&c1[0][0]);
	}
	kmul(&m1[0],2,1/N1,&m1[0]);
	kmul(&c1[0][0],4,1/N1,&c1[0][0]);
	bmul(&m1[0],&m1[0],2,1,2,&m1r[0][0]);
	subb(&c1[0][0],&m1r[0][0],4,&c1[0][0]);

	for(i=0;i<N2;i++)
	{
		temp=&x2[i][0];
		x[0]=*temp;
		x[1]=*(temp+1);
		madd(&m2[0],&x[0],2,&m2[0]);
		bmul(&x[0],&x[0],2,1,2,&xr[0][0]);
		madd(&c2[0][0],&xr[0][0],4,&c2[0][0]);
	}
	kmul(&m2[0],2,1/N2,&m2[0]);
	kmul(&c2[0][0],4,1/N2,&c2[0][0]);
	bmul(&m2[0],&m2[0],2,1,2,&m2r[0][0]);
	subb(&c2[0][0],&m2r[0][0],4,&c2[0][0]);

	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			if(c1[i][j]==c2[i][j]) k++;

	if(k==4)
	{
		subb(&m1[0],&m2[0],2,&m[0]);
		inv(&c1[0][0],2);
		bmul(&m[0],&c1[0][0],1,2,2,&result[0]);
		bmul(&m1[0],&c1[0][0],1,2,2,&mc1[0]);
		bmul(&mc1[0],&m1[0],1,2,1,&con1);
		con1=con1/2;
		bmul(&m2[0],&c2[0][0],1,2,2,&mc2[0]);
		bmul(&mc2[0],&m2[0],1,2,1,&con2);
		con2=con2/2;

		cout<<"the solution is:"<<endl;
		cout<<result[0]<<" "<<result[1]<<" "<<con2-con1<<endl;
	}
	else cout<<"error!"<<endl;
}