ラグランジュ補完 (Lagrange Polynomial)

概要

\(N\) 個の点 \((x_0, y_0), (x_1, y_1), ..., (x_{N-1}, y_{N-1})\) を通る多項式を求める。

\(\displaystyle g_i(x) = \frac{(x - x_0) \ldots (x - x_{i-1})(x - x_{i+1}) \ldots (x - x_{N-1})}{(x_i - x_0) \ldots (x_i - x_{i-1})(x_i - x_{i+1}) \ldots (x_i - x_{N-1})}\)

として、

\(\displaystyle f(x) = \sum_{i=0}^{N-1} y_i g_i(x)\)

とする。

すると、 \(g_i(x_j) = 0\) \((i \not = j)\), \(g_i(x_i) = 1\) より、 \(f(x_i) = y_i g_i(x_i) = y_i\) を満たす。

実装

# ラグランジュ補完の多項式の係数計算
# N 個の点 (x_0, y_0), (x_1, y_1), ..., (x_{N-1}, y_{N-1}) について
# X = [x_0, x_1, ..., x_{N-1}]
# Y = [y_0, y_1, ..., y_{N-1}]
def build(X, Y):
    A = []
    for x in X:
        res = 1
        for xi in X:
            if x == xi:
                continue
            res *= x - xi
        A.append(Y[x] / res)
    return A
 
# 係数の値から f(k) の値を計算
# (全ての i について k != x_i とする)
def calc(X, A, k):
    # assert i not in X
    base = 1
    for xi in X:
        base *= k - xi
    return sum(base / (k - x) * a for x, a in zip(X, A))

Verified

  • AOJ: "1328: Find the Outlier": source (Python3, 0.36sec)

参考

戻る