Lagrange Interpolation

One particularly useful property of polynomials is that a polynomial of degree $m$ is completely determined on $m + 1$ evaluation points, which implies that we can uniquely derive a polynomial of degree $m$ from a set $S$:

$$S = \{(x_0, y_0), (x_1, y_1), . . . , (x_m, y_m) | x_i \neq x_j \text{ for all indices i and j}\}$$

Polynomials therefore have the property that $m + 1$ pairs of points $(x_i, y_i) \text{ for } x_i \neq x_j$ are enough to determine the set of pairs $(x, P(x)) \text{ for all } x ∈ R$.

If the coefficients of the polynomial we want to find have a notion of multiplicative inverse, it is always possible to find such a polynomial using a method called Lagrange Interpolation, which works as follows.

The set of polynomials $l_j(x)$is called the Lagrange basis. Each of the polynomials has degree $m$ and take values $l_j(x_i) = 0$ if $i \neq j$ and $l_j(x_j) = 1$. Using the Kronecker delta this can be written as $l_j(x_i) = \delta_{ij}$.

Notice that the numerator $\prod_{i \neq j}(x - x_i)$ has m roots at the nodes $\{x_i\}_{i \neq j}$ while the denominator $\prod_{i \neq j}(x_j-x_i)$ scales the resulting polynomial so that $l_j(x_j) = 1$.

The naive implemetation of Lagrange Interpolation can be found below:

def lagrange_polynomial(evals):
    Qx = QQ[x]
    m = len(evals)
    P = Qx(0)
    for j in range(0, m):
        l_j = Qx(1)
        for i in range(0, m):
            if i == j:
                continue
            l_j = l_j * Qx(x - evals[i][0]) / (evals[j][0] - evals[i][0])
        P = P + l_j * evals[j][1]
    return P   

Sage example:

Qx = QQ[x]
Qx.lagrange_polynomial([(0, 4), (-2, 1), (2, 3)])