Algorithm Implementation/Mathematics/Polynomial interpolation

Wikipedia has related information at Lagrange polynomial.

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 $n$ points $(x 0, y 0), ..., (x n -1, y n -1)$ , compute the Lagrange polynomial $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 $i th$ term in the sum, $y_{i}{\frac {\prod _{j\neq i}(x-x_{j})}{\prod _{j\neq i}(x_{i}-x_{j})}}$ is constructed so that when $x j$ is substituted for $x$ to have a value of zero whenever $j \neq i$ , and a value of $y j$ whenever $j$ = $i$ . The resulting Lagrange polynomial is the sum of these terms, so has a value of $p (x j)$ = $0 + 0 + ... + y j + ... + 0$ = $y j$ for each of the specified points $(x j, y j)$ .

In both the pseudocode and each implementation below, the polynomial $p (x)$ = $a 0 + a 1 x + a 2 x 2 + ... + a n -1 x n -1$ is represented as an array of it's coefficients, $(a 0, a 1, a 2, ..., a n -1)$ .

Pseudocode

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

In sample implementations below, the polynomial $p (x)$ = $a 0 + a 1 x + a 2 x 2 + ... + a n -1 x n -1$ is represented as an array of it's coefficients, $(a 0, a 1, a 2, ..., a n -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