拉格朗日插值学习笔记

合集 - 学习笔记(20)

1.数论总结2024-08-192.平衡树总结2024-08-203.扫描线总结2024-08-244.组合数学总结2024-08-235.卡特兰数、Prufer 序列、BSGS 总结2024-08-236.概率和期望总结2024-08-237.基环树和笛卡尔树总结2024-11-24

8.拉格朗日插值学习笔记2024-11-27

9.虚树01-2710.同余最短路01-2111.CDQ 分治 && 整体二分01-2012.点分治&&点分树01-1513.差分约束 + 01BFS2024-11-2914.二分图2024-10-2315.网络流02-1316.0/1 分数规划02-0417.后缀数组 SA02-2318.Manacher02-2319.回文自动机 PAM02-2420.后缀自动机 SAM03-09

收起

拉格朗日插值学习笔记

插值

什么是插值？插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。

插值的一般形式如下：

已知 $n$ 个点 $P_1(x_1,y_1),P_2(x_2,y_2),\dots ,P_n(x_n,y_n)$，求 $n-1$ 次多项式 $f(x)$ 满足

$$f(x_i)=y_i,\mathrm{\forall }i\in [1,n].$$

初次学习时，可以理解为待定系数法求解析式。只不过解方程（即高斯消元）是 $O(N^3)$ 的，而插值可以在 $O(N^2)$ 的时间内计算。并且，插值可以求逆元。

拉格朗日插值

拉格朗日插值是众多插值方式中的一种。设第 $i$ 个点 $P_i(x_i,y_i)$ 在 $x$ 轴上的投影为 $P_i^{\prime }(x_i,0)$。人话：作垂直。

构造 $n$ 个函数 $f_1(x),f_2(x),\dots ,f_n(x)$，使得对于第 $i$ 个函数 $f_i(x)$，其图像过 $\{\begin{array}{l}{P}_i(x_i,y_i)\\ P_j^{\prime }(x_j,0)& j\ne i\end{array}$ ，则题目所求的函数为 $f(x)=\sum f_i(x)$。

如此构造的原因：

每个 $f_i(x)$ 都需要过 $P_i$，但不能过 $\mathrm{\forall }{P}_j(j\ne i)$。也就是在 $x_i$ 处点值为 $y_i$，而在 $\mathrm{\forall }{P}_j(j\ne i)$ 处的值都为 $0$。这样构造出来的 $f_i(x)$ 全部加起来才能符合要求。

于是设

$$f_i(x)=a\prod _{j\ne i}(x-x_j),$$

将点 $P_i(x_i,y_i)$ 代入得

$$a=\frac{y_i}{\prod _{j\ne i}(x_i-x_j)},$$

则

$$f_i(x)=y_i\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}=y_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}.$$

所以最终

$$\begin{array}{}\text{(1)}& \boxed{{\displaystyle f(x)=\sum _{i=1}^n\left(y_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\right).}}\end{array}$$

这个公式特别重要，你可以什么都不会，但只需要记住这个公式就可以做题了。要求 $f(k)$ 的时候，直接将 $k$ 代入上式中就行。

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

constexpr int MAXN=2005,MOD=998244353;
int n,k,x[MAXN],y[MAXN];
int power(int a,int b){
	int res=1;
	for(;b;a=a*a%MOD,b>>=1)
		if(b&1)
			res=res*a%MOD;
	return res;
}
void add(int&x,int y){
	x=x+y>=MOD?x+y-MOD:x+y;
}
int lagrange(int k){
	int res=0;
	for(int i=1,s1,s2;i<=n;++i){
		s1=y[i],s2=1;
		for(int j=1;j<=n;++j){
			if(i==j) continue;
			s1=s1*(k-x[j]+MOD)%MOD;
			s2=s2*(x[i]-x[j]+MOD)%MOD;
		}
		add(res,s1*power(s2,MOD-2)%MOD);
	}
	return res;
}

signed main(){
	n=read(),k=read();
	for(int i=1;i<=n;++i) x[i]=read(),y[i]=read();
	printf("%lld\n",lagrange(k));
	return 0;
}

在 ${\mathit{x}}_{\mathit{i}}$ 连续时的做法

例题：[CF622F] The Sum of the k-th Powers。

题意：求

$$\sum _{i=1}^n{i}^k(\mathrm{mod}{10}^9+7).$$

其中 $1\le n\le {10}^9$，$0\le k\le {10}^6$。

这道题的前半部分是证明原式是一个关于 $\mathit{n}$ 的 $\mathit{k}\mathbf{+}\mathbf{1}$ 次多项式。证法不是我们关心的内容，我们只关心这道题的后半部分：$x$ 取值连续时的拉格朗日插值。

插一个 $k+1$ 次多项式需要 $k+2$ 个点。我们先将 $n,k+2$ 代入 $(1)$ 式中的 $x,n$ 得

$$f(n)=\sum _{i=1}^{k+2}(f(i)\prod _{j\ne i}\frac{n-n_j}{n_i-n_j}),$$

由于 $n$ 和 $f(n)$ 的连续取值都已知，所以直接代入所有的 $n_i$ 和 $n_j$ 得

$$\begin{array}{rl}f(n)& =\sum _{i=1}^{k+2}(f(i)\prod _{j\ne i}\frac{n-j}{i-j})\\ & =\sum _{i=1}^{k+2}[f(i)\left(\prod _{j=1}^{i-1}\frac{n-j}{i-j}\right)\left(\prod _{j=i+1}^{k+2}\frac{n-j}{i-j}\right)]\\ & =\sum _{i=1}^{k+2}f(i)\frac{n^{\underset{―}{i}}}{i^{\underset{―}{i}}}\frac{(n-i-1{)}^{\underset{―}{k+2-i}}}{(-1{)}^{\underset{―}{k+2-i}}}\\ & =\sum _{i=1}^{m}f(i)\frac{n!}{i!(n-i)!}\frac{(n-i-1)!}{(n-m-1)!(-1{)}^{m-i}(m-i)!},\text{ 其中 }m=k+2\\ & =\frac{n!}{(n-m-1)!}\sum _{i=1}^{m}f(i)\frac{(-1{)}^{m-i}(n-i-1)!}{i!(n-i)!(m-i)!}\\ & =\frac{n!}{(n-m-1)!}\sum _{i=1}^{m}f(i)\frac{(-1{)}^{m-i}}{i!(n-i)(m-i)!}.\end{array}$$

所以

$$\boxed{{\displaystyle f(n)=\frac{n!}{(n-m-1)!}\sum _{i=1}^{m}f(i)\frac{(-1{)}^{m-i}}{i!(n-i)(m-i)!}.}}$$

复杂度 $O(m\log V)$，预处理逆元和 $f(i)$ 即可做到 $O(m)$。

实际上还有一种不用这么复杂的解决方案。将原式化为

$$\boxed{{\displaystyle f(n)=\sum _{i=1}^{m}f(i)(-1{)}^{m-i}\frac{{\text{pre}}_{i-1}\times {\text{suf}}_{i+1}}{(i-1)!(m-i)!}.}}$$

其中

$$\begin{array}{rl}{\text{pre}}_i& =\prod _{j=1}^i(n-j)\\ {\text{suf}}_i& =\prod _{j=i}^n(n-j)\end{array}$$

这是利用了 $\prod _{j\ne i}(n-j)=[\prod _{j<i}(n-j)][\prod _{j>i}(n-j)]$ 的原理。这种方法的复杂度也是 $O(m\log V)$，预处理逆元和 $f(i)$ 也可以做到 $O(m)$。

#include<bits/stdc++.h>
#define int long long
using namespace std;
 
constexpr int MAXN=1e6+5,MOD=1e9+7;
int n,k,pre[MAXN],suf[MAXN],fac[MAXN];
int power(int a,int b){
	int res=1;
	for(;b;a=a*a%MOD,b>>=1)if(b&1)res=res*a%MOD;
	return res;
}
int lagrange(int k){
	pre[0]=suf[k+3]=fac[0]=1;
	for(int i=1;i<=k+2;++i) pre[i]=pre[i-1]*(n-i)%MOD;
	for(int i=k+2;i;--i) suf[i]=suf[i+1]*(n-i)%MOD;
	for(int i=1;i<=k+2;++i) fac[i]=fac[i-1]*i%MOD;
	int res=0;
	for(int i=1,y=0,s1,s2;i<=k+2;++i){
		y=(y+power(i,k))%MOD;
		s1=pre[i-1]*suf[i+1]%MOD;
		s2=fac[i-1]*((k-i)&1?-1:1)*fac[k+2-i]%MOD;
		res=(res+y*s1%MOD*power(s2,MOD-2)%MOD)%MOD;
	}
	return res;
}
 
signed main(){
	scanf("%lld%lld",&n,&k);
	printf("%lld\n",(lagrange(k)+MOD)%MOD);
	return 0;
}

重心拉格朗日插值

拉格朗日插值还有一种特殊情况，那就是动态加点，随时询问。例题：LOJ #165. 拉格朗日插值。

在朴素做法中，加点是 $O(1)$ 的，查询是 $O(n^2)$ 的，再乘上 $O(n)$ 次询问，总体复杂度达到了 $O(n^3)$。发现问题在于：加点太快了，查询太慢了。考虑平均这两种操作。

将原式变形

$$\begin{array}{rl}f(x)& =\sum _{i=1}^n\left(y_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\right)\\ & =\sum _{i=1}^n{y}_i(\prod _{j\ne i}(x-x_j)\prod _{k\ne i}\frac{1}{x_i-x_k})\\ & =\sum _{i=1}^n(\prod _{j\ne i}(x-x_j)\prod _{k\ne i}\frac{y_i}{x_i-x_k})\\ & =\prod _{j=1}^n(x-x_j)(\sum _{i=1}^n\frac{y_i}{\prod _{k\ne i}(x_i-x_k)}\times \frac{1}{x-x_i}).\end{array}$$

注意最后一步中，我们为了让最外层的 $\mathrm{\Pi }$ 脱离 $i$ 的限制，在 $\mathrm{\Sigma }$ 里面又乘上了一个 ${\displaystyle \frac{1}{x-x_i}}$。

令

$$\begin{array}{c}{\displaystyle \ell (x)=\prod _{j=1}^n(x-x_j)}\\ {\displaystyle w_i=\frac{y_i}{\prod _{k\ne i}(x_i-x_k)}}\end{array}$$

则得到最终的式子：

$$\boxed{{\displaystyle f(x)=\ell (x)\sum _{i=1}^n\frac{w_i}{x-x_i}.}}$$

在新增加一个点之后，我们只需 $O(n\log V)$ 更新 $w_i$；在询问时，$\ell (x)$ 只需 $O(n)$ 计算，加上计算逆元的复杂度也是 $O(n\log V)$。

constexpr int MAXN=3005,MOD=998244353;
int n,x[MAXN],y[MAXN],w[MAXN],cnt;
int power(int a,int b){
	int res=1;
	for(;b;a=a*a%MOD,b>>=1)if(b&1)res=res*a%MOD;
	return res;
}

signed main(){
	n=read();
	while(n--){
		int opt=read(),X=read();
		if(opt==1){
			int Y=read();
			x[++cnt]=X,y[cnt]=w[cnt]=Y;
			for(int i=1;i<cnt;++i){
				w[i]=w[i]*power(x[i]-x[cnt],MOD-2)%MOD;
				w[cnt]=w[cnt]*power(x[cnt]-x[i],MOD-2)%MOD;
			}
		}else{
			int s1=1,s2=0;
			for(int i=1;i<=cnt;++i)
				if(X==x[i]){
					write(y[i]);
					goto byby;
				}
			for(int i=1;i<=cnt;++i){
				s1=s1*(X-x[i])%MOD;
				s2=(s2+w[i]*power(X-x[i],MOD-2)%MOD)%MOD;
			}
			write((s1*s2%MOD+MOD)%MOD);
			byby:;
		}
	}
	return fw,0;
}