拉格朗日插值
拉格朗日插值
插值
什么是插值？插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。
插值的一般形式如下：
已知 \(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~,\quad\forall i\in[1,n]~. \]
初次学习时，可以理解为待定系数法求解析式。只不过解方程（即高斯消元）是 \(O(N^3)\) 的，而插值可以在 \(O(N^2)\) 的时间内计算。并且，插值可以求逆元。
拉格朗日插值
拉格朗日插值是众多插值方式中的一种。设第 \(i\) 个点 \(P_i(x_i,y_i)\) 在 \(x\) 轴上的投影为 \(P_i'(x_i,0)\)。人话：作垂直。
构造 \(n\) 个函数 \(f_1(x),f_2(x),\dots,f_n(x)\)，使得对于第 \(i\) 个函数 \(f_i(x)\)，其图像过 \(\begin{cases}P_i(x_i,y_i)\\P_j'(x_j,0)&j\ne i\end{cases}\) ，则题目所求的函数为 \(f(x)=\sum f_i(x)\)。
如此构造的原因：
每个 \(f_i(x)\) 都需要过 \(P_i\)，但不能过 \(\forall P_j~(j\ne i)\)。也就是在 \(x_i\) 处点值为 \(y_i\)，而在 \(\forall P_j~(j\ne i)\) 处的值都为 \(0\)。这样构造出来的 \(f_i(x)\) 全部加起来才能符合要求。
于是设
将点 \(P_i(x_i,y_i)\) 代入得
则
所以最终
这个公式特别重要，你可以什么都不会，但只需要记住这个公式就可以做题了。要求 \(f(k)\) 的时候，直接将 \(k\) 代入上式中就行。
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;
}
求系数时的做法
用于需要多次求解函数值时的情况。最终的 \(f\) 就是各项系数。
void lagrange(int n){
for(int i=1;i<=n;i++){
t[i]=1;
for(int j=1;j<=n;j++) if(i!=j) t[i]=t[i]*(x[i]-x[j])%MOD;
t[i]=y[i]*power(t[i],MOD-2)%MOD;
}
fg[1]=1;
for(int i=1;i<=n;i++)
for(int j=n;j;j--)
fg[j]=(fg[j-1]+fg[j]*(MOD-x[i])%MOD)%MOD;
for(int i=1;i<=n;i++){
ll inv=power(MOD-x[i],MOD-2);
for(int j=1;j<=n;j++){
g[j]=(fg[j]-g[j-1])*inv%MOD;
f[j]=(f[j]+g[j]*t[i]%MOD+MOD)%MOD;
}
}
}
在 \(\boldsymbol{x_i}\) 连续时的做法
例题：[CF622F] The Sum of the k-th Powers。
题意：求
其中 \(1\le n\le10^9\)，\(0\le k\le10^6\)。
这道题的前半部分是证明原式是一个关于 \(\boldsymbol{n}\) 的 \(\boldsymbol{k+1}\) 次多项式。证法不是我们关心的内容，我们只关心这道题的后半部分：\(x\) 取值连续时的拉格朗日插值。
插一个 \(k+1\) 次多项式需要 \(k+2\) 个点。我们先将 \(n,k+2\) 代入 \((1)\) 式中的 \(x,n\) 得
由于 \(n\) 和 \(f(n)\) 的连续取值都已知，所以直接代入所有的 \(n_i\) 和 \(n_j\) 得
所以
复杂度 \(O(m\log V)\)，预处理逆元和 \(f(i)\) 即可做到 \(O(m)\)。
实际上还有一种不用这么复杂的解决方案。将原式化为
其中
这是利用了 \(\prod_{j\ne i}(n-j)=\left[\prod_{j<i}(n-j)\right]\hspace{-0.2em}\left[\prod_{j>i}(n-j)\right]\) 的原理。这种方法的复杂度也是 \(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)\)。发现问题在于：加点太快了，查询太慢了。考虑平均这两种操作。
将原式变形
注意最后一步中，我们为了让最外层的 \(\Pi\) 脱离 \(i\) 的限制，在 \(\Sigma\) 里面又乘上了一个 \(\dfrac1{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;
}
应用
拉格朗日插值要么用在直接解决某些可以转化为多项式的问题，要么用于优化可以写成多项式形式的 DP。
需要知道：将一个函数前缀和之后会升次，将一个函数差分之后会降次。一些题目经常用到。
还能解决像 CF622F 这样 \(\sum_{i=1}^ni^k\) 的问题，一些题目也会用到。
浙公网安备 33010602011771号