User:Rcampbell1/sandbox

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₀, y₀), ..., (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≠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₀ + a₁x + a₂x² + ... + a_n-1x^n-1 is represented as an array of it's coefficients, (a₀, a₁, a₂, ..., a_n-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 := $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₀ + a₁x + a₂x² + ... + a_n-1x^n-1 is represented as an array of it's coefficients, (a₀, a₁, a₂, ..., 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>

void interpolate(int numpts, float* xval, float* yval, float* thepoly){
   float theterm[10];
   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];};
   };
}
         
int main(){
   float xval[10], yval[10];
   float thepoly[10], tmpx, tmpy;
   int numpts, i;
   printf("number of points: "); scanf("%d",&numpts);
   printf("Enter %d points (x,y):\n",numpts);
   for (i=0; i<numpts; i++){
      printf("x y: "); scanf("%f %f",&tmpx,&tmpy);
      xval[i] = tmpx; yval[i] = tmpy;
   };
   interpolate(numpts,xval,yval,thepoly);
   printf("p(x) = (%6.2f)",thepoly[0]);
   if (thepoly[1] != 0.0) {printf(" + (%6.2f)x",thepoly[1]);};
   for (i=2; i<numpts; i++){
      if (thepoly[i] != 0.0){
         printf(" + (%6.2f)x^%d",thepoly[i],i);
      };
   };
   printf("\n");
   exit(0);
}

Python

def interpolate(inpts):
    n = len(inpts)
    thepoly = n*[0]
    for i in range(len(inpts)):
        prod = 1
        # 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)