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

• • • • 什么是拉格朗日插值法
      • 模板
      • 例题

什么是拉格朗日插值法

本文转载自 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;}
    //..........................
   // 已知a_i 的d次多项式，求第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;
    }
// 已知a_i，并且知道a_i是m次多项式
   ll polysum(ll m,ll *a,ll n) { // a[0].. a[m] \sum_{i=0}^{n} a[i]
        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;
}