学习笔记:拉格朗日插值

 


建议看: cmd:从拉插到快速插值求值

拉格朗日插值

首先我们知道一个事是: n+1 个不同点可以唯一确定一个 n 次多项式

那么当你知道了这 n+1 个点之后就可以通过拉插插出这个多项式是啥。

定义:

对于 n 个点 (xi,yi) ,我们希望找到 fi 满足:

fi(j)={yij=xi0j=xk(ki)anythingotherwise

fi 也是 n1 次的,那么我们可以得到多项式 F=ifi ,这样 F 满足他是一个经过这 n 个点的 n1 次多项式。

那么对于每个 fi ,我们尝试用因式构造,简单来说,为了满足 fi(xj)=0,ji ,我们使他含有 (xxj),ji 的因式,此时再去满足 fi(xi)=yi ,我们先把 xi 带进去,就会得到 j,ji(xixj) ,这时我们给整个式子乘上 yij,ji(xixj) ,就满足了所有条件。

所以我们获得了拉插公式:

F=iyij,jixxjxixj

例题:

P4781 【模板】拉格朗日插值

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=998244353,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){
char c=getchar();ll x=0;bool f=0;
while(!isdigit(c))f=c=='-'?1:0,c=getchar();
while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
return f?-x:x;
}
inline ll qpow(ll x,ll y){
ll ans=1;
while(y){
if(y&1)ans=ans*x%mod;
x=x*x%mod;y>>=1;
}
return ans;
}
int n,K,x[N],y[N];
signed main(){
#ifndef ONLINE_JUDGE
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
n=read(),K=read();
for(int i=1;i<=n;++i)
x[i]=read(),y[i]=read();
ll ans=0;
for(int i=1;i<=n;++i){
ll i1=1,i2=1;
for(int j=1;j<=n;++j){
if(j==i)continue;
i1=i1*(K-x[j])%mod;i2=i2*(x[i]-x[j])%mod;
}
i1=(i1+mod)%mod;i2=qpow((i2+mod)%mod,mod-2);
ans+=i1*i2%mod*y[i]%mod;
}
cout<<(ans%mod+mod)%mod;
}

xi 取值连续时的插值

我们观察式子

F=iyij,jixxjxixj

对于 xxj ,我们做一个前缀积和后缀积就可以 O(n) 解决,对于 xixj 因为他们取值连续所以就变成了 k1!k2! 的形式,同样可以 O(n) 做。

例题:

CF622F The Sum of the k-th Powers

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int ll
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){
char c=getchar();ll x=0;bool f=0;
while(!isdigit(c))f=c=='-'?1:0,c=getchar();
while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
return f?-x:x;
}
inline ll qpow(ll x,ll y){
ll ans=1;
while(y){
if(y&1)ans=ans*x%mod;
x=x*x%mod;y>>=1;
}
return ans;
}
int y[N],n,K,pre[N],hou[N],jc[N];
signed main(){
#ifndef ONLINE_JUDGE
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
n=read(),K=read();
y[0]=0;jc[0]=1;
for(int i=1;i<=K+2;++i){
y[i]=(y[i-1]+qpow(i,K))%mod;
jc[i]=jc[i-1]*i%mod;
}
if(n<=K+2)return cout<<y[n],0;
pre[0]=n;
for(int i=1;i<=K+2;++i){
pre[i]=pre[i-1]*(n-i)%mod;
}
hou[K+3]=1;
for(int i=K+2;i>=0;--i)
hou[i]=hou[i+1]*(n-i)%mod;
int ans=0;
for(int i=0;i<=K+2;++i){
ans+=(i==0?1ll:pre[i-1])*hou[i+1]%mod*qpow(jc[i],mod-2)%mod*qpow(jc[K+2-i],mod-2)%mod*y[i]%mod*((K+2-i)&1?-1:1);
}
cout<<(ans%mod+mod)%mod;
}

重心拉格朗日插值法

感觉没什么意义,就是瞪了两眼式子省去了不必要的计算。

还是看式子:

F=iyij,jixxjxixj

我们每加入一个点对于每个 ixxj 变化可以 O(1) 算, xixj 也可以 O(1) 算。

插多项式系数

直接插一个值出来是 O(n2) 的,但是拉格朗日插值同样可以在 O(n2) 的时间复杂度内插出这个多项式的系数,而且支持取模。

还是得改写式子:

首先对于整个式子

F=iyij,jixxjxixj

yiji1xixj 是常数,姑且写为 ci ,可以 O(n2) 求出。对于 jixxj ,写为 jxxjxxi ,我们把 jxxj 写为 m

现在式子长成了

F=icimxxi

首先 ci 是常数,不会对 x 次数造成影响,主要造成影响的是 m ,它包含 n 个形如 xxj 的式子,最后又除了一个 xxi ,学过二项式定理的应该知道 x 的次数就是由每个 xxj 的式子提供来的,考虑 dp ,设 dpi 表示 xi 的系数,那么每次考虑一个 xxj 的时候,要么选择 x ,有 dpidpi1 ,如果选择 xxi 的话,就有 dpi=xidpi ,所以最终有 dpi=dpi1xidpi

如果只考虑 m 的话,直接 dp 是 O(n2) 的,但是有了 1xxi 之后朴素实现就是 O(n3) 的了,考虑有 1xxi 之后相当于一个退背包的过程,所以可以对于每个 i 直接 O(n) 做,所以也是 O(n2) 的。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int ll
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){
char c=getchar();ll x=0;bool f=0;
while(!isdigit(c))f=c=='-'?1:0,c=getchar();
while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
return f?-x:x;
}
inline ll qpow(ll x,ll y){
ll ans=1;
while(y){
if(y&1)ans=ans*x%mod;
x=x*x%mod;y>>=1;
}
return ans;
}
int n,x[N],y[N],c[N],f[N],g[N],ff[N];
signed main(){
#ifndef ONLINE_JUDGE
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
n=read();
for(int i=1;i<=n;++i)
x[i]=read(),y[i]=read();
for(int i=1;i<=n;++i){
c[i]=1;
for(int j=1;j<=n;++j)
if(j!=i)c[i]=c[i]*(x[i]-x[j])%mod;
c[i]=(qpow(c[i],mod-2)*y[i]%mod+mod)%mod;
}
f[0]=1;
for(int i=1;i<=n;++i){
for(int j=n-1;j;--j)
f[j]=((f[j-1]+f[j]*(-x[i]))%mod+mod)%mod;
f[0]=(f[0]*(-x[i])%mod+mod)%mod;
}
for(int i=1;i<=n;++i){
ll ny=qpow(((-x[i])%mod+mod)%mod,mod-2);
g[0]=f[0]*ny%mod;ff[0]=(ff[0]+c[i]*g[0])%mod;
for(int j=1;j<n;++j){
g[j]=((f[j]-g[j-1])*ny%mod+mod)%mod;
ff[j]=(ff[j]+c[i]*g[j])%mod;
}
}
for(int i=0;i<n;++i)
cout<<(ff[i]+mod)%mod<<' ';
}
posted @   lzrG23  阅读(39)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 上周热点回顾(2.17-2.23)
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
点击右上角即可分享
微信分享提示