UVALive 6185 高斯消元浮点型模版

题意：有一个未知的多项式，给你一个d，这个多项式最高次项最大就为d，d<=5，然后给出从0到d+2这个多项式的值，其中有一个是错误的（误差大于1），其余的误差均小于1e6，问错误的是那个

链接：点击打开链接

代码：///就是网上找的模版。。模版出处

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;

///高斯消元模板
const double eps = 1e-3;
const int Max_M = 20;       ///m个方程,n个变量
const int Max_N = 20;
int m, n;
double Aug[Max_M][Max_N+1]; ///增广矩阵
bool free_x[Max_N];         ///判断是否是不确定的变元
double x[Max_N];            ///解集

int sign(double x){ return (x>eps) - (x<-eps);}

/**
返回值:
-1 无解
0 有且仅有一个解
>=1 有多个解,根据free_x判断哪些是不确定的解
*/
int Gauss()
{
    int i,j;
    int 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行以下全为0(包括row行)
        {
            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;
}
///模板结束

double fx[20];///记录f（x）正确的值
int main()
{
    int d,i;
    while(cin>>d)
    {
        if(d==0)break;
        for(int i=0;i<=d+2;i++)
            cin>>fx[i];
        int no;
        for(no = 0;no <= d + 2;no ++){///no和no+1这次不使用
            memset(Aug,0.0,sizeof(Aug));
            memset(x,0.0,sizeof(x));
            memset(free_x,1,sizeof(free_x));
            int now = 0;
            for(i = 0;i <= d + 2;i ++){
                if(i==no)continue;
                Aug[now][0] = 1;
                for(int j = 1;j <= d;j ++)
                    Aug[now][j] = Aug[now][j-1] * i;
                Aug[now][d+1] = fx[i];
                now++;
            }
            m = d+2;
            n = d+1;
            int ans = Gauss();//cout<<x[0]<<endl;
            //cout<<ans<<endl;
            if(ans == 0)break;
        }
        cout<<no<<endl;
    }
    return 0;
}