拉格朗日插值学习笔记
拉格朗日插值学习笔记
应用
众所周知，在平面直角坐标系中，对于任意的 \(n\) 个点，都一定有一个不超过 \(n-1\) 次的函数与之相对应。拉格朗日插值适用于求解这 \(n\) 个点对应的函数。
思路
考虑给定的 \(n\) 个点的坐标表示为 \((x_i,y_i)\)，不难构造出如下函数：
那么此时只需要构造出符合要求的 \(g_i\) 即可。不难发现，为了使 \(f\) 符合条件，\(g_i\) 应该满足：
其中，\(x\) 表示所有 \(x_i\) 构成的集合。不难发现，当 \(n\ne x_i\) 时函数值为 \(0\)，说明该函数 \(g_i(x)\) 一定有因式 \((x-x_j)(i\ne j)\)，因此不难看出 \(g_i(x)\) 应当有因式 \(\prod_{j\ne i}(x-x_j)\)。又因为当 \(n=x_i\) 时，函数值为 \(1\)，不难构造出：
如此就构造出了符合条件的函数，此时可以 \(O(n^2)\) 求解。在单次求解的时候可以采用。### \(O(n^2)\) 写法
代码
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=2e3+5,M=998244353;
int n,k;
int x[N],y[N];
ll QuickPow(ll a,int b){
ll res=1;
while(b>0){
if(b&1)res=res*a%M;
a=a*a%M;
b>>=1;
}
return res;
}
ll Lagrange(int x[],int y[],int X){
ll res=0;
for(int i=1;i<=n;i++){
ll l_i=1;
for(int j=1;j<=n;j++){
if(i!=j)l_i=l_i*(X-x[j])%M*QuickPow(x[i]-x[j],M-2)%M;
}
res=(res+l_i*y[i]%M)%M;
}
return res;
}
int main(){
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++)scanf("%d%d",&x[i],&y[i]);
printf("%lld",(Lagrange(x,y,k)+M)%M);
return 0;
}
那么可能有人就要说了，我要是需要多次求解，那我不就似了吗，没错，你似了。所以不妨进行进一步转换，令 \(t_i=\frac{y_i}{\prod_{j\ne i}(x_i-x_j)}\)，不难看出，\(f(x)=\sum_{i=1}^{n} t_i\prod_{j\ne i}(x-x_j)\)。不妨计算出 \(\prod_{i=1}^{n}(x-x_i)\) 后将对应项除掉，接着利用多项式乘法算出每一次项前的系数后，就可以按照多项式加法的方式 \(O(n^2)\) 预处理，\(O(n)\) 求解。这样的方式适用于多组查询或需要求出函数对应项系数的情况。
代码
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=2e3+5,M=998244353;
int n,k;
ll t[N],x[N],y[N],g[N],fs[N],f[N];
ll QuickPow(ll a,int b){
ll res=1;
while(b>0){
if(b&1)res=res*a%M;
a=a*a%M;
b>>=1;
}
return res;
}
void Init(){
for(int i=0;i<n;i++){
t[i]=1;
for(int j=0;j<n;j++){
if(i==j)continue;
t[i]=t[i]*(x[i]-x[j])%M;
}
t[i]=QuickPow(t[i],M-2)*y[i]%M;
}
fs[0]=1;
for(int i=0;i<n;i++){
for(int j=n-1;j>0;j--)fs[j]=(fs[j-1]+fs[j]*(M-x[i])%M)%M;
fs[0]=fs[0]*(M-x[i])%M;
}
for(int i=0;i<n;i++){
ll inv=QuickPow(M-x[i],M-2);
g[0]=fs[0]*inv%M;
for(int j=1;j<n;j++)g[j]=(fs[j]-g[j-1])*inv%M;
for(int j=0;j<n;j++)f[j]=(f[j]+t[i]*g[j]%M)%M;
}
return ;
}
ll Lagrange(ll X){
ll res=0,Pow=1;
for(int i=0;i<n;i++){
res=(res+Pow*f[i]%M)%M;
Pow=Pow*X%M;
}
return res;
}
int main(){
scanf("%d%d",&n,&k);
for(int i=0;i<n;i++)scanf("%lld%lld",&x[i],&y[i]);
Init();
printf("%lld",(Lagrange(k)+M)%M);
return 0;
}
但也许题目中给出的 \(n\) 不接受 \(n^2\) 的预处理，那么你可以在题目不强制给定点坐标的题里大展身手。更具体的，我们可以利用一些特殊值，即尝试代入连续的点值，那么此时 \(x_i\) 是连续的，分母就变成了类似 \(\prod_{i=1}^{j-1}i\prod_{i=1}^{n+1-j}(-i)\) 的形式，这种情况只需要预处理阶乘就可以 \(O(1)\) 求解，分子也可以通过 \(O(n)\) 预处理前缀和后缀积达到 \(O(1)\) 求解的目的，那么单次求解的复杂度是 \(O(n)\)。但是这种求解方式常数较大，因为前缀后缀积也需要在每次求解时重新处理，并且要求能够快速求出点值对应的函数值，所以推荐在复杂度允许的情况下选择第二种做法。下方代码已求解函数 \(f(x)=\sum_{i=1}^{x}i^k\) 为例。
代码
void Init(){
f[0]=1;
pw[0]=0;
for(int i=1;i<=k+2;i++){
f[i]=f[i-1]*i%P;
pw[i]=(pw[i-1]+QuickPow(i,k,P))%P;
}
return ;
}
ll Lagrange(ll x,const int M){
pre[0]=nxt[k+3]=1;
for(int i=1;i<=k+2;i++)pre[i]=pre[i-1]*(x-i)%M;
for(int i=k+2;i>=1;i--)nxt[i]=nxt[i+1]*(x-i)%M;
ll res=0,y=0;
for(int i=1;i<=k+2;i++)res=(res+(((k-i)&1)?-1:1)*pw[i]*pre[i-1]%M*nxt[i+1]%M*QuickPow(f[i-1]*f[k+2-i],M-2,M))%M;
return res;
}
