拉格朗日插值学习笔记

基本内容

考虑一个最高次项为\(n\)的多项式函数\(f(x)\),它可以写成:

\[f(x)=\sum_{i=0}^{n}a_ix^i \]

拉格朗日插值要解决的问题是:

告诉你\(n+1\)个\(x_i\)(\(0\leq i\leq n\)),和它们对应的\(y_i=f(x_i)\)的值。让你回答询问:对某个给定的\(k\),求出\(f(k)\)的值。

朴素的想法是把\(a_0,a_1\dots a_n\)看做\(n+1\)个未知数,给定的\(n+1\)组\(x_i\)和\(y_i\)看做\(n+1\)个方程。然后用高斯消元解方程。时间复杂度是\(O(n^3)\)。

考虑更快的做法。也就是拉格朗日插值法。它的结论是:

\[f(k)=\sum_{i=0}^{n}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j} \]

可以感性理解其正确性:当\(k\)等于某个\(x_t\)时:若\(i\neq t\),则必会出现\(j=t\)的情况,也就是说\(\prod\)中必有一项的分子为\(0\),因此\(i\neq t\)的项对求和式就没有贡献。当\(i=t\)时,\(\prod\)的每一项分子分母都相等。因此\(f(k)=y_t\),这与题目条件是相符的。

至此,我们已经求出了\(f(k)\)。时间复杂度就是按这个式子暴力模拟的复杂度:\(O(n^2)\)。

求系数

我们不满足于求某个特定\(f(k)\)的值。我们还想要知道\(a_0\dots a_n\)具体是多少。

把上面的式子转化一下:\(f(k)=\sum_{i=0}^{n}\frac{y_i}{\prod_{j\neq i}x_i-x_j}\prod_{j\neq i}(k-x_j)\)。请务必注意:这个式子里只有\(k\)是变量,其他值都是常量。也就是说,我们的目标是要把它写成\(f(k)=a_0+a_1k+a_2k^2+\dots a_nk^n\)的这种形式。

其中\(\frac{y_i}{\prod_{j\neq i}x_i-x_j}\),作为常量,可以直接求出来。

那么关键的部分就是要把\(\prod_{j\neq i}(k-x_j)\)“展开”成一个关于\(k\)的求和的形式。我们先预处理出\(\prod_{j=0}^{n}(k-x_j)\):把它写成一个以\(k\)为自变量的最高次数为\(n+1\)次的多项式。然后对于每个\(i\),我们在预处理出的这个\(n+1\)次多项式的基础上,除以\((k-x_i)\)。

乘以\((k-x_i)\)和除以\((k-x_i)\)的操作都可以\(O(n)\)实现。具体来说,考虑两个生成函数:\(H(k)\)和\(G(k)\)。其中\(H(k)=G(k)(k-c)\)(请注意:\(c\)是在此处是一个常数,也就是前面式子里的\(x_i\)。而\(k\)是变量)。考虑生成函数的\(k^i\)项,可知:\(h_i=g_{i-1}-cg_i\)。那么,乘以\((k-c)\)相当于由\(G\)推\(H\),除以\((k-c)\)相当于由\(H\)推\(G\)。根据\(h_i=g_{i-1}-cg_i\)直接递推即可。

总时间复杂度\(O(n^2)\)。

乘/除以\((k-x_i)\)的代码:

void mul(int f[],int c){//mul(k-c)
	for(int i=n;i>=0;--i){
		f[i]=(i?f[i-1]:0)-a*f[i];
	}
}
void div(int f[],int c){//div(k-c)
	int lst=0;
	for(int i=n;i>=0;--i){
		lst+=a*f[i+1];
		swap(f[i],lst);
	}
}

LOJ165: 动态加点,随时询问

题目链接

朴素的做法中,加点是\(O(1)\)的;询问时按拉格朗日插值的公式模拟一遍,复杂度\(O(n^2)\)。因为有\(O(n)\)次询问,所以总时间复杂度\(O(n^3)\)。

发现朴素做法中,加点、查询的复杂度非常不平均:加点太快了、查询太慢了。于是我们想,能不能在加点的同时维护一些东西,使得两种操作的复杂度都变为\(O(n^2\log x)\)或\(O(n^2)\)。

考虑在询问时还是枚举\(i\)。我们预先维护好求积式中分母、分子的值。

  • 分母:对于每个\(i\),维护一个\(\text{fm}[i]=\frac{1}{\prod_{j\neq i}(x_i-x_j)}\),这个可以在加点时顺便维护出来。
  • 分子:在询问时,先特判\(k\)等于其中某个\(x_t\)的情况(直接输出\(y_t\))。否则就求一个变量\(\text{fz}=\prod_{j=0}^{n}(k-x_j)\),然后在枚举\(i\)时除以\((k-x_i)\)即可。

于是,一次询问的答案可以写为:\(\sum_{i=0}^{n}y_i\cdot \text{fm}[i]\cdot \text{fz}\cdot\frac{1}{k-x_i}\)。

如果在线回答询问,就需要\(O(\log x)\)求逆元。时间复杂度\(O(n^2\log x)\)。如果把询问离线下来,预处理所有逆元,则时间复杂度\(O(n^2)\)。两种做法的参考代码:在线(TLE 70分)离线预处理逆元(AC 100分)

\(x\)取值连续时

在\(x\)取值为连续的\(n+1\)个自然数时,我们可以把拉格朗日插值算法的时间复杂度优化为\(O(n)\)。

回到最开始的式子,即:\(f(k)=\sum_{i=0}^{n}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j}\)。

首先,我们可以预处理\(k-x_i\)的前缀积、后缀积。具体地,预处理出:\(\text{pre}[i]=\prod_{j=0}^{i}(k-x_j)\), \(\text{suf}[i]=\prod_{j=i}^{n}(k-x_j)\)。这样,上式中的分子就等于\(\prod_{j\neq i}(k-x_j)=\text{pre}[i-1]\cdot \text{suf}[i+1]\)。

然后考虑分母,即:\(\prod_{j\neq i}(x_i-x_j)\)。因为\(x\)取值连续,所以这个式子就等于\(\prod_{j\neq i}(i-j)\)。分\(j<i\)、\(j>i\)讨论,这个式子可以化为:\(i!\cdot(-1)^{n-i}\cdot(n-i)!\)。

综上所述,此时:

\[f(k)=\sum_{i=0}^{n}y_i\frac{\text{pre}[i-1]\cdot \text{suf}[i+1]}{i!\cdot(-1)^{n-i}\cdot(n-i)!} \]

按这个式子直接模拟,时间复杂度\(O(n)\)。

自然数幂求和

给定\(n\), \(k\),求:\(S(n,k)=\sum_{i=1}^{n}i^k\)。一般\(n\)很大(\(10^9\)级别)而\(k\)较小(\(10^5\)级别)。

其中一个经典的方法是把\(S(n,k)\)中\(n\)视为定值,然后写出其指数生成函数:\(\sum_{i=0}^{\infty}S(n,i)\frac{x^i}{i!}\)。再根据\(e^x=\sum_{i=0}^{\infty}\frac{x^i}{i!}\)进行一系列代换。最后用多项式求逆求出答案。时间复杂度\(O(k\log k)\)。这种方法的好处是可以一次性求出所有\(S(n,0\dots k)\)。

本文要介绍的拉格朗日插值法,一次只能求出一个\(S(n,k)\),但是时间复杂度只有\(O(k)\),而且更好写。

这种方法基于一个结论:

\(S(n,k)\)是以\(n\)为自变量的\(k+1\)次多项式。(请注意,这里与上一种方法恰恰相反,它把\(n\)看做变量、\(k\)看做定值)。

例如,\(\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6}\),这就是一个\(3\)次多项式。

因此我们可以带入\(k+2\)个点,然后用拉格朗日插值法直接求出\(S(n,k)\)。如果带入的\(k+2\)个点为\(1\dots k+2\),则时间复杂度为\(O(k)\)。

至于预处理\(1^k,2^k,\dots ,3^k\),可以用线性筛。在质数处暴力快速幂。因为质数只有\(O(\frac{k}{\log k})\)个。所以复杂度就是\(O(\frac{k}{\log k}\cdot \log k)=O(k)\)。

例题:LOJ2026「JLOI / SHOI2016」成绩比较

上一篇:前端链接


下一篇:[题解] ABC172E NEQ