拉格朗日插值法

拉格朗日插值法

简介

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

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

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

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

自然数的幂的前缀和

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

对于给定的nnkk,要求求出:
a=1nak \sum_{a = 1}^n a^k

通常nn会比较大,而kk只有几千或者几万。

k=1k = 1时,我们知道公式为n(n+1)/2n(n + 1) / 2,是一个二次多项式。
更一般的,答案一定是k+1k + 1次多项式。通过差分可以处理这种多项式,具体的可以参考我之前的博客或者有关第二类斯特林数的资料。
现在我们来考虑使用拉格朗日插值法来获得答案多项式。
首先如果我们得知了n=0,1,,k+1n = 0, 1, \dots, k + 1处的答案f(n)f(n),那么给定的nn处的答案可以写成:
f(n)=i=0k+1f(i)(n0)(n1)[n(i1)][(n(i+1)][n(k+1)](i0)(i1)[i(i1)][(i(i+1)][i(k+1)]=i=0k+1f(i)j=0i1(nj)j=i+1k+1(nj)i!(1)ki+1(k+1i)!=i=0k+1(1)ki+1f(i)j=0i1(nj)j=i+1k+1(nj)i!(k+1i)! \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}

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


赞 • 3 人点赞 1 条评论 Issue 页面