AOJ 1328 Find the Outlier (ICPC アジア 2012 D) (450 点)

連立方程式の練習。多項式補間でも解ける。誤差がきつい。

問題へのリンク

問題概要

次多項式 がある (係数はわからない)。 個の点 の値がわかっている...はずだったが、そのうちの 1 個については間違った値が出力された状態で入力が与えられる。どれが間違っているかを特定せよ。

なお、間違っている値は、正解と 1.0 以上の差があるとしてよい。答えは までの相対誤差または絶対誤差が許容される。

制約

• 
• 

考えたこと

個の点を与えれば、多項式は一意に決まる。具体的には多項式を

• 

としたときに、 を代入したときの値 がわかっていれば

• 

という についての方程式が立てられる。これが 本あれば解が一意に求められる。

実際は 個の値が与えられているので、そのうちの 個の値を用いて連立一次方程式を解くことを繰り返す。そうして残りの 2 個の値のうち、1 個が正解で 1 個が間違っているならば、その間違っている値が答えであるとわかる。

なお、間違っている値を 個の値に含めてしまったときは、必ず 2 個とも値が外れることがわかる。

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

// Gauss Jordan
using D = long double;
const D EPS = 1e-5;
template<class T> struct Matrix {
    vector<vector<T> > val;
    Matrix(int n, int m, T x = 0) : val(n, vector<T>(m, x)) {}
    void init(int n, int m, T x = 0) {val.assign(n, vector<T>(m, x));}
    size_t size() const {return val.size();}
    inline vector<T>& operator [] (int i) {return val[i];}
};

template<class T> int GaussJordan(Matrix<T> &A, bool is_extended = false) {
    int m = A.size(), n = A[0].size();
    int rank = 0;
    for (int col = 0; col < n; ++col) {
        if (is_extended && col == n-1) break;
        int pivot = -1;
        T ma = EPS;
        for (int row = rank; row < m; ++row) {
            if (abs(A[row][col]) > ma) {
                ma = abs(A[row][col]);
                pivot = row;
            }
        }
        if (pivot == -1) continue;
        swap(A[pivot], A[rank]);
        auto fac = A[rank][col];
        for (int col2 = 0; col2 < n; ++col2) A[rank][col2] /= fac;
        for (int row = 0; row < m; ++row) {
            if (row != rank && abs(A[row][col]) > EPS) {
                auto fac = A[row][col];
                for (int col2 = 0; col2 < n; ++col2) {
                    A[row][col2] -= A[rank][col2] * fac;
                }
            }
        }
        ++rank;
    }
    return rank;
}

template<class T> vector<T> linear_equation(Matrix<T> A, vector<T> b) {
    int m = A.size(), n = A[0].size();
    Matrix<T> M(m, n + 1);
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) M[i][j] = A[i][j];
        M[i][n] = b[i];
    }
    int rank = GaussJordan(M, true);
    vector<T> res;
    for (int row = rank; row < m; ++row) if (abs(M[row][n]) > EPS) return res;
    res.assign(n, 0);
    for (int i = 0; i < rank; ++i) res[i] = M[i][n];
    return res;
}




D pow(D a, int n) {
    D res = 1.0;
    for (int i = 0; i < n; ++i) res *= a;
    return res;
}

D func(const vector<D> &coef, int i) {
    D res = 0.0;
    for (int p = 0; p < coef.size(); ++p) res += coef[p] * pow(i, p);
    return res;
}

int main() {
    int d;
    while (cin >> d, d) {
        vector<D> v(d+3);
        for (int i = 0; i < d+3; ++i) cin >> v[i];

        int res = 0;
        for (int i = 0; i < d+3; ++i) {
            for (int j = i+1; j < d+3; ++j) {
                Matrix<D> A(d+1, d+1);
                vector<D> b(d+1);
                for (int k = 0, iter = 0; k < d+3; ++k) {
                    if (k == i || k == j) continue;
                    for (int p = 0; p < d+1; ++p) {
                        A[iter][p] = pow(k, p);
                        b[iter] = v[k];
                    }
                    ++iter;
                }
                vector<D> ans = linear_equation(A, b);
                if (ans.empty()) continue;
                D vi = func(ans, i), vj = func(ans, j);
                int num = 0;
                if (fabs(vi - v[i]) > EPS) res = i, ++num;
                if (fabs(vj - v[j]) > EPS) res = j, ++num;
                if (num == 1) goto end;
            }
        }
    end:
        cout << res << endl;
    }
}