拉格朗日插值
拉格朗日插值法
在数值分析中，拉格朗日插值法是以法国 \(18\) 世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日（插值）多项式。
设 \(n-1\) 次多项式 \(f(x)\) 的函数图像过 \(n\) 个已知点 \((x_1,y_1),(x_2,y_2),(x_3,y_3),\cdots,(x_n,y_n)\)，且对于 \(\forall i\neq j,x_i\neq x_j\)，那么就可以确定 \(f(x)\) 为一个 \(n-1\) 次多项式是满足条件的。
因为可以设 \(f(x)=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\cdots+a_1x+a_0\)。
那么就可以根据这 \(n\) 个点列出来 \(n\) 个不等价的方程，解出来即可确定唯一 \(f(x)\)。
但是使用高斯消元法解方程的时间复杂度为 \(\mathcal O\left(n^3\right)\) 的，有没有更高效的方法呢？
那就需要使用拉格朗日（Lagrange）插值法。
拉格朗日插值法给出了 \(f(x)\) 的表达式：
正确性
考虑到拉格朗日插值法正确当且仅当 \(\forall x_i,f(x_i)=y_i\)。
将 \(x_i\) 代入：
若 \(j\neq i\)，则在 \(\displaystyle y_j\prod\limits_{k\neq j}\dfrac{x_i-x_k}{x_j-x_k}\) 中必然存在 \(k=i\)，使得 \(x_i-x_k=0\)，从而使整个乘积为 \(0\)。
因此：
故，拉格朗日插值法是正确的。
横坐标为连续整数
这种特殊情况可以做到 \(\mathcal O(n)\)。
设 \(n-1\) 多项式 \(f(x)\)，且已知 \(f(1)=y_1,f(2)=y_2,f(3)=y_3,\cdots,f(n)=y_n\)。
则有：
对于 \(\displaystyle\prod\limits_{j=1}^n(x-j)\) 显然是可以 \(\mathcal O(n)\) 处理出来的，除以 \(x-i\) 即可得到分子。
而对于分母 \(\displaystyle\prod\limits_{j\neq i}(i-j)\)，有：
证明
可以先写一下例子：
$$ \begin{aligned} \prod_{j\neq1}(1-j)&=(1-2)(1-3)(1-4)\cdots(1-n)\\ &=(-1)^{n-1}(2-1)(3-1)(4-1)\cdots(n-1)\\ &=(-1)^{n-1}\times1\times2\times3\times\cdots(n-1)\\ &=(-1)^{n-1}(n-1)!\\ \prod_{j\neq2}(2-j)&=(2-1)(2-3)(2-4)\cdots(2-n)\\ &=(-1)^{n-2}(2-1)(3-2)(4-2)\cdots(n-2)\\ &=(-1)^{n-2}\times1\times1\times2\times\cdots(n-2)\\ &=(-1)^{n-2}(n-2)!\\ \end{aligned} $$
不难发现，$i-j$ 可以转化为 $j-i$ 计算再给整体乘 $-1$。则 $i-j<0$ 有 $n-i$ 个 $j$，即整体乘 $(-1)^{n-i}$。
因此就可以拆成阶乘，$i<j$ 部分是 $(n-i)!$，而 $i>j$ 部分是 $(i-1)!$，于是便得到上式。
只需要预处理 \(1\sim n\) 的阶乘及其逆元、\(x-i\) 的逆元与 \(\displaystyle\prod\limits_{j=1}^n(x-j)\) 即可做到单次计算复杂度 \(\mathcal O(n)\)。如果预处理 \(x-j\) 前后缀积， 则可以避免求 \(x-i\) 的逆元。
这常用于求解自然数幂次和，即 \(\displaystyle\sum\limits_{i=1}^ni^k\)。当 \(n\) 较小时可以使用线性筛做到 \(\mathcal O\left(\dfrac{n}{\ln n}\log k\right)\) 求解。（若 \(n,k\) 同阶，则复杂度为 \(\mathcal O(n)\)）
\(n\) 较大时，则可以使用拉格朗日插值求解。可以证明，\(\displaystyle\sum\limits_{i=1}^ni^k\) 一定是一个关于 \(n\) 的 \(k+1\) 次多项式，因此线性筛找 \(k+2\) 个值，插值即可。时间复杂度 \(\mathcal O(k)\)。
参考代码
这是经典题目CF622F The Sum of the k-th Powers的参考代码。使用维护前后缀积的写法。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=1e9,K=1e6,P=1e9+7;
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int y[K+2+1],fact[K+2+1],invFact[K+2+1];
void pre(int k){
fact[0]=1;
for(int i=1;i<=k+2;i++){
fact[i]=1ll*fact[i-1]*i%P;
}
invFact[k+2]=qpow(fact[k+2],P-2);
for(int i=k+1;i>=0;i--){
invFact[i]=invFact[i+1]*(i+1ll)%P;
}
static int vis[K+2+1],prime[K+2+1],size;
y[1]=1;
for(int i=2;i<=k+2;i++){
if(!vis[i]){
vis[i]=i;
prime[++size]=i;
y[i]=qpow(i,k);
}
for(int j=1;j<=size&&i*prime[j]<=k+2;j++){
vis[i*prime[j]]=prime[j];
y[i*prime[j]]=1ll*y[i]*y[prime[j]]%P;
if(i%prime[j]==0){
break;
}
}
}
for(int i=1;i<=k+2;i++){
y[i]=(y[i-1]+y[i])%P;
}
}
int f(int n,int k){
if(n<=k+2){
return y[n];
}
int ans=0;
static int pre[K+2+1],suf[K+2+1+1];
pre[0]=suf[k+3]=1;
for(int i=1;i<=k+2;i++){
pre[i]=1ll*pre[i-1]*(n-i)%P;
}
for(int i=k+2;1<=i;i--){
suf[i]=1ll*suf[i+1]*(n-i)%P;
}
for(int i=1;i<=k+2;i++){
int w=1;
if(k+2-i&1){
w=-1;
}
ans=(ans+w*1ll*y[i]*pre[i-1]%P*suf[i+1]%P*invFact[i-1]%P*invFact[k+2-i]%P)%P;
}
if(ans<0){
ans+=P;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,k;
cin>>n>>k;
pre(k);
cout<<f(n,k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
例题 AC 代码
普通的拉格朗日插值即可。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=2e3,P=998244353;
struct node{
int x,y;
}a[N+1];
int n,k;
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int f(int k){
int ans=0;
for(int i=1;i<=n;i++){
int pl=a[i].y;
for(int j=1;j<=n;j++){
if(i==j){
continue;
}
pl=1ll*pl*(k-a[j].x)%P*qpow(a[i].x-a[j].x,P-2)%P;
}
ans=(1ll*pl+ans)%P;
}
if(ans<0){
ans+=P;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n>>k;
for(int i=1;i<=n;i++){
cin>>a[i].x>>a[i].y;
}
cout<<f(k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
浙公网安备 33010602011771号