重心插值（Barycentric interpolation ）是拉格朗日插值的变换。

有时候需要通过一组多项式点值直接计算另一个不同点处的值，例如，$p(x)$ 是一个度为 2 的多项式，可以在 O(N) 时间内用 $p(0), p(1), p(2)$ 直接计算 $p(100)$ 处的值，而无须重构多项式。

回顾下拉格朗日插值 $$ P(x) = \sum^{n}_{i=1}y_iL_i(x) $$ 其中 $L_i(x)$ 是拉格朗日基 $$ li(x)=\prod^{n-1}{j=1,j \ne i}\frac{x-x_j}{x_i-x_j}=\frac{x-x_0}{x_i-x0}...\frac{x-x{i-1}}{xi-x{i-1}}\frac{x-x_{i+1}}{xi-x{i+1}}...\frac{x-x_{n-1}}{xi-x{n-1}} $$ 分子在除了 $x_i$ 外的每个评估点都等于 0，分母在 $x_i$ 处等于分子。所以 $L_i(x)$ 在 $x_i$ 处等于 1，非 $x_i$ 处等于 0

令 $M(x)=\prod_{j=0}^{n-1}(x-x_j)=(x-x_0)(x-x1)...(x-x{n-1})$，$l_i(x)$ 的分母用 $d_i$（重心权重的倒数） 表示，得到 $$ L_i(x)=\frac{M(x)}{d_i(x-x_i)} $$ 得到 $$ P(x) = M(x) * \sum_i\frac{y_i}{d_i(x-x_i)} $$ 其中，$d_i$ 可以预计算，$M(x)$ 的时间复杂度为 O(N)，所以 $L_i(x)$ 的时间复杂度为 O(N)，所以 $P(x)$ 的时间复杂度为 O(N)

我们可以选择求值点以避免预计算 $d_i$，比如 $\left{1, w, w^2, ...,w^{n-1}\right}$，$n$ 是 2 的幂，即 $n=2^k$。有 $$ di=\prod{0 \le j \ne i < n}(w^i-w^j) =(w^i-w^{i+\frac{n}{2}})\prod_{0 \le j \ne i < \frac{n}{2}}(w^i-w^j)(w^i-w^{j+\frac{n}{2}}) $$ 因为 $w^{\frac{n}{2}}=-1$，所以有 $$ di=2w^i\prod{0 \le j \ne i < \frac{n}{2}}(w^i-w^j)(w^i+w^{j})=2w^i\prod_{0 \le j \ne i < \frac{n}{2}}(w^{2i}-w^{2j}) $$ 重复以上步骤，得到 $\text{log}_2n$ 项乘积 $$ d_i=2w^i 2w^{2i} 2w^{4i} ... 2w^{n/2i}=n * w^{(n-1)i} $$ 因为 $w^n=1$，所以有 $$ d_i=\frac{n}{w^i} $$ 另外，对于 2-adic 子群，有 $$ zH(x) = \prod{i=0}^{n-1}(x-w^i) =x^n-1 $$ 证明：因为 $w^{\frac{n}{2}}=-1$，当 $n=4$ 时，有 $$ \begin{align} & (x-w^0)(x-w^1)(x-w^2)(x-w^3) \ &=(x-1)(x-w)(x+1)(x-w^3)\ &=(x^2-1)(x-w)(x+w)\ &=(x^2-1)(x^2-w^2)\ &=(x^2-1)(x^2+1)\ &=x^4-1 \end{align} $$ 所以有 $$ M(x)=\prod_{j=0}^n(x-x_j)=(x-x_0)(x-x1)...(x-x{n})=x^n-1\ $$

最终得到 $$ P(x)=\frac{x^n-1}{n} \sum_i\frac{y_i w^i}{x-w^i} $$

参考

Barycentric interpolation

A quick barycentric evaluation tutorial