拉格朗日插值法 在ACM竞赛中的一些应用

最新推荐文章于 2019-11-10 12:09:45 发布
最新推荐文章于 2019-11-10 12:09:45 发布
阅读量3.3k 收藏 24
点赞数 4
分类专栏： 数学
原文链接：https://riteme.github.io/blog/2017-3-18/lagrange-interpolation.html
版权

文章目录

什么是拉格朗日插值法

本文转载自 https://riteme.github.io/blog/2017-3-18/lagrange-interpolation.html
大佬写的太好了

这里写图片描述
这里写图片描述
当然求自然幂数和有很多种方法 Acdreamer

模板
namespace polysum {
    #define rep(i,a,n) for (int i=a;i<n;i++)
    #define per(i,a,n) for (int i=n-1;i>=a;i--)
    const int D=1e6+10;
    ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
    ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
    //..........................
   // 已知f(0),f(1)..f(d) 求 f(n)
    ll calcn(int d,ll *a,ll n) { // a[0].. a[d]  a[n] 
        if (n<=d) return a[n];
        p1[0]=p2[0]=1;
        rep(i,0,d+1) {
            ll t=(n-i+mod)%mod;
            p1[i+1]=p1[i]*t%mod;
        }
        rep(i,0,d+1) {
            ll t=(n-d+i+mod)%mod;
            p2[i+1]=p2[i]*t%mod;
        }
        ll ans=0;
        rep(i,0,d+1) {
            ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
            if ((d-i)&1) ans=(ans-t+mod)%mod;
            else ans=(ans+t)%mod;
        }
        return ans;
    }
    // 初始化，初始化的时候记得将D的值
    void init(int M) {
        f[0]=f[1]=g[0]=g[1]=1;
        rep(i,2,M+5) f[i]=f[i-1]*i%mod;
        g[M+4]=powmod(f[M+4],mod-2);
        per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;
    }
// 已知 f(0),f(1)...f(m),求  \sum_{i=0 }^{n} f[i]
   ll polysum(ll m,ll *a,ll n) { // a[0].. a[m] 
        ll b[D];
        ll b[D];
        for(int i=0;i<=m;i++) b[i]=a[i];
        b[m+1]=calcn(m,b,m+1);
        rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
        return calcn(m+1,b,n);// m次多项式的和是m+1 次多项式
    }
 
    ll qpolysum(ll R,ll n,ll *a,ll m) {
     // a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
        if (R==1) return polysum(n,a,m);
        a[m+1]=calcn(m,a,m+1);
        ll r=powmod(R,mod-2),p3=0,p4=0,c,ans;
        h[0][0]=0;h[0][1]=1;
        rep(i,1,m+2) {
            h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
            h[i][1]=h[i-1][1]*r%mod;
        }
        rep(i,0,m+2) {
            ll t=g[i]*g[m+1-i]%mod;
            if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
            else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
        }
        c=powmod(p4,mod-2)*(mod-p3)%mod;
        rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
        rep(i,0,m+2) C[i]=h[i][0];
        ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
} // polysum::init();
例题
1 CF622F

1这里写图片描述


LL qpow(LL a,LL b){
	LL ans = 1;
	while( b > 0){
		 if(b & 1) ans = ans*a%mod;
		 a = a*a%mod;
		 b >>= 1;
	}
	return ans;
}
namespace polysum {
    #define rep(i,a,n) for (int i=a;i<n;i++)
    #define per(i,a,n) for (int i=n-1;i>=a;i--)
    const int D=1e6+10;
    ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
    ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
    ll calcn(int d,ll *a,ll n) { // a[0].. a[d]  a[n]
        if (n<=d) return a[n];
        p1[0]=p2[0]=1;
        rep(i,0,d+1) {
            ll t=(n-i+mod)%mod;
            p1[i+1]=p1[i]*t%mod;
        }
        rep(i,0,d+1) {
            ll t=(n-d+i+mod)%mod;
            p2[i+1]=p2[i]*t%mod;
        }
        ll ans=0;
        rep(i,0,d+1) {
            ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
            if ((d-i)&1) ans=(ans-t+mod)%mod;
            else ans=(ans+t)%mod;
        }
        return ans;
    }
    void init(int M) {
        f[0]=f[1]=g[0]=g[1]=1;
        rep(i,2,M+5) f[i]=f[i-1]*i%mod;
        g[M+4]=powmod(f[M+4],mod-2);
        per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;
    }
// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
        ll b[D];
        for(int i=0;i<=m;i++) b[i]=a[i];
        b[m+1]=calcn(m,b,m+1);
        rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
        return calcn(m+1,b,n-1);
    }
    ll qpolysum(ll R,ll n,ll *a,ll m) {
     // a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
        if (R==1) return polysum(n,a,m);
        a[m+1]=calcn(m,a,m+1);
        ll r=powmod(R,mod-2),p3=0,p4=0,c,ans;
        h[0][0]=0;h[0][1]=1;
        rep(i,1,m+2) {
            h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
            h[i][1]=h[i-1][1]*r%mod;
        }
        rep(i,0,m+2) {
            ll t=g[i]*g[m+1-i]%mod;
            if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
            else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
        }
        c=powmod(p4,mod-2)*(mod-p3)%mod;
        rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
        rep(i,0,m+2) C[i]=h[i][0];
        ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
} // polysum::init();
const int maxn = 1e6+100;
LL a[maxn],b[maxn];
int main(void)
{


   LL n,k;
   cin>>n>>k;
   if(k == 0)
      return 0*printf("%I64d",n);
   polysum::init(k+100);
   b[0] = 0;
   for(int i = 1; i <= k+1; ++i)
      b[i] = qpow(i,k);
   LL ans = polysum::polysum(k,b,n+1);
   cout<<ans<<endl;
   return 0;
}
BZOJ2665
const int maxn = 500+10;
//LL qpow(LL a,LL b){LL s=1;while(b>0){if(b&1)s=s*a%mod;a=a*a%mod;b>>=1;}return s;}
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
LL f[maxn*2][maxn*2];
LL qpow(LL a,LL b,LL m){
	LL ans = 1;
	while(b > 0){
		if(b&1) ans = ans*a%m;
		a = a*a%m;
		b >>= 1;
	}
	return ans;
}
int main(void)
{
	
     LL A,n,d;
     cin>>A>>n>>d;
     me(f);
     f[0][0] = 1;
     for(int i = 1;i <= 2*n; ++i){
     	f[i][0] = 1;
     	for(int j = 1;j <=2 * n; ++j)
     	  f[i][j] = (f[i-1][j-1]*i%d*j+f[i-1][j])%d;
	 }
	 LL ans  = 0;
     if(A <= 2*n)
        ans = f[A][n];
     else{
     	for(LL i = 0;i <= 2*n; ++i){
     		LL tmp1 = 1,tmp2 = 1;
     		for(LL j = 0;j <= 2*n; ++j){
     			if(j != i)
     			tmp1 = tmp1*(A-j)%d,tmp2 = tmp2*(i-j)%d;
			 }
			ans = ans+f[i][n]*(tmp1*qpow(tmp2,d-2,d)%d)%d;
			ans = ans%d;
		 }
	 }
     cout<<(ans%d+d)%d<<endl;
     return 0;
}
BZOJ4559

参考博客

https://blog.csdn.net/cdsszjj/article/details/78778488

int R[maxn];
LL  U[maxn];
LL  polysum[maxn];
LL  dp[maxn][maxn],C[maxn][maxn];
LL  N,M,K;
LL qpow(LL a,LL b){
	LL ans = 1;
	while(b > 0){
		if(b & 1) ans =ans*a%mod;
		a = a*a%mod;
		b >>= 1;
	}
	return ans;
}
LL CC(LL x,LL y){
	if( x < 0||y < 0||x < y)
	   return 0;
	return C[x][y];
} 
LL poly(LL u, LL r){
	LL ans = 0; 
	int F[maxn];
	for(int i = 1;i <=  N + 1; ++i){
		F[i] = 0;
		for(int j = 1;j <= i; ++j){
			F[i] = (F[i] + qpow(i-j,r-1)*qpow(j,N-r)%mod)%mod;
		}
	    if(i == u)
	      return F[i];
	}
   	for(int i = 1;i <= N+1; ++i){

   		LL tmp1 = 1,tmp2 = 1;
   		for(int j = 1;j <= N+1; ++j){
   			if(i != j) tmp1 = tmp1*(u-j)%mod,tmp2 = tmp2*(i-j)%mod;
		   }
		assert(tmp1 != 0&&tmp2 != 0);
		ans = (ans + F[i]*tmp1%mod*qpow(tmp2,mod-2)%mod)%mod;
    }
	return ans;
}
int main(void)
{
	C[0][0] = 1;
   for(int i = 1;i < maxn; ++i){
   	C[i][0] = 1;
   	for(int j = 1;j < maxn; ++j)
        C[i][j] = (C[i-1][j-1]+C[i-1][j])%mod;
   }
  
   cin>>N>>M>>K;
//   cout<<N<<M<<K<<endl;
   for(int i = 1;i <= M; ++i)
     scanf("%lld",&U[i]);
   for(int i = 1;i <= M; ++i)
     scanf("%d",&R[i]);
   
   for(int i = 1;i <= M; ++i){
   	 // 拉格朗日算法求解 \sum_{x = 1}{x = U[i]}(U[i]-x)^(R-1)*x^(N-R) = poly[i] 
    polysum[i] = poly(U[i],R[i]);
   }
  
//   for(int i = 1;i <= M; ++i)
//      cout<<polysum[i]<<endl;
   dp[0][0] = 1;
   for(int i = 1;i <= M; ++i){
   	for(int j = 0;j <= N; ++j){
   		dp[i][j] = 0;
   		for(int w = 0;w <= j; ++w){
   			dp[i][j] = (dp[i][j] + polysum[i]*dp[i-1][w]%mod*CC(w,R[i]-j+w-1)%mod*CC(N-w-1,j-w)%mod)%mod; 
		   }
	   }
   }
   cout<<(dp[M][N-K-1]%mod+mod)%mod<<endl;
  

   return 0;
}
zzuzxy
关注
  • 4
    点赞
  • 踩
  • 24
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
拉格朗日插值算法及应用
Jon416的博客
03-30 1170
上学读书的时候对于数值分析的课程没有好好学过，直到工作了才发现会用到很多数值分析的算法，真应了那句话“出来混的，迟早要还的”，上学时玩的有多疯狂，工作时都要加倍补回来。这一篇先介绍拉格朗日插值算法，后边还有牛顿-辛普森法、最小二乘、极大似然、卡尔曼滤波……写这些文章的原因一方面是让自己掌握这些方法，二是能够让自己以后忘记时看到这些文章立马就能重新捡起这些知识，三是希望帮到一些需要了解这些知识的同学。插值是一种定义函数的过程，而该被定义的函数一定要在指定的点集上具有指定的函数值。
【原创】拉格朗日插值法学习笔记
致上
03-09 866
为你撰写的故事拉格朗日插值说在前面 拉格朗日插值 说在前面 为什么我的输入法打不起中文标点了啊,重启试一下. !@#$%^&amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;&amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;*()_|}{}&amp;amp;amp;amp;amp;amp;amp;amp;amp;quot;&amp;amp;amp;amp;amp;amp;amp;amp;am
拉格朗日插值法实现
bingshen的专栏
08-23 3367
拉格朗日插值法,代码,数学,算法
拉格朗日插值法
Never give in.
10-06 957
Text对于一个k次多项式函数f(x)=∑i=0kaixif(x)=\sum\limits_{i=0}^{k}a_ix^i有两种表示其的方法。 可以用传统的每一项的系数来表示。同时由于只要有K+1个(x,fx)的点对，就能唯一确定这个多项式，因此把这种表示法叫点值表示法。系数表示法转换成点值表示法叫求值，反过来就叫插值。当然可以高斯消元，但是复杂度比较大。主要讨论拉格朗日插值法。假设已经求得了K+
拉格朗日插法的应用(ACM，2019南昌邀请赛B题)
weixin_43958964的博客
07-30 1027
首先来了解拉格朗日插值法。 可理解为： 个点为(xi,yi)，这个多项式可以看做n个下面的函数的和。其中，第i个多项式函数fi除了在fi(xi)=yi外，其他的位置取值都为0，即fi(xj)=0(i≠j)。 f(x)=a(x−x1)(x−x2) 显然在f(x)=0 时，必有x1,x2 两解。推而广之，一个有n个顶点的n次多项式可以写为， f(x)=a(x−x1)(x−x2)……(x−...
牛客多校第一场 F. Sum of Maximum（拉格朗日插值）
Lee_w_j__的博客
07-20 1186
题目描述 Given a1, a2, ..., an, find  modulo (109+7). 输入描述: The input consists of several test cases and is terminated by end-of-file. The first line of each test case contains an integer n. The se...
缺失值处理：拉格朗日插值法
12-21
在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面...
【拉格朗日插值法】拉格朗日插值法介绍
04-11
拉格朗日插值法就是构造一个多项式，使得恰好在每一个x处取到对应的y 首先，n+1个点(xi,yi)若xi不同，则可以确定唯一的最高幂次不超过n的多项式。而如果题目直接或是间接给出了n+1个点，让我们求由这些点构成的...
人工智能/机器学习+拉格朗日插值法
02-17
多维插值：拉格朗日插值法目前主要用于一维情况，而在实际应用中，需要对多维数据进行插值。因此，未来的发展方向将是如何将拉格朗日插值法扩展到多维情况。 与深度学习的结合：拉格朗日插值法是一种传统的数值方法...
数学建模之拉格朗日插值法
10-04
拉格朗日插值法
缺失值处理：拉格朗日插值法.pdf
05-06
拉格朗日插值法是一种常见的缺失值处理方法，它可以通过已有的数据点， 来估计缺失数据点的值。其基本思想是假设缺失值的取值可以通过已有的数据点 进行插值来估计，具体的步骤如下： 1. 假设已有的数据点为$(x_1, y...
ACM基本算法
11-11
ACM基本算法
ACM算法大全
12-03
文档包含各种经典C语言算法，各种排序算法，搜索算法，图算法。。
ACM超级经典算法大集合
08-19
ACM超级经典算法大集合。大家都知道acm大赛对计算机专业的大学生说是很有用的
从拉格朗日插值法还原多项式
CRZbulabula的博客
03-10 1506
这里有份实现该功能的代码int Mul(const LL &x,const LL &y) {return x * y % mo;} int Add(const int &x,const int &y) {return (x + y) % mo;} int Dec(const int &x,const int &y) {return (x - y + mo) % mo;} void Work(in
NYOJ 题目178 找规律（Lagrange插值）
Jogging_Clown
02-26 864
找规律 时间限制：1000 ms  |  内存限制：65535 KB 难度：4 描述 大家一定见过这种题目：给你一些数请找出这些数之间的规律，写出下一个满足该规律的数。 比如：2 5 10 17 26，则可以看出这些数符合n*n+1这个通项公式，则下一个数为37。 这种通项公式不只一个，所以答案是不唯一的。但如果已知了N个数，且已知其通项公式是一个次数小
NYOJ 178 找规律(Lagrange插值公式)
herongwei 的 BLOG
02-24 2193
链接：click here 题意： 描述 大家一定见过这种题目：给你一些数请找出这些数之间的规律，写出下一个满足该规律的数。 比如：2 5 10 17 26，则可以看出这些数符合n*n+1这个通项公式，则下一个数为37。 这种通项公式不只一个，所以答案是不唯一的。但如果已知了N个数，且已知其通项公式是一个次数小于N的多项式，则答案就唯一确定了。 现在给你一个数列，请找出规律
拉格朗日插值法复习小计
Facico的博客
02-26 5197
拉格朗日插值法概述因为n次方的函数图像可以有n+1个点确定(比如说y=kx+b就只用两个点确定，y=ax2+bx+cax^2+bx+c，只用三个点就可以确定)，所以个n次方的方程给你n+1个x或者y就能确定这个方程。拉格朗日插值法就可以把这n+1个对应的值插进一个方程，然后给一个x，可以求出y。比如y=3x+2，给你(x=1,y=5;x=4,y=14)那么随意询问一个x，用拉格朗日插值法就可以求出对
拉格朗日与牛顿插值c++
mvpas的博客
11-10 283
#include using namespace std; typedef struct data{ double x; double y; }Data;//构造插值点 Data d[20];//插值点序列，最多20个 /拉格朗日/ double lagrange(int n, Data d,double z) { double result = 0;//结果，初始是0 for(int i=0;i...
拉格朗日插值法的应用
最新发布
09-12
拉格朗日插值法是一种用于数据插值的方法。它适用于需要根据给定的数据点来估计未知数据点的场景。这种方法的基本原理是使用一个多项式函数，通过将多个已知数据点带入该函数进行插值计算，从而得到一个拟合曲线。 拉格朗日插值法可以用于数据清洗中。例如，当我们有一组数据，并且其中某些数据点缺失或者存在异常值时，可以利用拉格朗日插值法来填充缺失值或校正异常值。这样可以更好地还原数据的真实情况，使得数据分析结果更加准确可靠。具体的应用过程可以参考相关文献中的代码示例和算法流程图。 总结起来，拉格朗日插值法是一种常用的数据插值方法，可以在数据分析和处理中广泛应用，特别适用于数据清洗和补全。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值