Algorithm Implementation/Mathematics/Polynomial interpolation

Lagrange polynomial 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 . 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 j≠i, 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 := 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 #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