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

文章目录

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

1 CF622F
BZOJ2665
BZOJ4559

什么是拉格朗日插值法

本文转载自 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


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;
}

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

文章目录

什么是拉格朗日插值法

模板

例题

1 CF622F

BZOJ2665

BZOJ4559

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