BZOJ 4689 Find the Outlier 高斯消元

Description

Abacus教授刚刚完成了一个制作数表的计算引擎的设计。它被设计用于同时计算一个多项式在许多点的取值。例如

对于多项式 f(x)=x^2+2x+1 ，一种可能的计算结果是 f(0)=1,f(1)=4,f(2)=9.f(3)=16,f(4)=25 。不幸的是，引

擎存在一个故障使得计算出的值总有一个是错的，例如对于上述多项式，它可能输出 1,4,12,16,25 而不是 1,4,9

,16,25 。请你帮教授找出发生故障的是哪个点值。

Input

输入包含多组测试数据。

每组数据第一行包含一个正整数 d 表示多项式的度数，即多项式最高次项的项数，保证 d≤5 。

接下来 d+3 行，每行一个实数，第 i 行表示输出的 f(i)  的值，保证-100.0≤f(i)≤100.0 。

你可以认为恰好只有一个点值出故障，且与实际值的误差超过 1.0 。

由于不可避免的误差，其他数字与精确值的误差不超过10^(-6)  。

输入以一个零作为结束。

Output

对于每组数据，输出一个非负整数 i 表示 f(i)  的值发生故障。

Sample Input

2
1.0
4.0
12.0
16.0
25.0
1
-30.5893962764
5.76397083962
39.3853798058
74.3727663177
4
42.4715310246
79.5420238202
28.0282396675
-30.3627807522
-49.8363481393
-25.5101480106
7.58575761381
5
-21.9161699038
-48.469304271
-24.3188578417
-2.35085940324
-9.70239202086
-47.2709510623
-93.5066246072
-82.5073836498
0

Sample Output

2
1
1
6

HINT

传送门

……什么鬼。。

D如此之小。。

给出D+3个方程，有一个错掉的，问错了哪个。

那当然是暴力枚举了。。然后剩下(D+2)个方程，

由于只有(D+1)个元，所以会有一个多余。

那么只要任意取(D+1)个方程，解出它们的解，然后判断剩下的方程满不满足即可。

假如满足，根据题意，当前枚举到的不合法方程一定就是答案了。

好久没打高斯消元都快不会了= =

只好手推一遍2333

#include<bits/stdc++.h>
using namespace std;
const double 
	eps=1e-3;
const int 
	N=15;
int D,ksm[9][6];
double a[N][N],tG[N];
bool gauss(int tt){
	double x[N];
	for (int k=1;k<=D+1;k++){
		for (int i=1;i<=D+1;i++) x[i]=a[i][k]/a[k][k];
		x[k]=0.0;
		for (int i=k+1;i<=D+1;i++){
			tG[i]-=tG[k]*x[i];
			for (int j=1;j<=D+1;j++) a[i][j]-=a[k][j]*x[i];
		}
	}
	double ANS[N];
	ANS[D+1]=tG[D+1]/a[D+1][D+1];
	for (int i=D;i;i--){
		double tmp=0.0;
		for (int j=i+1;j<=D+1;j++)
			tmp+=ANS[j]*a[i][j];
		ANS[i]=(tG[i]-tmp)/a[i][i];
	}

	double tmp=0.0;
	for (int i=1;i<=D+1;i++) tmp+=ANS[i]*a[D+2][i];
	if (fabs(tmp-tG[D+2])<=eps) return 1;
	return 0;
}
void solve(){
	double G[N];
	for (int i=1;i<=D+3;i++) scanf("%lf",&G[i]);
	for (int i=1;i<=D+3;i++){
		int n=0;
		for (int j=1;j<=D+3;j++)
			if (i!=j){
				tG[++n]=G[j];
				for (int k=D;~k;k--) a[n][k+1]=1.0*ksm[j][k];
			}
		if (gauss(i)){printf("%d\n",i-1);return;}
	}
}
int main(){
	for (int i=1;i<9;i++){
		ksm[i][0]=1;
		for (int j=1;j<6;j++) ksm[i][j]=ksm[i][j-1]*i;
	}
	while (scanf("%d",&D),D) solve();
	return 0;
}