AOJ 1328 Find the Outlier

問題

xのd次式f(x)がある。
f(x)に、x = 0, 1, 2, ..., d + 2を代入した結果v[0], v[1], ..., v[d + 2]が与えられる。

このうち、どれか1つは間違った数字になっている。
このとき、その数字を見つけよ。

制約条件

d≦5

v[i]

≦100

方針

2つを取り除いて連立式を作り、f(x)の形を具体的に求める。
このf(x)を使って逆にv[0], ..., v[d + 2]を計算し、外れている値がちょうど1つであるかどうかを見る。

EPSを適当にすると通らないので、10^-9から、答えがユニークになるまで
EPSを10倍しながら同じ解法を繰り返したら通った。

ソースコード

連立方程式の解はspaghetti sourceより。

int n;
double v[10];

vector<double> solve2(int i1, int i2){
	matrix M(n + 1, array(n + 1));
	array b(n + 1);
	int r = 0;
	rep(i, n + 3) if(i != i1 && i != i2){
		rep(j, n + 1) M[r][j] = j ? pow(i * 1.0, j) : 1.0;
		b[r++] = v[i];
	}
	LUinfo data = LU_decomposition(M);
	b = LU_backsubstitution(data, b);
	
	vector<double> res;
	rep(i, n + 3){
		double tmp = 0;
		rep(j, n + 1) tmp += (j ? pow(i * 1.0, j) : 1.0) * b[j];
		res.pb(tmp);
	}
	return res;
}
double EPS;
double solve(int idx){
	double err;
	rep(it, n + 3) if(it != idx){
		vector<double> res = solve2(idx, it);
		if(abs(res[it] - v[it]) > EPS) return 1e12;
		err = abs(res[it] - v[idx]);
	}
	return err;
}
int main(){
	while(cin >> n, n){
		rep(i, n + 3) cin >> v[i];
		for(EPS = 1e-9; ; EPS *= 10){
			double mnerr = 1e12; int mni = -1;
			rep(i, n + 3){
				double err = solve(i);
				if(err < mnerr) mnerr = err, mni = i;
			}
			if(mni >= 0){
				cout << mni << endl;
				break;
			}
		}
	}
	return 0;
}