拉格朗日插值法

最新推荐文章于 2022-11-30 10:13:34 发布

_shivering_

最新推荐文章于 2022-11-30 10:13:34 发布

阅读量878

点赞数 1

分类专栏：数学

本文链接：https://blog.csdn.net/shiveringkonnyaku/article/details/83014056

版权

数学专栏收录该内容

12 篇文章 0 订阅

订阅专栏

~~推公式利器~~

如果要求一个未知系数的n次多项式，必然要给出n+1个线性无关的方程。

把系数看成未知数，可以高斯消元；高斯消元是n^3的，对于比较大的数据就gg了。

然后就要开心地上拉格朗日插值法。。

仔细观察多项式，多项式也可以看成是一个函数，经过点(x0,y0)..(xn,yn)。

是否可以构造出这样的函数式呢？（可以证明，一旦构造出来，两者是等价的）

可以的。设多项式为L(x)。

我们期望当x=xi时，L(x)=yi；那么设L(x)=sum(wi(x)*yi)，则当x=xi时，wi(x)=1，wj(x)=0(i!=j)。

->wi(x)= $\prod (x-xj)/(xi-xj) (i!=j)$ 。

->L(x)= $\sum_{i=0}^{n}yi*\prod_{j=0,j!=i}^{n}((x-xj)/(xi-xj))；$ ；

这样显然地满足了n+1个线性无关的方程，L(x)就是所求的多项式。

然而，我们还希望它支持快速地插入一对(xn,yn)。

考虑变形。G(x)=(x-x0)*(x-x1)*...*(x-xn)。P(x,i)=(xi-x0)*(xi-x1)*..*(xi-xn); （没有(xi-xi)）

于是L(x)=sum(G(x)/((x-xi)*P(x,i))*yi)=G(x)*sum(inv((x-xi)*P(x,i))*yi)。

新插入一对(xn,yn)，则G(x)乘上(x-xn)->O(1)；

对于每个i!=n，P(x,i)乘上(xi-xn)->O(n)；

对于新插入的(xn,yn)，P(x,n)=(xn-x0)*..*(xn-x(n-1))；->O(n)；

于是求解一个n次多项式就变成了O(n^2)。

按照这个模拟就行了。（虽然不能直接按照写的存，但只要把不是x的部分存下来就可以代表它了）

这样虽然难以得到它展开后的系数，但是如果只要求对求出的多项式求值的话也足够了。

听说还有牛顿插值法。但是数学太差的我看得心态爆炸。。于是对此在这里没有展开。

放题目：

bzoj5339 TJOI2018教科书般的亵渎（求自然数幂和）

吐槽：模拟考时yy了一个m^2 log n分治（？）/倍增(?)不知道是啥的做法，比拉格朗日插值法多一个log。

非常伤心地T掉了。然后一看，果然是我太懒，写了递归方式，复杂度直接爆成了m^3~m^4。然后考后改成非递归。

跑得非常快。。。就以m^2 log n强行在bzoj上卡过了。其实这种做法优点在于式子直观好推（组合数一波），写起来方便，细节不多（拉格朗日要什么线性无关啊逆元啊之类的，存的东西也有点多，不熟悉的话会写很久），但是要求常数优秀且适用范围非常窄。如果时限不紧的话其实也可以尝试一下(逃。。

这是拉格朗日插值法：

#include<bits/stdc++.h>
using namespace std;
#define rep(x,y,z) for (int x=y; x<=z; x++)
#define downrep(x,y,z) for (int x=y; x>=z; x--)
#define ms(x,y,z) memset(x,y,sizeof(z))
#define LL long long
#define repedge(x,y) for (int x=hed[y]; ~x; x=edge[x].nex)
inline LL read(){
	LL x=0; int w=0; char ch=0;
	while (ch<'0' || ch>'9') w|=ch=='-',ch=getchar();
	while (ch>='0' && ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return w? -x:x;
}
const LL P=1e9+7;
const int maxm=55;
const int maxn=1e7+5; 
LL n,a[maxm],K;
int m;
LL pw(LL a,LL b){
	if (!b) return 1;
	LL s,res; s=a%P; res=1;
	for (LL i=1; i<=b; i<<=1){
		if (i&b) res=(res*s)%P;
		s=(s*s)%P;
	}
	return res;
}
void ex_gcd(LL a,LL b,LL c,LL &x,LL &y){
	if (!a){ x=0; y=c/b; return; }
	if (!b){ y=0; x=c/a; return; }
	ex_gcd(b,a%b,c,y,x);
	y-=(a/b)*x;
}
LL inv(LL t){
	LL x,y; ex_gcd(t,P,1,x,y);
	return (x%P+P)%P; 
}
LL x[maxm*2],y[maxm*2][maxm*2],w[maxm*2][maxm*2];
void pre(int k){
	x[k+2]=k+2; rep(i,1,k+2) y[k][i]=(y[k][i-1]+pw(x[i],k))%P;
	rep(i,1,k+1) w[k][i]=(w[k-1][i]*inv(x[i]-k-2))%P; 
	w[k][k+2]=1; rep(i,1,k+1) w[k][k+2]=(w[k][k+2]*inv(k+2-x[i]))%P;	 
}
LL solve(LL n,LL k){
	if (!n) return 0; 
	rep(i,1,k+2) if (!((n-x[i])%P)) return (y[k][i]%P+P)%P;
	LL ans=0; rep(i,1,k+2) ans=(ans+((inv((n-x[i])%P)*w[k][i])%P*y[k][i])%P)%P; 
    LL res=1; rep(i,1,k+2) res=(res*((n-x[i])%P))%P; 
	ans=(ans*res)%P; 
	return (ans%P+P)%P;
}
int main(){
	x[1]=1; x[2]=2; y[0][1]=1; y[0][2]=2; 
	w[0][1]=inv(-1); w[0][2]=inv(1); 
	rep(i,1,51) pre(i); 
	int cas=read(); 
	rep(cl,1,cas){
	   	n=read(); m=read(); K=m+1; 
	   	rep(i,1,m) a[i]=read(); sort(a+1,a+m+1);
	   	LL ans=0; LL nw=n; 
		rep(i,1,m){
			ans=(ans+solve(nw,K))%P; 
			rep(j,i,m) ans=(ans-pw(a[j]-n+nw,K))%P; 
		    nw-=a[i]-a[i-1];
		} 
		ans=(ans+solve(nw,K))%P;
		printf("%lld\n",(ans%P+P)%P);
	}   
	return 0;
}

这是倍增做法：

#include<bits/stdc++.h>
using namespace std;
#define rep(x,y,z) for (int x=y; x<=z; x++)
#define downrep(x,y,z) for (int x=y; x>=z; x--)
#define ms(x,y,z) memset(x,y,sizeof(z))
#define LL long long
#define repedge(x,y) for (int x=hed[y]; ~x; x=edge[x].nex)
inline LL read(){
	LL x=0; int w=0; char ch=0;
	while (ch<'0' || ch>'9') w|=ch=='-',ch=getchar();
	while (ch>='0' && ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return w? -x:x;
}
const LL P=1e9+7;
const int maxm=55;
const int maxn=1e7+5; 
LL n,a[maxm],K;
int m;
LL pw(LL a,LL b){
	if (!b) return 1;
	LL s,res; s=a%P; res=1;
	for (LL i=1; i<=b; i<<=1){
		if (i&b) res=(res*s)%P;
		s=(s*s)%P;
	}
	return res;
}
LL C[maxm*4][maxm*4],b[maxm],f[maxm][maxm],g[maxm][maxm];
LL solve(LL n,LL k){
	if (!n) return 0;
    int cnt=0; LL x=n;
    for (LL i=n; i; i>>=1) ++cnt;
    rep(i,1,cnt) { b[cnt-i+1]=x; x>>=1; } 
    rep(i,1,cnt) rep(j,0,k) g[i][j]=pw(b[i],j); 
	rep(i,1,cnt) rep(j,0,k){
    	f[i][j]=f[i-1][j]; 
		rep(p,0,j) { f[i][j]+=(((g[i-1][j-p]*C[j][p])%P)*f[i-1][p])%P; f[i][j]%=P; }
	    if (b[i]&1) { f[i][j]+=pw(b[i],j); f[i][j]%=P; } 
	} 
	return f[cnt][k]; 
}
int main(){
	rep(i,0,maxm) C[i][0]=1;
	rep(i,1,maxm) rep(j,1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
	int cas=read(); 
	rep(cl,1,cas){
	   	n=read(); m=read(); K=m+1; 
	   	rep(i,1,m) a[i]=read(); sort(a+1,a+m+1);
	   	LL ans=0; LL nw=n; 
		rep(i,1,m){
			ans=(ans+solve(nw,K))%P; 
			rep(j,i,m) ans=(ans-pw(a[j]-n+nw,K))%P;
		    nw-=a[i]-a[i-1];
		} 
		ans=(ans+solve(nw,K))%P;
		printf("%lld\n",(ans%P+P)%P);
	}   
	return 0;
}

听说此题有四五种做法，害怕。想找个时间补一波伯努利数的做法。留坑代填（不存在的。。）（逃

_shivering_

关注

1
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
拉格朗日插值法

推公式利器如果要求一个未知系数的n次多项式，必然要给出n+1个线性无关的方程。把系数看成未知数，可以高斯消元；高斯消元是n^3的，对于比较大的数据就gg了。然后就要开心地上拉格朗日插值法。。仔细观察多项式，多项式也可以看成是一个函数，经过点(x0,y0)..(xn,yn)。是否可以构造出这样的函数式呢？（可以证明，一旦构造出来，两者是等价的）可以的。设多项式为L(x)。...
复制链接

扫一扫