拉格朗日插值学习笔记

逛知乎时偶然看到了一个很经典的找规律填数问题，然后下面的回答基本都是 114514恶臭，突然想知道大伙是如何构造出这种能填入恶臭数字的函数的，于是就去了解了一波插值，于是就学了一波拉格朗日插值，于是就有了这篇博客。

引入

众所周知，\(n+1\)​​ 个点 \((x_i,y_i)\)​​ 可以唯一地确定一个 \(n\)​​ 次多项式，给你 \(n+1\)​​ 个点，让你求出所对应的那个多项式，怎么做？

常规的套路是待定系数法列出方程组，然后用高斯消元 \(O(n^3)\)​​​​​ 求解，显然时间复杂度过大，还容易出精度问题，还有什么更优的方法吗？有，拉格朗日插值，时间复杂度 \(O(n^2)\)​​​​​，且代码实现起来基本不会有精度问题。

说到拉格朗日插值，就不得不先提一下到底什么是插值了。简单来说，插值就是给你几个离散的点 \(\big(x,f(x)\big)\)​​​​​​​​​​​​​，让你根据这些点求出一个函数，满足这些点都在这个函数上，称为插值函数；然后一般会给你几个新的 \(x\)​​​​​​，让你求出这些 \(x\)​​​​ 所对应的 \(f(x)\)​​​​​​​​​​​​​ 是多少，称为内插，一般多用于数值分析领域。

显然，插值的方法有很多种，根据不同的插值法，能求到不同的插值函数，一般把插值法按照它们所求得的插值函数类型分类。比如你用某个插值法求出的插值函数是一个多项式，那么就称该插值法为多项式插值法，拉格朗日插值就是一个经典的多项式插值法，它可以在给定 \(n+1\) 个点的条件下求出一个 \(n\) 次多项式函数，满足这些点都在这个函数上。

举个栗子，给出 \((1,1),~(2,3),~(3,5),~(4,7)\)，问你 \(f(5)\) 是多少。

答案是 \(114514\)​​。

因为当这个函数是 \(\displaystyle f(x)=\frac {18111}2x^4-90555~x^3+\frac{633885}2x^2-452773~x+217331\) 时，有：

\[f(1)=1\\f(2)=3\\f(3)=5\\f(4)=7\\f(5)=114514 \]

这插值怎么这么臭啊（恼）

拉格朗日插值

目的是构造函数 \(F(x)\)​​​​​​，使得对于给定的点 \(x_i\)​​​​​​，都有 \(F(x_i)=y_i\)​​​​​​。怎么构造呢？我们可以让 \(F(x)=y_1f_1(x)+y_2f_2(x)+\dots+y_{n+1}f_{n+1}(x)\)​​，​且 \(f_i(x)\) 只有当 \(x=x_i\) 时取值为 1，否则取值为 0，显然这样构造出来的 \(F(x)\) 是满足条件的。

那如何构造 \(f_i(x)\)​​ 呢？简单，只需让 \(f_i(x)=\prod\limits_{j\neq i}\dfrac{x-x_j}{x_i-x_j}\)​​​ 即可。就拿上面的很臭的例子来说吧，我们有：

\[f_1(x)=\frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}\\ f_2(x)=\frac{(x-1)(x-3)(x-4)}{(2-1)(2-3)(2-4)}\\ f_3(x)=\frac{(x-1)(x-1)(x-4)}{(3-1)(3-2)(3-4)}\\ f_4(x)=\frac{(x-1)(x-2)(x-3)}{(4-1)(4-2)(4-3)}\\ \]

显然当 \(x\neq x_i\)​​ 时，\(f_i(x)\)​​ 的分子肯定有一项为 0，值为 0；当 \(x=x_i\) 时，\(f_i(x)\) 的分子分母相等，值为 1，即 \(f_i(x)\) 刚好满足上文所说的条件。

此时 \(F(x)=f_1(x)+3f_2(x)+5f_3(x)+7f_4(x)\)，试着把几个点值都带进去，发现都对的上114514无视即可，于是满足条件的插值函数 \(F(x)\)​​​​ 就成功构造出来了~~

应用&&例题

在算法竞赛中，往往会遇到这样一类问题，它们会直接或间接地给你几个点，然后让你求另一个点的值，这类问题往往会跟多项式问题挂上钩。现在我们知道了拉格朗日插值的公式，就可以利用它解决一些经典的插值问题了，这里上几个例题帮助说明。

【模板】拉格朗日插值 - 洛谷

题目大意

给出 \(n\)​​​ 个点和一个数 \(k\)​​​，让你求出根据这 \(n\)​​​ 个点确定的多项式 \(y=F(x)\)​​​ 在 \(k\)​​​ 上的值，即 \(F(k)\)​​​，答案对 998244353 取模。

\((1\leq2e3,~1\leq x_i,y_i,k<998244353),~x_i~\)两两不同。

题目分析

模板题，根据拉格朗日插值，我们有：

\[F(x)=\sum_{i=1}^n\left(y_i\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}\right) \]

直接把 \(x=k\)​​​​ 带进去暴力算就行，时间复杂度 \(O(n^2)\)​​​​，记得分母先做连乘最后再求逆元，不然复杂度多个 \(\log\)​​​，只能得 60 分。我就是每次都对分母求逆元然后过不了

代码实现

#include<bits/stdc++.h>
using namespace std;
const long long mod=998244353;
const int maxn=2e3+10;
long long x[maxn],y[maxn];
int n,k;
void exgcd(long long a,long long b,long long &x,long long &y){
	if(b==0){
		x=1,y=0;
		return ;
	}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
long long inv(long long a){
	long long x,y;
	exgcd(a,mod,x,y);
	return x;
}
int main(){
	cin >> n >> k;
	for(int i=1;i<=n;i++){
		cin >> x[i] >> y[i];
	}
	long long ans=0;
	for(int i=1;i<=n;i++){
		long long fz=y[i];
		long long fm=1;
		for(int j=1;j<=n;j++){
			if(i==j) continue;
			long long t=k-x[j];
			fm*=x[i]-x[j];
			fm%=mod;
			fz=fz*t%mod;
		}
		fm<0?fm+=mod:0;
		fm=inv(fm);
		ans+=fz*fm%mod;
		ans%=mod;
	}
	ans<0?ans+=mod:0;
	cout << ans;
} 

---

CF622F The Sum of the k-th Powers - 洛谷

题目大意

给定 \(n\) 和 \(k\)，让你求 \(\sum_{i=1}^ni^k\mod {(1e9+7)}\)。

\((1\leq n\leq 1e9,~0\leq k\leq1e6)\)​​

题目分析

乍一看好像跟拉格朗日插值没啥关系，但众所周知，次方数的前缀和是有公式的，而 \(k\)​​​​ 次方数的前缀和公式一定是个 \(k+1\)​​​​​ 次多项式。简单的证明的话请看这里，当然方法不止一种，有兴趣的同学可以自行去了解一下。

假设 \(k\)​​​​​​​​​ 次方数的前缀和公式 \(F_k(n)=\sum_{i=1}^ni^k\)​​​​​​​​​，则我们可以计算出 \(k+2\)​​​​​​​ 个点 \(\big(x,F_k(x)\big)\)​​​​​​​，然后利用拉格朗日插值求得 \(F_k(n)\)​​​​。

计算 \(k+2\)​​​​​​​​ 个点可以直接用快速幂 \(O(n\log n)\)​​​​​​​​ 求，也可以线性筛 \(x^k\)​​​​​​​​ 然后做前缀和 \(O(n)\)​​​​​ 求，问题不大，主要问题是拉格朗日插值的时间复杂度是 \(O(n^2)\)​​​​​ 的，\(\bf n\)​​​​​ 在这里是点的数量，即 \(k+2\)​​​​​，直接套板子必 T。

注意到横坐标 \(x\)​ 的值是连续的，所以可以简化运算至 \(O(n)\) 时间复杂度，推导如下：

\[F(n)=\sum_{i=1}^{k+2}\Bigg({\large y_i*f_i(n)}\Bigg)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-x_j}{x_i-x_j}\right)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\right) \]

我们可以 \(O(1)\)​ 求 \(f_i(n)=\displaystyle\prod\limits_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\)​​ ：

\[分子=\prod_{j=1,j\neq i}^{k+2}(n-j)=\prod_{j=1,j\neq i}^{i-1}(n-j)*\prod_{j=i+1,j\neq i}^{k+2}(n-j)\\ 分母=\prod_{j=1,j\neq i}^{k+2}(i-j)=\prod_{j=1,j\neq i}^{i-1}(i-j)*\prod_{j=i+1,j\neq i}^{k+2}(i-j)=(-1)^{k+2-i}*(i-1)!*(k+2-i)! \]

显然可以 \(O(n)\)​​​ 预处理分子和分母，然后 \(O(1)\)​​​ 求解 \(y_i*f_i(n)\)​​​，从而将总时间复杂度降为了 \(O(n)\)​​。

代码实现

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int maxn=1e6+6;
long long n,k;
long long prime[maxn];
bool isp[maxn];
long long y[maxn],f[maxn],invf[maxn],pre[maxn],suf[maxn];
long long qsm(long long a,long long b){
	long long ret=1;
	while(b){
		if(b&1){
			ret*=a,ret%=mod;
		}
		b>>=1,a*=a,a%=mod;
	}
	return ret;
}
void inip(){
	f[0]=invf[0]=f[1]=invf[1]=1;
	y[1]=1;
	for(int i=2;i<maxn;i++){
		f[i]=f[i-1]*i%mod;
		if(!isp[i]){
			y[i]=qsm(i,k);
			prime[++prime[0]]=i;
		}
		for(int j=1;j<=prime[0]&&i*prime[j]<maxn;j++){
			isp[i*prime[j]]=true;
			y[i*prime[j]]=y[i]*y[prime[j]]%mod;
			if(i%prime[j]==0) break;
		}
	}
}
int main(){
	cin >> n >> k;
	inip();
	pre[0]=1,suf[k+3]=1,invf[k+3]=qsm(f[k+3],mod-2);
	for(int i=1;i<=k+2;i++){
		pre[i]=pre[i-1]*(n-i)%mod;
		y[i]+=y[i-1];
		y[i]%=mod;
	}
	for(int i=k+2;i>=1;i--){
		invf[i]=invf[i+1]*(i+1)%mod;
		suf[i]=suf[i+1]*(n-i)%mod;
	}
	if(n<=k+2){
		cout << y[n];
	}
	else {
		long long ans=0;
		for(int i=1;i<=k+2;i++){
			long long fz=pre[i-1]*suf[i+1]%mod;
			long long fm=invf[i-1]*((k+2-i&1)?mod-invf[k+2-i]:invf[k+2-i])%mod;
			ans+=y[i]*fz%mod*fm%mod;
			ans%=mod;
		}
		ans<0?ans+=mod:ans;
        cout << ans;
	}
} 

参考资料

• 插值（离散数学名词）_百度百科 (baidu.com)
• 拉格朗日插值法(图文详解) - Angel_Kitty - 博客园 (cnblogs.com)
• 从拉插到快速插值求值 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
• 题解 P4781 拉格朗日插值 - attack 的博客 - 洛谷博客 (luogu.com.cn)