学习笔记:拉格朗日插值

news/2025/2/7 7:46:44/文章来源:https://www.cnblogs.com/GGrun-sum/p/18701974

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

拉格朗日插值

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

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

定义:

对于 $ n $ 个点 $ (x_i,y_i) $ ,我们希望找到 $ f_i $ 满足:

\[f_i (j) = \begin{cases} y_i &j = x_i \\ 0 &j = x_k(k \not = i) \\ anything &otherwise \end{cases} \]

且 $ f_i $ 也是 $ n-1 $ 次的,那么我们可以得到多项式 $ F = \sum_i f_i $ ,这样 $ F $ 满足他是一个经过这 $ n $ 个点的 $ n-1 $ 次多项式。

那么对于每个 $ f_i $ ,我们尝试用因式构造,简单来说,为了满足 $ f_i(x_j) = 0 ,j \not = i $ ,我们使他含有 $ (x-x_j) , j \not = i $ 的因式,此时再去满足 $ f_i(x_i) = y_i $ ,我们先把 $ x_i $ 带进去,就会得到 $ \prod_{j,j \not = i} (x_i - x_j) $ ,这时我们给整个式子乘上 $ \large \frac{y_i}{\prod_{j,j \not = i} (x_i - x_j)} $ ,就满足了所有条件。

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

\[F = \sum_i y_i \prod_{j,j \not = i} \frac{ x - x_j }{ x_i - x_j } \]

例题:

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_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::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;
}

$ x_i $ 取值连续时的插值

我们观察式子

\[F = \sum_i y_i \prod_{j,j \not = i} \frac{ x - x_j }{ x_i - x_j } \]

对于 $ \prod x - x_j $ ,我们做一个前缀积和后缀积就可以 $ O(n) $ 解决,对于 $ \prod x_i - x_j $ 因为他们取值连续所以就变成了 $ k_1 ! k_2 ! $ 的形式,同样可以 $ 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_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::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 = \sum_i y_i \prod_{j,j \not = i} \frac{ x - x_j }{ x_i - x_j } \]

我们每加入一个点对于每个 $ i $ , $ \prod x - x_j $ 变化可以 $ O(1) $ 算, $ \prod x_i - x_j $ 也可以 $ O(1) $ 算。

插多项式系数

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

还是得改写式子:

首先对于整个式子

\[F = \sum_i y_i \prod_{j,j \not = i} \frac{ x - x_j }{ x_i - x_j } \]

有 $ y_i \prod_{j \not = i} {\frac{1}{x_i - x_j} } $ 是常数,姑且写为 $ c_i $ ,可以 $ O(n^2) $ 求出。对于 $ \prod_{j \not = i} x - x_j $ ,写为 $ \frac { \prod_j x - x_j } { x - x_i } $ ,我们把 $ \prod_j x - x_j $ 写为 $ m $ 。

现在式子长成了

\[F = \sum_i c_i \frac{m}{ x - x_i } \]

首先 $ c_i $ 是常数,不会对 $ x $ 次数造成影响,主要造成影响的是 $ m $ ,它包含 $ n $ 个形如 $ x - x_j $ 的式子,最后又除了一个 $ x - x_i $ ,学过二项式定理的应该知道 $ x $ 的次数就是由每个 $ x - x_j $ 的式子提供来的,考虑 dp ,设 $ dp_i $ 表示 $ x_i $ 的系数,那么每次考虑一个 $ x - x_j $ 的时候,要么选择 $ x $ ,有 $ dp_i \gets dp_{i-1} $ ,如果选择 $ x - x_i $ 的话,就有 $ dp_i = - x_i dp_i $ ,所以最终有 $ dp_i = dp_{i-1} - x_i dp_i $ 。

如果只考虑 $ m $ 的话,直接 dp 是 $ O(n^2) $ 的,但是有了 $ \frac{1}{x - x_i} $ 之后朴素实现就是 $ O(n^3) $ 的了,考虑有 $ \frac{1}{x - x_i} $ 之后相当于一个退背包的过程,所以可以对于每个 $ i $ 直接 $ O(n) $ 做,所以也是 $ O(n^2) $ 的。

#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_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::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<<' ';
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/879916.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

BOM最全基础信息:标准件、通用件、替换件、必选件

在生产制造领域,物料清单(BOM)是产品设计、生产计划和供应链管理的核心基础。本文系统梳理了BOM中各类零部件的分类方法,供大家参考。在生产制造的复杂领域中,我们会与各式各样的产品组成部分打交道。清晰、准确地对它们进行分类,并实施有效的管理,对于提升生产效率、保…

人工智能辅助芯片设计

芯片设计:一个近乎无限的问题空间 设计复杂性呈指数级增长 设计复杂性的含义 一连串棘手的问题 贯穿整个流程优化 HDL生成研究 使用GCN加速设计评估 人工智能辅助验证 参考文献链接https://www.hc2024.hotchips.org/assets/program/tutorials/3-HC24.synopsys.SteliosDiaman…

应用随机过程 | 期末 cheat sheet

出分后发布笔记……这篇博客汇总了「应用随机过程」2018 - 2022 的期末试题,并根据题型分类总结。 本站相关博客:应用随机过程 | 期末知识点总结特别鸣谢:知乎 | 九一居士 |《应用随机过程》课程笔记系列目录1 马尔可夫链计算题2 常返的马尔可夫链3 连续时间参数的马尔可夫链…

推荐《AI芯片开发核心技术详解》、《智能汽车传感器:原理设计应用》、《TVM编译器原理与实践》、《LLVM编译器原理与实践》4本书,非常感谢

4本书推荐《AI芯片开发核心技术详解》、《智能汽车传感器:原理设计应用》、《TVM编译器原理与实践》、《LLVM编译器原理与实践》由清华大学出版社资深编辑赵佳霓老师策划编辑的新书《AI芯片开发核心技术详解》已经出版,京东、淘宝天猫、当当等网上,相应陆陆续续可以购买。该…

OpenVX基本原理与历史

OpenVX基本原理 2.1 引言 2.1.1 摘要 OpenVX 是一个低级编程框架域,用于支持软件开发人员,可高效访问计算机视觉硬件加速功能和性能的可移植性。OpenVX 旨在支持现代硬件架构,例如,移动和嵌入式 SoC 以及桌面系统。其中许多系统是并行和异构的:多个处理器类型包括多核 CPU…

L4D2自制角色Mod - HUI篇

如何以相对简易的思路自制求生之路2求生者头像Mod本文是笔者尝试制作 求生之路2 角色 Mod 的过程中编写的笔记,笔者的背景是有基础的计算机知识和图像处理软件的使用经验,相信大多数读者朋友都有同样的水平。本文面向希望能快速简单地自定义游戏内角色图像/模型,但对更深层次…

使用Netty与前端请求进行交互实现实时通讯

引言因为不满足与一般的SpringBoot CRUD开发(太无聊了)所以去学一下网络编程,第一站就是通过B站老罗的EasyChat项目了解到了Netty这个网络框架,在学习这个项目之前也是去学习了一下Netty框架的使用以及相关的原理知识所以是有一定了解的,但是只是一味的学习不去实践总感觉是空中…

如何使用 Filebeat 8 连接 Easysearch

在日志场景,还是有很多小伙伴在使用 Filebeat 采集日志的。今天我来实战下使用 Filebeat 8 连接 Easysearch 。本次使用 Easysearch-1.9.0 版本和 Filebeat-8.17.0 版本做演示,也适用 Filebeat-oss-8.17.0 版本。 Easysearch 不开启兼容参数的情况 Easysearch 默认情况下未开…

25.2.7小记

异常捕捉 try{} catch{} import java.util.Scanner;public class ArrayIndex {public static void main(String[] args) {int [] a =new int[10];int idx;Scanner in = new Scanner(System.in);idx = in.nextInt();try{a[idx] = 10;System.out.println("hello");}cat…

《笨办法学Python3》PDF、EPUB免费下载

本书是一本Python入门书,适合对计算机了解不多,没有学过编程,但对编程感兴趣的读者学习使用。这本书以习题的方式引导读者一步一步学习编程,从简单的打印一直讲到完整项目的实现,让初学者从基础的编程技术入手,最终体验到软件开发的基本过程。本书是基于Python 3.6版本编…

RocketMQ实战—5.消息重复+乱序+延迟的处理

大纲 1.根据RocketMQ原理分析为什么会重复发优惠券 2.引入幂等性机制来保证数据不会重复 3.如何用死信队列处理优惠券系统数据库宕机 4.基于RocketMQ的订单库同步为什么会消息乱序 5.如何解决RocketMQ的消息乱序问题 6.RocketMQ的顺序消息机制的代码实现 7.基于RocketMQ的数据过…