[BZOJ2655]calc(DP+拉格朗日插值)

[BZOJ2655]calc(DP+拉格朗日插值)

题面

一个序列a1,...,an是合法的，当且仅当：a1,...,an都是[1,A]中的整数且a1,...,an互不相等。
一个序列的值定义为它里面所有数的乘积，即a1a2...an。
求所有不同合法序列的值的和。两个序列不同当且仅当他们任意一位不一样。
输出答案对一个数mod取余的结果。

分析

设(dp_{i,j})表示前(j)个数,第(j)个数(leq i)的合法单调递增序列的值之和,那么转移方程为:
(dp_{i,j}=dp_{i-1,j-1}cdot i+dp_{i-1,j})

答案为(dp_{A,n}),暴力DP复杂度(O(nA)).不妨把(dp_{i,n})看成一个多项式,次数为(2i+1),用前面个点插出(A)处点值。注意(A<n)时有DP值可能为0,要去掉,因此DP的第二维要算到(3n)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#define maxn 3000
using namespace std;
typedef long long ll;
ll mod;
inline ll fast_pow(ll x,ll k){
	ll ans=1;
	while(k){
		if(k&1) ans=ans*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return ans;
}
inline ll inv(ll x){
	return fast_pow(x,mod-2);
}
ll fact(int n){
	ll ans=1;
	for(int i=1;i<=n;i++) ans=ans*i%mod;
	return ans;
}

ll lagrange(ll *x,ll *y,int n,ll k){
	ll ans=0;
	for(int i=1;i<=n;i++){
		ll up=y[i],dn=1;
		for(int j=1;j<=n;j++){
			if(j!=i){
				up=up*(k-x[j]+mod)%mod;//x[j]=j
				dn=dn*(x[i]-x[j]+mod)%mod;
			}
		}
		ans=(ans+up*inv(dn)%mod)%mod;
	} 
	ans=(ans+mod)%mod;
	return ans;
} 
ll A;
int n;
ll dp[maxn+5][maxn+5]; //dp[i][j]表示前j个数,第j个数<=i的递增序列值之和 
ll x[maxn+5],y[maxn+5];
int main(){
	scanf("%lld %d %lld",&A,&n,&mod);
	dp[0][0]=1;
	for(int i=1;i<=3*n;i++){
		dp[i][0]=1;
		for(int j=1;j<=n;j++) dp[i][j]=(dp[i-1][j]+dp[i-1][j-1]*i%mod)%mod;
	}
	int tot=0;
	for(int i=1;i<=3*n;i++){//dp[i][n]是一个2*n次多项式 
		if(dp[i][n]){
			x[++tot]=i;
			y[tot]=dp[i][n];//去掉A<n时可能等于0的dp值 
		}
	}
	printf("%lld
",lagrange(x,y,tot,A)*fact(n)%mod);//求dp[A][n],由于dp是递增序列,要乘上n! 
}