Lagrange Interpolation
(Lagrange_Interpolation.py)

• View this file on GitHub
• Last update: 2025-07-12 16:16:16+09:00

Outline

$(N+1)$ 個の条件 $y_i=f(x_i)$ $(0 \leq i \leq N)$ を満たす高々 $N$ 次多項式 $f(x)$ や, $f$ のある点での値 $f(X)$ を求める.

Theory

$x_0, \dots, x_N$ は全て異なるとする. このとき, $(N+1)$ 個の条件式

\[y_i=f(x_i) \quad (0 \leq i\leq N)\]

を満たす高々 $N$ 次の多項式を求める. $x_0, \dots, x_N$ が全て異なるという条件下では $f$ は必ず一意に存在する.

$i=0,1,2, \dots, N$ に対して, Lagrange 基底多項式 $\ell_i(x)$ を

\[\ell_i(x):=\prod_{\substack{0 \leq k \leq N \\ k \neq i}} \dfrac{x-x_k}{x_i-x_k}\]

とする. このとき, $j=0,1, \dots, N$ に対して,

\[\ell_i(x_j)=[i=j]\]

が成り立つ.

よって,

\[f:=\sum_{i=0}^N y_i \ell_i\]

がその求めるべき多項式になる.

ここで, $x_0, \dots, x_N$ が等差数列であったとする. つまり, $x_i=ai+b$ であるとする. このとき,

\[x_i-x_k=(ai+b)-(ak+b)=a(i-k)\]

が成り立つから, Lagrange 基底多項式は

\[\begin{align*} \ell_i(x) &=\prod_{\substack{0 \leq k \leq N \\ k \neq i}} \dfrac{x-x_k}{x_i-x_k} \\ &=\prod_{\substack{0 \leq k \leq N \\ k \neq i}} \dfrac{x-(ai+b)}{a(i-k)} \\ &=\dfrac{1}{(-a)^N} \cdot \dfrac{ \mathrm{Left}_{i-1}(x) \cdot \mathrm{Right}_{i+1}(x)}{(-1)^{i} i! (N-i)!} \\ \end{align*}\]

である. ただし,

\[\mathrm{Left}_i(x):=\prod_{k=0}^i (x-(ak+b)), \quad \mathrm{Right}_i(x):=\prod_{k=i}^N (x-(ak+b))\]

である.

よって, 評価点の $x$ 座標が等差数列であるとき, ある点 $X$ での値 $f(X) \pmod{P} $ を $O(N+\log P)$ 時間で求めることができる.

Code

def Lagrange_Interpolation_Point(F: list[tuple[int, int]], X: int, P: int) -> int:
    """ F の情報から, 多項式 f に関する P を法とする Lagrange 補完を行い, f(X) を求める.

    Args:
        F (list[tuple[int, int]]): [(x0, y0), ..., (xN,yN)] の形のリスト, i 番目の要素は, f が f(x[i]) = y[i] であることを意味する.
        X (int): 求める X の値
        P (int): 法

    Returns:
        int: f(X) mod P

    Complexity:
        O(N^2 + log P)
    """

    N = len(F) - 1

    x = [p[0] for p in F]
    y = [p[1] for p in F]

    X %= P
    z = 0
    for i in range(N + 1):
        a = b = 1
        for j in range(N + 1):
            if i == j:
                continue

            a *= X - x[j]; a %= P
            b *= x[i] - x[j]; b %= P

        c = a * pow(b, -1, P) % P
        z += y[i] * c % P

    return z % P

def Lagrange_Interpolation_Polynomial(T,P):
    """ 法が P の下でのLagrange 補間を行い, 多項式の係数リストを返す.

    [Input]
    T: [(x_0,y_0), ..., (x_N, y_N)]: F(x_n)=y_n  (i !=j => x_i != x_j (mod P))
    P: 法

    [Output]
    [[X^0]F, [X^1]F, ..., [X^{N-1}]F ]

    [Complexity]
    O(N^2)

    [Thanks]
    hamayanhamayan
    """

    N=len(T)
    X=[0]*N; Y=[0]*N
    for i in range(N):
        X[i]=T[i][0]
        Y[i]=T[i][1]

    Poly=[0]*(N+1); Poly[0]=1
    for x,y in zip(X,Y):
        tmp=[0]*(N+1)
        for i in range(N):
            tmp[i+1]=Poly[i]
        for i in range(N):
            tmp[i]=(tmp[i]-x*Poly[i])%P
        Poly=tmp

    res=[0]*N
    for i,(x,y) in enumerate(zip(X,Y)):
        if y==0:
            continue

        Q=1
        for j in range(N):
            if j!=i:
                Q=Q*(x-X[j])%P
        Q=pow(Q, -1, P)

        tmp=Poly.copy()

        for j in  range(N,0,-1):
            res[j-1]=(res[j-1]+(tmp[j]*Q)%P*y)%P
            tmp[j-1]=(tmp[j-1]+tmp[j]*x)%P
    return res

def Lagrange_Interpolation_Point_Arithmetic(a: int, b: int, Y: list[int], X: int, P: int) -> int:
    """ 法が P の下での Lagrange 補間を行い, x = X での値を返す. ただし, x_i = a i + b

    Args:
        a (int): x の値の傾き
        b (int): x の値の切片
        Y (list[int]): f(a * i + b) = Y[i] となる Y
        X (int): 求める X の値
        P (int): 剰余

    Returns:
        int: f(X) mod P

    Complexity:
        O(N + log P)
    """

    d = len(Y) - 1

    X %= P
    left = [1] * (d + 1)
    for i in range(d + 1):
        alpha = X - a * i + b
        if i > 0:
            left[i] = left[i - 1] * alpha % P
        else:
            left[i] = alpha % P

    right = [1] * (d + 1)
    for i in range(d, -1, -1):
        alpha = X - a * i + b
        if i < d:
            right[i] = right[i + 1] * alpha % P
        else:
            right[i] = alpha % P

    fact = 1
    for i in range(1, d + 1):
        fact = (fact * i) % P

    fact_inv = [1] * (d + 1)
    fact_inv[-1] = pow(fact, -1, P)
    for i in range(d - 1, -1, -1):
        fact_inv[i] = fact_inv[i + 1] * (i + 1) % P

    y = 0
    coef = pow(-a, -d, P)

    for i in range(d + 1):
        alpha = fact_inv[i] * fact_inv[d - i] % P
        if i == 0:
            eta = right[i + 1]
        elif i == d:
            eta = left[i - 1]
        else:
            eta = left[i - 1] * right[i + 1] % P

        beta = eta * alpha % P
        gamma = Y[i] * beta % P
        y += coef * gamma % P
        coef = -coef

    return y % P

Traceback (most recent call last):
  File "/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/onlinejudge_verify/documentation/build.py", line 71, in _render_source_code_stat
    bundled_code = language.bundle(stat.path, basedir=basedir, options={'include_paths': [basedir]}).decode()
                   ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/python.py", line 96, in bundle
    raise NotImplementedError
NotImplementedError

Back to top page