几道拉格朗日插值法题及杜教模板应用
Posted __Archer
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了几道拉格朗日插值法题及杜教模板应用相关的知识，希望对你有一定的参考价值。
今天做了几道拉格朗日插值的题，这个知识点虽然很早就知道，但是却很少遇到相关题目，今天记录一下相关套路。
首先是无敌的杜教模板：
#include <bits/stdc++.h> #define rep(ii,a,b) for(int ii=a;ii<=b;++ii) #define per(ii,a,b) for(int ii=b;ii>=a;--ii) using namespace std; typedef long long ll; const int maxn=1e5+10,maxm=2e6+10; const int INF=0x3f3f3f3f; const ll mod=1000000007; //head ll pow_mod(ll a,ll b,ll c=mod,ll ans=1){while(b){if(b&1) ans=(a*ans)%c;a=(a*a)%c,b>>=1;}return ans;} namespace polysum { // 注意该模板在计算0次多项式时可能有问题，需要特判!!!!! const int maxn=1010000; const ll mod=1000000007; ll a[maxn],f[maxn],g[maxn],p[maxn],p1[maxn],p2[maxn],b[maxn],h[maxn][2],C[maxn]; ll calcn(int d,ll *a,ll n) {//d次多项式(a[0-d])求第n项 if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d) { 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 maxm) {//初始化预处理阶乘和逆元(取模乘法) f[0]=f[1]=g[0]=g[1]=1; rep(i,2,maxm+4) f[i]=f[i-1]*i%mod; g[maxm+4]=pow_mod(f[maxm+4],mod-2); per(i,1,maxm+3) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *a,ll m) { // a[0].. a[m] \\sum_{i=0}^{n-1} a[i] // m次多项式求第n项前缀和 a[m+1]=calcn(m,a,m+1); rep(i,1,m+1) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,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=pow_mod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0; h[0][1]=1; rep(i,1,m+1) { 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+1) { 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=pow_mod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+1) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+1) C[i]=h[i][0]; ans=(calcn(m,C,n)*pow_mod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } }
第一题 : $cf 622 F. The Sum of the k-th Powers$
题意:求$\\sum\\limits_{i=1}^{n}i^k \\mod \\; 1000000007$
$solution:$设函数$f(n)=\\sum\\limits_{i=1}^{n}i^k$,这是一个$k+1$次多项式，直接插值求解即可。
#include <bits/stdc++.h> #define rep(ii,a,b) for(int ii=a;ii<=b;++ii) #define per(ii,a,b) for(int ii=b;ii>=a;--ii) using namespace std; typedef long long ll; const int maxn=1e6+10,maxm=2e6+10; const int INF=0x3f3f3f3f; const ll mod=1e9+7; //head int casn,n,m,k; int num[maxn]; ll a[maxn]; ll pow_mod(ll a,ll b,ll c=mod,ll ans=1){while(b){if(b&1) ans=(a*ans)%c;a=(a*a)%c,b>>=1;}return ans;} namespace polysum { // 注意该模板在计算0次多项式时可能有问题，需要特判!!!!! const int maxn=1010000; const ll mod=1e9+7; ll a[maxn],f[maxn],g[maxn],p[maxn],p1[maxn],p2[maxn],b[maxn],h[maxn][2],C[maxn]; ll calcn(int d,ll *a,ll n) {//d次多项式(a[0-d])求第n项 if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d) { 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 maxm) {//初始化预处理阶乘和逆元(取模乘法) f[0]=f[1]=g[0]=g[1]=1; rep(i,2,maxm+4) f[i]=f[i-1]*i%mod; g[maxm+4]=pow_mod(f[maxm+4],mod-2); per(i,1,maxm+3) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *a,ll m) { // a[0].. a[m] \\sum_{i=0}^{n-1} a[i] // m次多项式求第n项前缀和 a[m+1]=calcn(m,a,m+1); rep(i,1,m+1) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,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=pow_mod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0; h[0][1]=1; rep(i,1,m+1) { 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+1) { 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=pow_mod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+1) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+1) C[i]=h[i][0]; ans=(calcn(m,C,n)*pow_mod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } } int main() { ll n;int k;scanf("%lld%d",&n,&k); if(k==0) { printf("%lld\\n",n); return 0; } polysum::init(k+5); for(int i=1;i<=k;i++) a[i]=pow_mod(i,k,mod); ll ans = polysum::polysum(n+1,a,k); printf("%lld\\n",ans); return 0; }
第二题 ： $2019\\;icpc\\;$南昌邀请赛 $B. Polynomial$
$solution:$已知这是一个$n$次多项式，又给了$n+1$个点值，那么已经可以把这个多项式$f(x)$插值出来了，要求$\\sum\\limits_{i=l}^{r}f(i)$,设$g(x)=\\sum\\limits_{i=0}^{x}f(i)$,只需要计算$g(r)-g(l-1)$,可以看出$g$是一个$n+1$次多项式，所以再次插值即可
#include <bits/stdc++.h> #define rep(ii,a,b) for(int ii=a;ii<=b;++ii) #define per(ii,a,b) for(int ii=b;ii>=a;--ii) using namespace std; typedef long long ll; const int maxn=1e6+10,maxm=2e6+10; const int INF=0x3f3f3f3f; const ll mod=9999991; //head int casn,n,m,k; int num[maxn]; ll a[maxn],b[maxn]; ll pow_mod(ll a,ll b,ll c=mod,ll ans=1){while(b){if(b&1) ans=(a*ans)%c;a=(a*a)%c,b>>=1;}return ans;} namespace polysum { // 注意该模板在计算0次多项式时可能有问题，需要特判!!!!! const int maxn=1010000; const ll mod=9999991; ll a[maxn],f[maxn],g[maxn],p[maxn],p1[maxn],p2[maxn],b[maxn],h[maxn][2],C[maxn]; ll calcn(int d,ll *a,ll n) {//d次多项式(a[0-d])求第n项 if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d) { 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 maxm) {//初始化预处理阶乘和逆元(取模乘法) f[0]=f[1]=g[0]=g[1]=1; rep(i,2,maxm+4) f[i]=f[i-1]*i%mod; g[maxm+4]=pow_mod(f[maxm+4],mod-2); per(i,1,maxm+3) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *a,ll m) { // a[0].. a[m] \\sum_{i=0}^{n-1} a[i] // m次多项式求第n项前缀和 a[m+1]=calcn(m,a,m+1); rep(i,1,m+1) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,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=pow_mod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0; h[0][1]=1; rep(i,1,m+1) { 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+1) { 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=pow_mod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+1) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+1) C[i]=h[i][0]; ans=(calcn(m,C,n)*pow_mod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } } const int up = 9999990; int main() { int T;scanf("%d",&T); while(T--) { int n,m;scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lld",&a[i]); polysum::init(n+5); a[n+1] = polysum::calcn(n,a,n+1); //printf("fuck = %lld\\n",a[n+1]); b[0]=a[0]; for(int i=1;i<=n+1;i++) b[i] = (b[i-1]+a[i])%mod; for(int i=1;i<=m;i++) { int l,r;scanf("%d%d",&l,&r); ll ans = polysum::calcn(n+1,b,r) - polysum::calcn(n+1,b,l-1); //printf("r = %lld,l-1 = %lld\\n",polysum::calcn(n+1,b,r),polysum::calcn(n+1,b,l-1)); //ll ans = polysum::polysum(l,a,n) - polysum::polysum(r+1,a,n); ans = (ans%mod+mod)%mod; printf("%lld\\n",ans); } } return 0; }
第三题:$2019\\;icpc$南昌区域赛 $M$
$solution:$首先从偶数开始，连续四个数的异或和是$0$,然后发现要求的东西分为两部分:$\\sum\\limits_{k=1}^{t}\\sum\\limits_{i=\\frac{x+1}{2}}^{\\frac{y}{2}}(2*i)^k + \\sum\\limits_{k=1}^{t}\\sum\\limits_{i=x}^{y}[i^k\\mod\\;4\\;=2或0]$.
对前一部分，交换求和顺序后$\\sum\\limits_{i=\\frac{x+1}{2}}^{\\frac{y}{2}}\\frac{2*i-(2*i)^{t+1}}{1-2*i}$
设:$f(x) = \\sum\\limits_{i=0}^{x}\\frac{2*i-(2*i)^{t+1}}{1-2*i}$ ,f(x)是一个$t+2$次的多项式，插值即可。
#include <bits/stdc++.h> #define rep(ii,a,b) for(int ii=a;ii<=b;++ii) #define per(ii,a,b) for(int ii=b;ii>=a;--ii) using namespace std; typedef long long ll; const int maxn=1e5+10,maxm=2e6+10; const int INF=0x3f3f3f3f; const ll mod=1000000007; //head ll pow_mod(ll a,ll b,ll c=mod,ll ans=1){while(b){if(b&1) ans=(a*ans)%c;a=(a*a)%c,b>>=1;}return ans;} namespace polysum { // 注意该模板在计算0次多项式时可能有问题，需要特判!!!!! const int maxn=1010000; const ll mod=1000000007; ll a[maxn],f[maxn],g[maxn],p[maxn],p1[maxn],p2[maxn],b[maxn],h[maxn][2],C[maxn]; ll calcn(int d,ll *a,ll n) {//d次多项式(a[0-d])求第n项 if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d) { 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 maxm) {//初始化预处理阶乘和逆元(取模乘法) f[0]=f[1]=g[0]=g[1]=1; rep(i,2,maxm+4) f[i]=f[i-1]*i%mod; g[maxm+4]=pow_mod(f[maxm+4],mod-2); per(i,1,maxm+3) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *a,ll m) { // a[0].. a[m] \\sum_{i=0}^{n-1} a[i] // m次多项式求第n项前缀和 a[m+1]=calcn(m,a,m+1); rep(i,1,m+1) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,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=pow_mod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0; h[0][1]=1; rep(i,1,m+1) { 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+1) { 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=pow_mod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+1) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+1) C[i]=h[i][0]; ans=(calcn(m,C,n)*pow_mod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } } ll a[maxn]; ll cal(ll x,int y) { ll res = x/4; if(y==0) return res; res += (y<=(x%4)); res%=mod; return res; } //------------------------ ll x,y; ll solve(ll n,int t) { ll res=0; res = polysum::calcn(t+1,a,n); //printf("first = %lld\\n",polysum::calcn(t+1,a,n)); return res; } int main() { int t;scanf("%d%lld%lld",&t,&x,&y); a[0]=0; ll tmp=0; for(int i=1;i<=t+1;i++) { tmp = (2ll*i%mod - pow_mod(2ll*i,t+1,mod)+mod)%mod*pow_mod((1ll-2ll*i%mod+mod)%mod,mod-2,mod)%mod; a[i] = (a[i-1] + tmp) %mod; } polysum::init(maxn+10); ll ans = (solve(y/2,t) - solve((x+1)/2-1,t) +mod)%mod; ans = (ans + cal(y,1)*t%mod-cal(x-1,1)*t%mod+mod)%mod; ans = (ans + cal(y,2) - cal(x-1,2) + mod) %mod; ans = (ans + cal(y,3)*(t/2)%mod - cal(x-1,3)*(t/2)%mod + mod) %mod; // ll res=0; // res = (res + cal(y,1)*t%mod-cal(x-1,1)*t%mod+mod)%mod; // res = (res + cal(y,2) - cal(x-1,2) + mod) %mod; // res = (res + cal(y,3)*(t/2)%mod - cal(x-1,3)*(t/2)%mod + mod) %mod; //printf("second = %lld\\n",res); printf("%lld\\n",ans); return 0; }
第四题:$cf 955 F$
$solution:$这个$dp$转移方程好像很经典的，或许这一类转移最后的结果都是一个多项式的形式
#include <bits/stdc++.h> #define rep(ii,a,b) for(int ii=a;ii<=b;++ii) #define per(ii,a,b) for(int ii=b;ii>=a;--ii) using namespace std; typedef long long ll; const int maxn=3e3+10,maxm=2e6+10; const int INF=0x3f3f3f3f; const ll mod=1000000007; //head ll pow_mod(ll a,ll b,ll c=mod,ll ans=1){while(b){if(b&1) ans=(a*ans)%c;a=(a*a)%c,b>>=1;}return ans;} namespace polysum { // 注意该模板在计算0次多项式时可能有问题，需要特判!!!!! const int maxn=1010000; const ll mod=1000000007; ll a[maxn],f[maxn],g[maxn],p[maxn],p1[maxn],p2[maxn],b[maxn],h[maxn][2],C[maxn]; ll calcn(int d,ll *a,ll n) {//d次多项式(a[0-d])求第n项 if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d) { 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 maxm) {//初始化预处理阶乘和逆元(取模乘法) f[0]=f[1]=g[0]=g[1]=1; rep(i,2,maxm+4) f[i]=f[i-1]*i%mod; g[maxm+4]=pow_mod(f[maxm+4],mod-2); per(i,1,maxm+3) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *a,ll m) { // a[0].. a[m] \\sum_{i=0}^{n-1} a[i] // m次多项式求第n项前缀和 a[m+1]=calcn(m,a,m+1); rep(i,1,m+1) a[i]=(a[i-1以上是关于几道拉格朗日插值法题及杜教模板应用的主要内容，如果未能解决你的问题，请参考以下文章