Algorithm Implementation/Mathematics/Polynomial interpolation

From Wikibooks, open books for an open world
Jump to navigation Jump to search

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

Algorithm[edit | edit source]

Given the n points (x0, y0), ..., (xn-1, yn-1), compute the Lagrange polynomial . Note that the ith term in the sum, is constructed so that when xj is substituted for x to have a value of zero whenever ji, and a value of yj whenever j = i. The resulting Lagrange polynomial is the sum of these terms, so has a value of p(xj) = 0 + 0 + ... + yj + ... + 0 = yj for each of the specified points (xj, yj).

In both the pseudocode and each implementation below, the polynomial p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 is represented as an array of it's coefficients, (a0, a1, a2, ..., an-1).

Pseudocode[edit | edit source]

algorithm lagrange-interpolate is
    input: points (x0, y0), ..., (xn-1, yn-1) 
    output: Polynomial p such that p(x) passes through the input points and is of minimal degree

    for each point (xi, yi) do
        compute tmp := 
        compute term := tmp*          

    return p, the sum of the values of term

In sample implementations below, the polynomial p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 is represented as an array of it's coefficients, (a0, a1, a2, ..., an-1).

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[edit | edit source]

#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[edit | edit source]

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