拉格朗日插值学习笔记

拉格朗日插值学习笔记

插值

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

插值的一般形式如下：

已知 \(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)\) 全部加起来才能符合要求。

于是设

\[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}~. \]

所以最终

\[\boldsymbol{f(x)=\sum_{i=1}^n\left(y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\right)~.} \]

这个公式特别重要，你可以什么都不会，但只需要记住这个公式就可以做题了。要求 \(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;
}

在 \(\boldsymbol{x_i}\) 连续时的做法

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

题意：求

\[\sum_{i=1}^ni^k\pmod{10^9+7}~. \]

其中 \(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\) 得

\[f(n)=\sum_{i=1}^{k+2}\left(f(i)\prod_{j\ne i}\frac{n-n_j}{n_i-n_j}\right)~, \]

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

\[\begin{aligned} f(n)&=\sum_{i=1}^{k+2}\left(f(i)\prod_{j\ne i}\frac{n-j}{i-j}\right)\\ &=\sum_{i=1}^{k+2}\left[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)\right]\\ &=\sum_{i=1}^{k+2}f(i)\frac{n^{\underline i}}{i^{\underline i}}\frac{(n-i-1)^{\underline{k+2-i}}}{(-1)^{\underline{k+2-i}}}\\ &=\sum_{i=1}^mf(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}^mf(i)\frac{(-1)^{m-i}(n-i-1)!}{i!(n-i)!(m-i)!}\\ &=\boldsymbol{\frac{n!}{(n-m-1)!}\sum_{i=1}^mf(i)\frac{(-1)^{m-i}}{i!(n-i)(m-i)!}~.} \end{aligned} \]

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

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

\[\boldsymbol{f(n)=\sum_{i=1}^mf(i)(-1)^{m-i}\frac{\textbf{pre}_{i-1}\times\textbf{suf}_{i+1}}{(i-1)!(m-i)!}}~. \]

其中

\[\begin{aligned} \text{pre}_i&=\prod_{j=1}^i(n-j)\\ \text{suf}_i&=\prod_{j=i}^n(n-j) \end{aligned} \]

这是利用了 \(\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)\)。发现问题在于：加点太快了，查询太慢了。考虑平均这两种操作。

将原式变形

\[\begin{aligned} 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}^ny_i\left(\prod_{j\ne i}(x-x_j)\prod_{k\ne i}\frac1{x_i-x_k}\right)\\ &=\sum_{i=1}^n\left(\prod_{j\ne i}(x-x_j)\prod_{k\ne i}\frac{y_i}{x_i-x_k}\right)\\ &=\prod_{j=1}^n(x-x_j)\left(\sum_{i=1}^n\frac{y_i}{\prod_{k\ne i}(x_i-x_k)}\times\frac1{x-x_i}\right)~. \end{aligned} \]

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

令

\[\begin{matrix} \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{matrix} \]

则得到最终的式子：

\[\boldsymbol{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;
}