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 (xi, yi).

Algorithm

Given the n points (x0, y0), ..., (xn-1, yn-1), compute the Lagrange polynomial ${\displaystyle 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})}}}$. Note that the ith term in the sum, ${\displaystyle y_{i}{\frac {\prod _{j\neq i}(x-x_{j})}{\prod _{j\neq i}(x_{i}-x_{j})}}}$ 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

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 := ${\displaystyle y_{i}/\prod _{j\neq i}(x_{i}-x_{j})}$
compute term := tmp*${\displaystyle \left(\prod _{j\neq i}(x-x_{j})\right)}$

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

#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