其他分享
首页 > 其他分享> > 几道拉格朗日插值法题及杜教模板应用

几道拉格朗日插值法题及杜教模板应用

作者:互联网

今天做了几道拉格朗日插值的题,这个知识点虽然很早就知道,但是却很少遇到相关题目,今天记录一下相关套路。
首先是无敌的杜教模板:

#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