Algorithm Implementation/Mathematics/Polynomial interpolation

Lagrange interpolation is an algorithm which returns the polynomial of minimum degree which passes through a given set of points (x_i, y_i).

Algorithm

Given the Template:Math points Template:Math, compute the Lagrange polynomial $p(x)=\sum _{i=0,...,n-1}{y}_i\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}$. Note that the Template:Math term in the sum, $y_i\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}$ is constructed so that when Template:Math is substituted for Template:Math to have a value of zero whenever Template:Math, and a value of Template:Math whenever Template:Math = Template:Math. The resulting Lagrange polynomial is the sum of these terms, so has a value of Template:Math = Template:Math = Template:Math for each of the specified points Template:Math.

In both the pseudocode and each implementation below, the polynomial Template:Math = Template:Math is represented as an array of it's coefficients, Template:Math.

Pseudocode

algorithm lagrange-interpolate is
    input: points Template:Math 
    output: Polynomial p such that Template:Math passes through the input points and is of minimal degree

    for each point (xi, yi) do
        compute tmp := $y_i/\prod _{j\ne i}(x_i-x_j)$
        compute term := tmp*$(\prod _{j\ne i}(x-x_j))$          

    return p, the sum of the values of term

In sample implementations below, the polynomial Template:Math = Template:Math is represented as an array of it's coefficients, Template:Math.

While the code is written to expect points taken from the real numbers (aka floating point), returning a polynomial with coefficients in the reals, this basic algorithm can be adapted to work with inputs and polynomial coefficients from any field, such as the complex numbers, integers mod a prime or finite fields.

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