几道拉格朗日插值法题及杜教模板应用
作者：互联网
今天做了几道拉格朗日插值的题，这个知识点虽然很早就知道，但是却很少遇到相关题目，今天记录一下相关套路。
首先是无敌的杜教模板：
#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;
}
}
View Code
第一题 : $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;
}
View Code
第二题 ： $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;
}
View Code
第三题:$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;
}
View Code
第四题:$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]+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 dp[maxn][maxn];
int n;
vector<int>e[maxn<<1];
void add(int u,int v) {
e[u].push_back(v);
}
void dfs(int u,int fa) {
for(int i=1;i<=n;i++) dp[u][i]=1;
for(auto v : e[u]) {
if(v==fa) continue;
dfs(v,u);
for(int j=1;j<=n;j++) dp[u][j] = dp[u][j]*dp[v][j]%mod;
}
for(int i=1;i<=n;i++) dp[u][i] = (dp[u][i] + dp[u][i-1])%mod;
}
int main() {
int D;scanf("%d%d",&n,&D);
for(int i=1;i<n;i++) {
int p;scanf("%d",&p);
add(p,i+1);
}
dfs(1,-1);
polysum::init(maxn+10);
//for(int i=1;i<=n;i++) printf("i = %d,%lld\n",i,dp[1][i]);
ll ans = polysum::calcn(n,dp[1],D);
printf("%lld\n",ans);
return 0;
}
View Code
第五题:$2018$牛客多校第一场$F$
标签：题及,const,插值法,ll,杜教,ii,int,maxn,mod 来源： https://www.cnblogs.com/ArcherLwz/p/14879481.html