# A quick barycentric evaluation tutorial

One powerful technique when working with polynomials is taking a set of evaluations of that polynomial and using that _directly_ to compute an evaluation at a different point. For example, if $P(x)$ is a degree-100 polynomial, you can use the evaluations $P(0), P(1) ... P(100)$ to directly compute $P(101)$, or $P(1247130)$, in $O(N)$ time, without ever reconstructing the polynomial. This post describes how this is done.

_See also, this earlier and more detailed paper by Oxford researchers on the topic: https://people.maths.ox.ac.uk/trefethen/barycentric.pdf_

## General technique

Let $P(x)$ be the polynomial, and $x_1 ... x_N$ as the set of points for which we have evaluations of $P(x)$. Call these evaluations $y_1 ... y_N$.

Think of $P(x)$ as a linear combination $\sum_i y_i L_i(x)$, where $L_i(x)$ is the polynomial that equals 1 at $x_i$ and 0 at all the other x-coordinates in the set.

Now, let us explore how these $L_i$ can be computed. Each $L_i(x)$ can be expressed as:

$$\frac{(x - x_1)(x - x_2)...(x - x_{i-1})(x - x_{i+1}) ... (x - x_N)}{(x_i - x_1)(x_i - x_2)...(x_i - x_{i-1})(x_i - x_{i+1}) ... (x_i - x_N)}$$

To see why this is true, notice that the numerator of the fraction is a polynomial that equals 0 at every evaluation point except $x_i$, and the denominator is just whatever the numerator equals at $x_i$. Hence, the result equals 0 at every evaluation point except $x_i$, and 1 at $x_i$.

Let us denote this denominator $d_i$. Let us also denote the product $(x - x_1)(x - x_2) ... (x - x_N)$ (_without_ skipping any indices) $M(x)$. We can now express $L_i(x)$ as:

$$\frac{M(x)}{d_i(x - x_i)}$$

The $d_i$ values take $O(N^2)$ time to compute (because each one takes $O(N)$ time and there are $N$ of them), but their values depend only on the evaluation points, not the specific polynomial. Hence, they can be precomputed and stored and don't need to be re-computed each time. $M(x)$ takes $O(N)$ time to compute at any specific point. Given these values, each $L_i(x)$ can be computed in $O(N)$ time. Hence, the expression:

$$P(x) = M(x) * \sum_i \frac{y_i}{d_i(x - x_i)}$$

can be computed in $O(N)$ time.

## Special case: roots of unity

In many applications of barycentric evaluation, we can choose the evaluation points. One very convenient choice is the set $\{1, \omega, \omega^2 ... \omega^{N-1}\}$, for an $\omega$ where $\omega^N = 1$ and $N$ is a power of two. The main advantage of this choice is that it removes the need to pre-compute the $d_i$ values, because there is a simple closed-form value for them.

To see why, let us try to simplify the expression for $d_i$:

$$\prod_{0 \le j \ne i < N} (\omega^i - \omega^j)$$

Notice that $\omega^{\frac{N}{2}} = -1$. We can re-express the above as:

$$2 \omega^i * \prod_{0 \le j \ne i < \frac{N}{2}} (\omega^i - \omega^j)(\omega^i + \omega^j)$$

Each term inside the product groups together $(\omega^i - \omega^j)$ and $(\omega^i - \omega^{j + \frac{N}{2}})$, as $\omega^{j + \frac{N}{2}} = -\omega^j$. The term on the left is for $(\omega^i - \omega^{i + \frac{N}{2}})$, which simplifies to $2\omega^i$ for the same reason.

Now, notice that $(\omega^i - \omega^j)(\omega^i + \omega^j)$ combines into $(\omega^{2i} - \omega^{2j})$. So we have:

$$2 \omega^i * \prod_{0 \le j \ne i < \frac{N}{2}} (\omega^{2i} - \omega^{2j})$$

Repeat this, and we have:

$$2 \omega^i * 2 \omega^{2i} * \prod_{0 \le j \ne i < \frac{N}{4}} (\omega^{4i} - \omega^{4j})$$

Keep repeating, and we get:

$$2 \omega^i * 2 \omega^{2i} * 2\omega^{4i} * ... * 2\omega^{\frac{N}{2}i}$$

The product term on the right falls away because the product eventually becomes empty. The remaining term simply becomes $N * \omega^{(N-1)i}$, or:

$$d_i = \frac{N}{\omega^i}$$

Additionally, note that $M(x) = \prod_i (x - \omega^i)$ simply becomes $x^N - 1$, so the entire calculation simplifies down to:

$$P(x) = \frac{x^N - 1}{N} * \sum_i \frac{y_i * \omega^i}{x - \omega^i}$$