拉格朗日插值法

简介

众所周知，给定$n + 1$个横坐标不相同的点，可以唯一确定一个$n$次的多项式。那么如何求出这个多项式？最直观的做法就是列方程求解。但是这样需要$\Theta(n^3)$的时间来计算。而拉格朗日插值法则通过构造的方法，得到了一个经过$n + 1$个点的$n$次多项式。
具体的过程是这样的，假设现在我们得到了$n + 1$个点：
$$ (x_0, y_0), \;(x_1, y_1),\; \dots,\;(x_n, y_n) $$

设拉格朗日基本多项式为：
$$ \ell_j(x) = \prod_{i = 0, i \neq j}^n {x - x_i \over x_j - x_i} \tag{1.1} $$

这个基本多项式构造十分巧妙，因为注意到$\ell_j(x_j) = 1$，并且$\ell_j(x_i) = 0, \;\forall i \neq j$。那么，接着构造出这个$n$次多项式：
$$ P(x) = \sum_{i = 0}^n y_i\ell_i(x) \tag{1.2} $$

根据基本多项式的性质，我们可以知道$P(x_i) = y_i$，也就是经过了这$n + 1$个点。
通过简单的多项式乘法和多项式除法就可以在$\Theta(n^2)$的时间求出这个多项式的系数表达。
拉格朗日插值公式形式比较优美，接下来将介绍拉格朗日插值法另外一个应用。

自然数的幂的前缀和

这是一个十分经典的问题：

对于给定的$n$和$k$，要求求出：
$$ \sum_{a = 1}^n a^k $$

通常$n$会比较大，而$k$只有几千或者几万。

当$k = 1$时，我们知道公式为$n(n + 1) / 2$，是一个二次多项式。
更一般的，答案一定是$k + 1$次多项式。通过差分可以处理这种多项式，具体的可以参考我之前的博客或者有关第二类斯特林数的资料。
现在我们来考虑使用拉格朗日插值法来获得答案多项式。
首先如果我们得知了$n = 0, 1, \dots, k + 1$处的答案$f(n)$，那么给定的$n$处的答案可以写成：
$$ \begin{aligned} f(n) & = \sum_{i = 0}^{k + 1} f(i) {(n - 0)(n - 1)\cdots[n - (i - 1)][(n - (i + 1)]\cdots[n - (k + 1)] \over (i - 0)(i - 1)\cdots[i - (i - 1)][(i - (i + 1)]\cdots[i - (k + 1)]} \\ & = \sum_{i = 0}^{k + 1} f(i) {\prod_{j = 0}^{i - 1} (n - j) \prod_{j = i + 1}^{k + 1} (n - j) \over i!(-1)^{k - i + 1}(k + 1 - i)!} \\ & = \sum_{i = 0}^{k + 1} (-1)^{k - i + 1}f(i) {\prod_{j = 0}^{i - 1} (n - j) \prod_{j = i + 1}^{k + 1} (n - j) \over i!(k + 1 - i)!} \end{aligned} $$

注意到后面的分式中，分子是一个前缀积乘以一个后缀积，而分母是两个阶乘。这些都可以在$\Theta(k)$的时间内求出。
现在剩下的问题就是如何求出$f(0), f(1), \dots, f(k + 1)$了。由于$g(x) = x^k$是个完全积性函数，所以我们可以通过欧拉筛法求出$g$函数前面的一些值。具体的就是对于质数采取直接快速幂，合数则拆出任意一个因子来算，通常是欧拉筛法中可以顺便求得的最小质因子。根据素数定理，素数大约有$O(k/\ln k)$个。每次快速幂需要花费$\Theta(\log k)$的时间，因此总的时间复杂度可以估计为$\Theta(k)$，是一个非常优秀的算法。
上面的方法具有通用性，只要我们可以快速的求出某个$k$次多项式的前$k + 1$个值，那么剩下的部分可以使用拉格朗日插值法在$\Theta(k)$的时间内完成计算。