zoj 3645 高斯消元

题意是：给12个方程，形如 （a0-x0）^2 +(a1-x1)^2+ ·········+ （a10-x10）^2=d^2;

利用嘴一个方程和上面11个方程相减，编程一次方程组 yoga高斯消元法解方程；

#include<cstdio>
#include<cmath>
#include<cstring>
#include<queue>
#include<iostream>
using namespace std;
#define FF  freopen("Input.txt","r",stdin)
#define mem(x,y) memset(x,y,sizeof(x))
const double eps=1e-12;
int n,m; //m个方程，n个变量
double Aug[15][15]; //增广矩阵
bool free_x[15]; //判断是否是不确定的变元
double x[15];    //解集
int sign(double x) { return (x>eps)- (x<-eps) ;}
/*
返回值：-1 无解，
        0 有且紧一组解。
        >=1 有多个解。
*/
int gauss()
{
    int i,j,row,col,max_r;
    for(row=0,col=0;row<m&&col<n;row++,col++)
    {
        max_r=row;
        for(i=row+1;i<m;i++)
        {
            if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0)
               max_r=i;
        }
        if(max_r!=row)
        {
            for(j=row;j<n+1;j++)
              swap(Aug[max_r][j],Aug[row][j]);
        }
        if(sign(Aug[row][col])==0) { row--;continue;}
        for(i=row+1;i<m;i++)
        {
            if(sign(Aug[i][col])==0) continue;
            double ta=Aug[i][col]/Aug[row][col];
            for(j=col;j<n+1;j++)
               Aug[i][j]-=Aug[row][j]*ta;
        }
    }
    for(i = row; i < m; i++) ///col=n存在0...0,a的情况,无解
    {
        if(sign(Aug[i][col]))
            return -1;
    }
    if(row < n) //存在0...0,0的情况,有多个解,自由变元个数为n-row个
    {
        for(i = row-1; i >=0; i--)
        {
            int free_num = 0;   //自由变元的个数
            int free_index;     //自由变元的序号
            for(j = 0; j < n; j++)
            {
                if(sign(Aug[i][j])!=0 && free_x[j])
                    free_num++,free_index=j;
            }
            if(free_num > 1) continue; //该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元
            //只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的
            double tmp = Aug[i][n];
            for(j = 0; j < n; j++)
            {
                if(sign(Aug[i][j])!=0 && j!=free_index)
                    tmp -= Aug[i][j]*x[j];
            }
            x[free_index] = tmp/Aug[i][free_index];
            free_x[free_index] = false;
        }
        return n-row;
    }//有且仅有一个解,严格的上三角矩阵(n==m)
    for(i = n-1; i >= 0; i--)
    {
        double tmp = Aug[i][n];
        for(j = i+1; j < n; j++)
            if(sign(Aug[i][j])!=0)
                tmp -= Aug[i][j]*x[j];
        x[i] = tmp/Aug[i][i];
    }
    return 0;
}
int main()
{
    FF;
    int i,j,t;
    double a[12][12];
    scanf("%d",&t);
    while(t--)
    {
        mem(Aug,0.0);
        mem(x,0.0); mem(free_x,true);
        for(i=0;i<12;i++)
          for(j=0;j<12;j++)
            scanf("%lf",&a[i][j]);
        double sum=0;
        for(i=0;i<11;i++) sum+=a[11][i]*a[11][i];
        for(i=0;i<11;i++)
        {
            for(j=0;j<11;j++)
            {
                Aug[i][j]=2*( a[i][j]-a[11][j] );
                Aug[i][11]+=a[i][j]*a[i][j];
            }
            Aug[i][11]+=-a[i][11]*a[i][11]+a[11][11]*a[11][11]-sum;
        }
        m=n=11;
        gauss();
        for(i=0;i<n;i++)
           printf("%.2lf%c",x[i],i==n-1?'\n':' ');
    }
    return 0;
}