多项式快速插值学习小记

最新推荐文章于 2023-11-18 20:28:13 发布
最新推荐文章于 2023-11-18 20:28:13 发布
阅读量1k 收藏 1
点赞数 2
分类专栏： NTT 模板与算法 文章标签： NTT
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/liyizhixl/article/details/88872799
版权
  • 今天终于抽空把这个综（du）合（liu）知识点学了，心力交瘁……

多项式快速插值

  • 给出 n n n 个点 ( x i , y i ) (x_i,y_i) (xi​,yi​) ，要求一个次数为 n − 1 n-1 n−1 的多项式 F ( x ) F(x) F(x) 满足： F ( x i ) = y i F(x_i)=y_i F(xi​)=yi​

  • 显然这个多项式是唯一确定的。

  • 根据拉格朗日插值法，我们有： F ( x ) = ∑ i = 1 n ∏ j ̸ = i ( x − x j ) ∏ j ̸ = i ( x i − x j ) y i F(x)=\sum_{i=1}^{n}\frac{\prod_{j\not=i}(x-x_j)}{\prod_{j\not=i}(x_i-x_j)}y_i F(x)=i=1∑n​∏j̸​=i​(xi​−xj​)∏j̸​=i​(x−xj​)​yi​

  • 这样我们是可以 O ( n 2 ) O(n^2) O(n2) 求的，考虑优化。

  • 我们先考虑对于每个 i i i ，如何快速得到 ∏ j ̸ = i ( x i − x j ) \prod_{j\not=i}(x_i-x_j) ∏j̸​=i​(xi​−xj​) 。

  • 设 M ( x ) = ∏ i = 1 n ( x − x i ) M(x)=\prod_{i=1}^{n}(x-x_i) M(x)=∏i=1n​(x−xi​) ，我们即需求： M ( x ) x − x i \frac{M(x)}{x-x_i} x−xi​M(x)​

  • 根据洛必达法则，当 x x x 取 x i x_i xi​ 时，分子分母都等于0，一个0/0型，上下求导得： l i m x → x i M ( x ) x − x i = M ′ ( x ) lim_{x→x_i}\frac{M(x)}{x-x_i}=M'(x) limx→xi​​x−xi​M(x)​=M′(x)

  • 于是我们先分治NTT求出 M ( x ) M(x) M(x) ，再求导得到 M ′ ( x ) M'(x) M′(x) ，之后将 x 1 − n x_{1-n} x1−n​ 代入多点求值即可！

  • 这一部分的复杂度是 O ( n   l o g 2 n ) O(n\ log^2n) O(n log2n) 的。

  • 求导、多点求值这些前置知识我的博客里都有讲：

    多项式的求逆、取模和多点求值学习小记

    多项式的ln、exp和快速幂学习小记

  • 又设 V i = y i ∏ j ̸ = i ( x i − x j ) V_i=\frac{y_i}{\prod_{j\not=i}(x_i-x_j)} Vi​=∏j̸​=i​(xi​−xj​)yi​​ ，则此时 V i V_i Vi​ 已知，我们要求： F ( x ) = ∑ i = 1 n V i ∏ j ̸ = i ( x − x j ) F(x)=\sum_{i=1}^{n}V_i\prod_{j\not=i}(x-x_j) F(x)=i=1∑n​Vi​j̸​=i∏​(x−xj​)

  • 还是分治NTT，设 L ( x ) = ∑ i = 1 n / 2 ( x − x i ) L(x)=\sum_{i=1}^{n/2}(x-x_i) L(x)=∑i=1n/2​(x−xi​) ， R ( x ) = ∑ i = n / 2 + 1 n ( x − x i ) R(x)=\sum_{i=n/2+1}^{n}(x-x_i) R(x)=∑i=n/2+1n​(x−xi​) ，则有： F ( x ) = ∑ i = 1 n / 2 V i ∏ j ̸ = i , 1 ≤ j ≤ n / 2 ( x − x j ) R ( x ) + ∑ i = n / 2 + 1 n V i ∏ j ̸ = i , n / 2 + 1 ≤ j ≤ n ( x − x j ) L ( x ) F(x)=\sum_{i=1}^{n/2}V_i\prod_{j\not=i,1\leq j\leq n/2}(x-x_j)R(x)+\sum_{i=n/2+1}^{n}V_i\prod_{j\not=i,n/2+1\leq j\leq n}(x-x_j)L(x) F(x)=i=1∑n/2​Vi​j̸​=i,1≤j≤n/2∏​(x−xj​)R(x)+i=n/2+1∑n​Vi​j̸​=i,n/2+1≤j≤n∏​(x−xj​)L(x)

  • 递归即可求得，递归底层的 F ( x ) F(x) F(x) 就是 V i V_i Vi​ 。

  • 还有就是这里的 L ( x ) 、 R ( x ) L(x)、R(x) L(x)、R(x) 在多点求值中已经算过了，不用再算一遍啦。

  • 这一部分的复杂度也是 O ( n   l o g 2 n ) O(n\ log^2n) O(n log2n) 。

  • 故总时间复杂度即为 O ( n   l o g 2 n ) O(n\ log^2n) O(n log2n) ，常数很大很大。

  • 模板题：洛谷 P5158 【模板】多项式快速插值

Code

#include<cstdio>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long LL;
const int N=1e5+5,M=18,G=3,mo=998244353;
int tot;
int xx[N],yy[N],val[N];
int a[N],b[N],c[N],rr[N];//a=b*c+rr
int ra[N],rb[N<<2],irb[N<<2];
int f[N*M<<1],stf[N<<2],enf[N<<2];
int g[N*M<<1],stg[N<<2],eng[N<<2];
int h[N],sth[N],enh[N],f3[N];
int f1[N<<2],f2[N<<2],wn[N<<2],rev[N<<2];
inline int read()
{
	int X=0,w=0; char ch=0;
	while(!isdigit(ch)) w|=ch=='-',ch=getchar();
	while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
	return w?-X:X;
}
void write(int x)
{
	if(x>9) write(x/10);
	putchar(x%10+'0');
}
inline int ksm(int x,int y)
{
	int s=1;
	while(y)
	{
		if(y&1) s=(LL)s*x%mo;
		x=(LL)x*x%mo;
		y>>=1;
	}
	return s;
}
inline void NTT(int *y,int len,int ff)
{
	for(int i=0;i<len;i++)
		if(i<rev[i]) swap(y[i],y[rev[i]]);
	for(int h=2,d=len>>1;h<=len;h<<=1,d>>=1)
		for(int i=0,k=h>>1;i<len;i+=h)
			for(int j=0,cnt=0;j<k;j++,cnt+=d)
			{
				int u=y[i+j],t=(LL)wn[cnt]*y[i+j+k]%mo;
				y[i+j]=u+t>=mo?u+t-mo:u+t;
				y[i+j+k]=u-t<0?u-t+mo:u-t;
			}
	if(ff==-1)
	{
		reverse(y+1,y+len);
		int inv=ksm(len,mo-2);
		for(int i=0;i<len;i++) y[i]=(LL)y[i]*inv%mo;
	}
}
void make(int v,int l,int r)
{
	if(l==r)
	{
		g[stg[v]=++tot]=mo-xx[l];
		g[eng[v]=++tot]=1;
		return;
	}
	int mid=l+r>>1,ls=v<<1,rs=ls|1;
	make(ls,l,mid),make(rs,mid+1,r);
	int na=eng[ls]-stg[ls]+1,nb=eng[rs]-stg[rs]+1;
	int len=1,ll=0;
	while(len<na+nb) len<<=1,ll++;
	for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<ll-1;
	int w0=ksm(G,(mo-1)/len);
	for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
	for(int i=0;i<na;i++) f1[i]=g[stg[ls]+i];
	for(int i=na;i<len;i++) f1[i]=0;
	for(int i=0;i<nb;i++) f2[i]=g[stg[rs]+i];
	for(int i=nb;i<len;i++) f2[i]=0;
	NTT(f1,len,1),NTT(f2,len,1);
	for(int i=0;i<len;i++) f1[i]=(LL)f1[i]*f2[i]%mo;
	NTT(f1,len,-1);
	stg[v]=tot+1;
	na+=nb-1;
	for(int i=0;i<na;i++) g[++tot]=f1[i];
	eng[v]=tot;
}
void getinv(int len,int ll)
{
	if(len==1)
	{
		irb[0]=ksm(rb[0],mo-2);
		return;
	}
	getinv(len>>1,ll-1);
	for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<ll-1;
	int w0=ksm(G,(mo-1)/len);
	for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
	for(int i=0;i<len>>1;i++) f1[i]=rb[i];
	for(int i=len>>1;i<len;i++) f1[i]=0;
	NTT(f1,len,1),NTT(irb,len,1);
	for(int i=0;i<len;i++) irb[i]=(2-(LL)f1[i]*irb[i]%mo+mo)*irb[i]%mo;
	NTT(irb,len,-1);
	for(int i=len>>1;i<len;i++) irb[i]=0;
}
void solve(int v,int l,int r,int fa)
{
	int na=enf[fa]-stf[fa],nb=eng[v]-stg[v];
	if(na>=nb)
	{
		int nc=na-nb;
		for(int i=0;i<=na;i++) a[i]=f[stf[fa]+i];
		for(int i=0;i<=nb;i++) b[i]=g[stg[v]+i];
		for(int i=0;i<=nc;i++) ra[i]=a[na-i];
		for(int i=0;i<=nb;i++) rb[i]=b[nb-i];
		for(int i=nc+1;i<=nb;i++) rb[i]=0;
		int len=1,ll=0;
		while(len<=nc*2+1) len<<=1,ll++;
		for(int i=nb+1;i<len;i++) rb[i]=0;
		for(int i=0;i<len;i++) irb[i]=0,f1[i]=0;
		getinv(len,ll);
		for(int i=0;i<=nc;i++) f1[i]=ra[i],f2[i]=irb[i];
		for(int i=nc+1;i<len;i++) f1[i]=f2[i]=0;
		NTT(f1,len,1),NTT(f2,len,1);
		for(int i=0;i<len;i++) f1[i]=(LL)f1[i]*f2[i]%mo;
		NTT(f1,len,-1);
		for(int i=0;i<=nc;i++) c[nc-i]=f1[i];
		for(int i=nc+1;i<nb;i++) c[i]=0;
		len=1,ll=0;
		while(len<nb<<1) len<<=1,ll++;
		for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<ll-1;
		int w0=ksm(G,(mo-1)/len);
		for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
		for(int i=0;i<nb;i++) f1[i]=b[i],f2[i]=c[i];
		for(int i=nb;i<len;i++) f1[i]=0,f2[i]=0;
		NTT(f1,len,1),NTT(f2,len,1);
		for(int i=0;i<len;i++) f1[i]=(LL)f1[i]*f2[i]%mo;
		NTT(f1,len,-1);
		for(int i=0;i<nb;i++) rr[i]=(a[i]-f1[i]+mo)%mo;
		while(nb>1 && !rr[nb-1]) nb--;
		stf[v]=tot+1;
		for(int i=0;i<nb;i++) f[++tot]=rr[i];
		enf[v]=tot;
	}else
	{
		stf[v]=tot+1;
		for(int i=stf[fa];i<=enf[fa];i++) f[++tot]=f[i];
		enf[v]=tot;
	}
	if(l==r)
	{
		val[l]=f[stf[v]];
		return;
	}
	int mid=l+r>>1;
	solve(v<<1,l,mid,v);
	solve(v<<1|1,mid+1,r,v);
}
void work(int v,int l,int r)
{
	if(l==r) return;
	int mid=l+r>>1,ls=v<<1,rs=ls|1;
	work(ls,l,mid);
	work(rs,mid+1,r);
	int na=enh[l]-sth[l]+1,nb=eng[rs]-stg[rs]+1;
	int len=1,ll=0;
	while(len<na+nb) len<<=1,ll++;
	for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<ll-1;
	int w0=ksm(G,(mo-1)/len);
	for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;
	for(int i=0;i<na;i++) f1[i]=h[sth[l]+i];
	for(int i=na;i<len;i++) f1[i]=0;
	for(int i=0;i<nb;i++) f2[i]=g[stg[rs]+i];
	for(int i=nb;i<len;i++) f2[i]=0;
	NTT(f1,len,1),NTT(f2,len,1);
	for(int i=0;i<len;i++) f1[i]=(LL)f1[i]*f2[i]%mo;
	NTT(f1,len,-1);
	na+=nb-1;
	for(int i=0;i<na;i++) f3[i]=f1[i];
	for(int i=na;i<len;i++) f3[i]=0;
	
	na=enh[mid+1]-sth[mid+1]+1,nb=eng[ls]-stg[ls]+1;
	for(int i=0;i<na;i++) f1[i]=h[sth[mid+1]+i];
	for(int i=na;i<len;i++) f1[i]=0;
	for(int i=0;i<nb;i++) f2[i]=g[stg[ls]+i];
	for(int i=nb;i<len;i++) f2[i]=0;
	NTT(f1,len,1),NTT(f2,len,1);
	for(int i=0;i<len;i++) f1[i]=(LL)f1[i]*f2[i]%mo;
	NTT(f1,len,-1);
	na+=nb-1;
	for(int i=0;i<na;i++) h[sth[l]+i]=f3[i]+f1[i]>=mo?f3[i]+f1[i]-mo:f3[i]+f1[i];
	enh[l]=sth[l]+na-1;
}
int main()
{
	int n=read();
	for(int i=1;i<=n;i++) xx[i]=read(),yy[i]=read();
	make(1,1,n);
	int m=eng[1]-stg[1];
	for(int i=0;i<=m;i++) f1[i]=g[stg[1]+i];
	for(int i=0;i<m;i++) f1[i]=(LL)f1[i+1]*(i+1)%mo;
	f1[m--]=0;
	stf[tot=0]=1;
	for(int i=0;i<=m;i++) f[enf[0]=++tot]=f1[i];
	solve(1,1,n,0);
	for(int i=1;i<=n;i++) val[i]=(LL)yy[i]*ksm(val[i],mo-2)%mo;
	tot=0;
	for(int i=1;i<=n;i++)
	{
		sth[i]=enh[i]=++tot;
		h[tot]=val[i];
	}
	work(1,1,n);
	for(int i=0;i<n;i++) write(h[sth[1]+i]),putchar(' ');
	return 0;
}
Felix-Lee
关注
  • 2
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
拉格朗日插值(多项式快速插值)(多项式全家桶)
coding_farmer_ac的博客
01-03 424
最近在学多项式,想起来一个月前写的一道拉格朗日,特来水贴。
插值（一）——多项式插值（C++）
最新发布
m0_50793313的博客
02-08 1464
采用C++实现多项式插值，这篇文章只是做一个原理探讨，实际应用并不会使用这种方法进行插值。
【MATLAB】史上最全的5种数据插值算法全家桶
Lwcah的博客
11-18 4070
有意向获取代码，请转文末观看代码获取方式~大家吃一顿火锅的价格便可以拥有5种数据插值算法，绝对不亏，知识付费是现今时代的趋势，而且都是我精心制作的教程，有问题可随时反馈~也可单独获取某一算法的代码（见每一算法介绍后文）~
专业数据处理--如何根据已有储层数据实现多层数据的快速插值
hhue2007的专栏
03-24 250
大家在做储层数据可视化分析与展示过程中，由于人工解释的层位有限，导致显示不那么精细逼真，一般情况下，研究人员主要解释有代表性的目的层，层与层之间的地质情况一般通过专业软件的插值功能或算法来实现，最常见的插值方法有线性插值、邻近插值等，当解释的层位数量达到一定程度后，后续的动态三维可视化显示就相对连续和逼真了。因此工作中就需要通过已知两层或多层数据，通过各种插值算法实现中间地段的快速插值计算来获取数据，有时会插入一层，有时会插入多层，大部分专业软件只能一次插入一层，最终只要满足研究需要，效果达到预期。
知识点 - 多项式插值法
qq_41848675的博客
08-12 4759
知识点 - 多项式插值法 解决问题类型： 已知f(0),f(1),f(2)…f(n)，求一个次数界为 nnn 的多项式，满足这些取值。 求nkn^knk 的前缀和，或nkn^knk 的k阶前缀和公式（直接插值，不过会T。对于第一个问题，后面给出了组合数解法。） 公式blog 讲义 2.多项式插值算法 2.1 多项式插值的存在唯一性 多项式一直以来备受数学家们青睐，一方面它构造起来简单，...
多项式插值与样条插值的解释与示例（matlab）
callmeup的博客
07-16 7177
插值介绍 信号插值，就是用已知点的测量值估计未知点的近似值。信号插值算法的应用范围有： 提高显示效果； 节省硬件成本，以软代硬； 减少远距离、大量数据通信的需要； 进行数据、图像解压缩 求解微分方程、积分方程； 计算函数值、零点、极值点、导数以及积分。 插值与拟合的异同点： 相同点：插值已知一些离散点，在一定约束下，求取定义在连续集合上的未知连续函数。 在图像上面不同：插值在图像上是一定得通过这些点，而拟合是拟合在图像上是逼近这些点。 在几何含义上不同：插值是找到一个(或几个分片光滑的)连续曲面来穿过
基于Matlab实现多项式与插值拟合案例（源码+数据）.rar
05-20
1、资源内容：基于Matlab实现多项式与插值拟合案例（源码+数据）.rar 2、适用人群：计算机，电子信息工程、数学等专业的大学生课程设计、期末大作业或毕业设计，作为“参考资料”使用。 3、解压说明：本资源需要电脑...
基于matlab实现泰勒级数展开多项式插值
01-06
对于一维函数y=2*x^3-3*x^2+2*x+0.5，在计算域等距选取4个点，取其解析解作为插值基点，构造基于泰勒级数展开式的多项式函数并通过所构造的函数计算(0,1)区间内等距分布的100个离散点的值 最后计算构造函数与原函数...
MATLAB多项式及插值PPT学习教案.pptx
10-07
MATLAB 多项式及插值学习教案 本篇学习教案旨在介绍 MATLAB...本篇学习教案旨在为读者提供一个系统、全面的 MATLAB 多项式及插值学习指南，帮助读者快速掌握 MATLAB 中的多项式及插值知识点，并将其应用于实际问题中。
七次多项式插值
11-16
七次多项式插值是一种数学方法，用于找到一个七次多项式函数，该函数精确地通过给定的八个离散点。在MATLAB环境中，这种插值可以通过使用内置的插值函数来实现，例如`polyfit`或`interp1`。这种技术在运动规划领域...
newton interpolation (牛顿多项式插值法)
12-06
使用matlab编写的牛顿多项式插值法，运行Run即可。代码中给了两个函数的实例插值。使用了matlab的符号计算。
基于Matlab实现多项式与插值拟合（源码+数据）.rar
10-03
1、资源内容：基于Matlab实现多项式与插值拟合（源码+数据）.rar 2、适用人群：计算机，电子信息工程、数学等专业的大学生课程设计、期末大作业或毕业设计，作为“参考资料”使用。 3、解压说明：本资源需要电脑端...
多项式插值问题
决明子的博客
03-23 2990
我写博客没别的意思，就是做个笔记，方便回头复习。 起因是，在看MOOC的西交大的算法课，这个问题看了N遍也没看懂，这我能忍？就查了点儿别的大佬的分享，简单总结一下，若侵权，请别纠。 这个问题的意思是：根据已知数据点的值，来预测未知数据点的值。 详细的说就是：已知n+1个点的值，就可以求出一个n次的插值函数p(x)，p(x)是接近未知原函数f(x)的函数，根据插值函数p(x)求出f(x)的未知...
洛谷P5158 【模板】多项式快速插值
weixin_34293911的博客
03-28 128
题面 传送门 前置芝士 拉格朗日插值，多项式多点求值 题解 首先根据拉格朗日插值公式我们可以暴力\(O(n^2)\)插出这个多项式，然而这显然是\(gg\)的 那么看看怎么优化，先来看一看拉格朗日插值的公式 \[f(x)=\sum_{i = 1}^{n} y_i \prod_{i \not = j} \frac{x - x_j}{x_i - x_j}\] 转化一下 \[f(x)=\sum_{i =...
FFT多项式快速幂，对于x^num取模，顺便再模一个998244353
INCINCIBLE的博客
02-05 3868
多项式快速幂 时间限制 : 60000 MS 空间限制 : 524288 KB问题描述：给一个n次多项式，求它的k次方。没关系，随手模一个998244353就行了。没关系，再随手模一个xm就行了。输入格式：第一行n，意义如上。 第二行n+1个数，a0,a1,…,an，分别是0,1,…,n次项系数。 第三行k，意义如上。 第四行m，意义如上。输出格式一行，b0,b1,…,bm-1，分别是0
【数值分析】2 - 插值法
知不足而奋进，望远山而前行
10-22 1353
满足插值条件的pxp(x)px是否存在且唯一？如满足插值条件的pxp(x)px存在，如何构造pxp(x)px？如何估计用pxp(x)px近似替代fxf(x)fx产生的误差？在介绍牛顿插值公式之前，我们需要介绍**差商（均差）**的定义：首先给出一阶差商的定义式：然后基于一阶差商，给出二阶差商的定义式（其实就是一阶差商的一阶差商）：最后，写出k1k+1k1阶差商的定义式：需要注意的是，kkk阶差商必须由k1k+1k1个节点构成，k。
快速三线性插值
dao_lian的博客
05-21 9265
转载自https://lianera.github.io/post/2018/fast-trilinear-interpolation/ 快速三线性插值 最近需要对一个体素进行插值，并且应用到一张大图像上。这个本来用三线性插值很容易就实现了，但是体素的尺寸很小，长宽高大概20x15x10的大小，图像的尺寸非常大，差不多是4000x3000，等我实现完普通的三线性插值后，发现算法速度实在太慢，单是这...
多点求值与快速插值
Johnny817的博客
01-29 701
多点求值 给出 nnn 次多项式 A(x)A(x)A(x) ，求出 mmm 个 xix_ixi​ 对应的 A(xi)A(x_i)A(xi​) 考虑分治，设 L(x)=∏i=1n2(x−xi)L(x)=\prod_{i=1}^{\frac{n}{2}}(x-x_i)L(x)=∏i=12n​​(x−xi​) ， R(x)=∏i=n2+1n(x−xi)R(x)=\prod_{i=\frac{n}{2}+...
数值计算之 插值法（3）多项式插值法的解，范德蒙矩阵，龙格现象
qq_41035283的博客
11-18 3391
数值计算之 插值法（3）多项式插值法的解，范德蒙矩阵，龙格现象前言多项式插值法的解与范德蒙矩阵龙格现象 前言 上两篇分别是拉格朗日插值法和牛顿插值法。 拉格朗日插值的思想是将一个多项式拆分为多个多项式，每个多项式完成一个节点插值。 牛顿插值法的思想是一种递推法，当出现新的抽样点时，增量更新一次插值函数。 从形式上来看，两种插值法得到的多项式不同，牛顿插值因为使用增量更新，因此计算量相对较小。 多项式插值法的解与范德蒙矩阵 未知表达式的函数f(x)f(x)f(x)满足y0=f(x0),y1=f(x1),y2=
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值