演算法實現/數學/多項式插值

拉格朗日插值是一種演算法，它返回透過給定點集 (x_i, y_i) 的最小次數多項式。

演算法

給定 $n$ 個點 $(x 0, y 0), ..., (x n -1, y n -1)$ ，計算拉格朗日插值多項式 $p(x)=\sum _{i=0,...,n-1}y_{i}{\frac {\prod _{j\neq i}(x-x_{j})}{\prod _{j\neq i}(x_{i}-x_{j})}}$ 。請注意，總和中的 $i th$ 項， $y_{i}{\frac {\prod _{j\neq i}(x-x_{j})}{\prod _{j\neq i}(x_{i}-x_{j})}}$ 被構造為，當 $x j$ 替換為 $x$ 時，只要 $j \neq i$ ，其值為零，而只要 $j$ = $i$ ，其值為 $y j$ 。由此產生的拉格朗日插值多項式是這些項的總和，因此對於每個指定點 $(x j, y j)$ ，其值均為 $p (x j)$ = $0 + 0 + ... + y j + ... + 0$ = $y j$ 。

在下面的虛擬碼和每個實現中，多項式 $p (x)$ = $a 0 + a 1 x + a 2 x 2 + ... + a n -1 x n -1$ 表示為其係數陣列， $(a 0, a 1, a 2, ..., a n -1)$ 。

虛擬碼

algorithm lagrange-interpolate is
    input: points  $(x 0, y 0), ..., (x n -1, y n -1)$  
    output: Polynomial p such that  $p (x)$  passes through the input points and is of minimal degree

    for each point (x_i, y_i) do
        compute tmp :=  $y_{i}/\prod _{j\neq i}(x_{i}-x_{j})$ 
        compute term := tmp* $\left(\prod _{j\neq i}(x-x_{j})\right)$           

    return p, the sum of the values of term

在下面的示例實現中，多項式 $p (x)$ = $a 0 + a 1 x + a 2 x 2 + ... + a n -1 x n -1$ 表示為其係數陣列， $(a 0, a 1, a 2, ..., a n -1)$ 。

雖然程式碼被編寫為期望從實數（即浮點數）中獲取點，並返回一個具有實數係數的多項式，但這個基本演算法可以被調整為使用來自任何域（例如複數、模素數的整數或有限域）的輸入和多項式係數。

C語言

#include <stdio.h>
#include <stdlib.h>

// input: numpts, xval, yval
// output: thepoly
void interpolate(int numpts, const float xval[restrict numpts], const float yval[restrict numpts],
    float thepoly[numpts])
{
    float theterm[numpts];
    float prod;
    int i, j, k;
    for (i = 0; i < numpts; i++)
        thepoly[i] = 0.0;
    for (i = 0; i < numpts; i++) {
        prod = 1.0;
        for (j = 0; j < numpts; j++) {
            theterm[j] = 0.0;
        };
        // Compute Prod_{j != i} (x_i - x_j)
        for (j = 0; j < numpts; j++) {
            if (i == j)
                continue;
            prod *= (xval[i] - xval[j]);
        };
        // Compute y_i/Prod_{j != i} (x_i - x_j)
        prod = yval[i] / prod;
        theterm[0] = prod;
        // Compute theterm := prod*Prod_{j != i} (x - x_j)
        for (j = 0; j < numpts; j++) {
            if (i == j)
                continue;
            for (k = numpts - 1; k > 0; k--) {
                theterm[k] += theterm[k - 1];
                theterm[k - 1] *= (-xval[j]);
            };
        };
        // thepoly += theterm (as coeff vectors)
        for (j = 0; j < numpts; j++) {
            thepoly[j] += theterm[j];
        };
    };
}

Python

from typing import Tuple, List

def interpolate(inpts: List[Tuple[float, float]]) -> List[float]:
    n = len(inpts)
    thepoly = n * [0.0]
    for i in range(n):
        prod = 1.0
        # Compute Prod_{j != i} (x_i - x_j)
        for j in (j for j in range(n) if (j != i)):
            prod *= (inpts[i][0] - inpts[j][0])
        # Compute y_i/Prod_{j != i} (x_i - x_j)
        prod = inpts[i][1] / prod
        theterm = [prod] + (n - 1) * [0]
        # Compute theterm := prod*Prod_{j != i} (x - x_j)
        for j in (j for j in range(n) if (j != i)):
            for k in range(n - 1, 0, -1):
                theterm[k] += theterm[k - 1]
                theterm[k - 1] *= (-inpts[j][0])
        # thepoly += theterm
        for j in range(n):
            thepoly[j] += theterm[j]
    return thepoly