Lagrange Interpolation
(Lagrange_Interpolation.py)

View this file on GitHub
Last update: 2023-08-06 22:58:46+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(L,X,P):
    """ 法が P の下でのLagrange 補間を行い, x=X での値を返す.

    [Input]
    L: [(x_0,y_0), ..., (x_N, y_N)]: F(x_i)=y_i (mod P)
    X: F(X) を返す.
    P: 法

    [Output]
    F(X)

    [Complexity]
    O(N^2+log P)
    """

    N=len(L)-1

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

    X%=P
    Y=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
        Y+=y[i]*c; Y%=P
    return Y

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(L,a,b,X,P):
    """ 法が P の下でのLagrange 補間を行い, x=X での値を返す. ただし, x_i=ai+b

    [Input]
    L: [y_0, ..., y_n]: F(x_i)=y_i (mod P)
    X: F(X) を返す.
    P: 法

    [Output]
    F(X)

    [Complexity]
    O(N+log P)
    """

    d=len(L)-1

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

    Right=[1]*(d+1)
    for i in range(d,-1,-1):
        if i<d:
            Right[i]=(Right[i+1]*(X-(a*i+b)))%P
        else:
            Right[i]=(X-(a*i+b))%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):
        V_inv=(Fact_inv[i]*Fact_inv[d-i])%P
        if i==0:
            S=(Right[i+1]*V_inv)%P
        elif i==d:
            S=(Left[i-1]*V_inv)%P
        else:
            u=(Left[i-1]*Right[i+1])%P
            S=(u*V_inv)%P

        M=L[i]*S%P
        Y=(Y+coef*M)%P
        coef=-coef
    return Y

Traceback (most recent call last):
  File "/opt/hostedtoolcache/Python/3.13.3/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.3/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/python.py", line 96, in bundle
    raise NotImplementedError
NotImplementedError

Lagrange Interpolation (Lagrange_Interpolation.py)

Outline

Theory

Code

Lagrange Interpolation
(Lagrange_Interpolation.py)