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 {$(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&ne;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 (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_{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.