题解 nflsoj99 牛顿的烈焰激光剑（容斥，DP，数学）

题目大意

题目链接

给定(n), (k)。对于一个长度为(n)的正整数序列(a)，它是合法的当且仅当：

• (forall iin[1,n]:a_iin[1,k])。
• (forall 1leq i,jleq n,i eq j: a_i eq a_j)。

定义一个序列的权值就是所有元素的乘积。求所有长度为(n)的合法序列的权值之和。对一个给定的质数取模。

数据范围：(nleq 2000), (kleq 10^9)。

本题题解

因为序列里的数两两不同。我们先求严格上升的这些序列的权值之和，再乘以(n!)就是答案了。

DP。设(f[i])表示所有长度为(i)的，严格上升的序列，的权值和。

转移时，相当于在长度为(i-1)的序列后面新加入一个数。对新加的这个数两个要求：

1. 它不能在序列的前(i-1)位中出现过。
2. 它要严格大于序列的第(i-1)个数。

先只考虑第1个要求。转移完成后再把方案数除以(i)，就相当于满足了第2个要求。

对第1个要求，可以容斥。先令(f[i]=f[i-1]cdot(1+2+dots+k))，也就是说，我们认为第(i)个数可以选(1dots k)中任何一个数。但这样会把和前面某个位置上数相同的情况（不符合要求1）也算进去。所以要减去(f[i-2]cdot(1^2+2^2+dots +k^2))。然而这样又多减了一些情况（有三个数都相同的情况，这种情况之前没有计算，但也被减掉了）。所以要再加上(f[i-3]cdot (1^3+2^3+dots k^3))。以此类推，可知：

[f[i]=frac{sum_{j=1}^{i}(-1)^{j+1}f[i-j]cdot(1^j+2^j+dots+k^j)}{i} ]

其中，设(s[j]=1^j+2^j+dots+k^j) ((1leq jleq n))，也就是自然数幂求和，这理论上是可以(O(nlog n))求逆预处理出来的，但实际上由于模数不确定又没人想写任意模数NTT，所以直接暴力乘，就变成了(O(n^2))预处理（注意，根据主定理，把普通牛顿迭代求逆里的NTT换成暴力乘，时间复杂度是(O(n^2))而不是(O(n^2log n))）。

总时间复杂度(O(n^2))。

参考代码：

//problem:nflsoj99
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

const int MAXN=2005;
int MOD;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline void add(int& x,int y){x=mod1(x+y);}
inline void sub(int& x,int y){x=mod2(x-y);}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}

int fac[MAXN+5],ifac[MAXN+5],inv[MAXN+5];
inline int comb(int n,int k){
	if(n<k)return 0;
	return (ll)fac[n]*ifac[k]%MOD*ifac[n-k]%MOD;
}
void facinit(int lim=MAXN){
	fac[0]=1;
	for(int i=1;i<=lim;++i)fac[i]=(ll)fac[i-1]*i%MOD;
	ifac[lim]=pow_mod(fac[lim],MOD-2);
	for(int i=lim-1;i>=0;--i)ifac[i]=(ll)ifac[i+1]*(i+1)%MOD;
	for(int i=lim;i>=1;--i)inv[i]=(ll)ifac[i]*fac[i-1]%MOD;
}

namespace calcpow{
int f[MAXN+5],g[MAXN+5],h[MAXN+5];
void mul(int* a,const int* b,int n){
	//a=a*b
	//真*任意模数NTT
	static int res[MAXN+5];
	for(int i=0;i<n;++i){
		res[i]=0;
		for(int j=0;j<=i;++j){
			res[i]=((ll)res[i]+(ll)a[j]*b[i-j])%MOD;
		}
	}
	for(int i=0;i<n;++i)a[i]=res[i];
}
void get_inv(const int* a,int n,int* res){
	if(n==1){
		res[0]=pow_mod(a[0],MOD-2);
		return;
	}
	get_inv(a,(n+1)/2,res);
	static int tmp[MAXN+5];
	for(int i=0;i<n;++i)tmp[i]=res[i];
	mul(res,res,n);
	mul(res,a,n);
	for(int i=0;i<n;++i)res[i]=mod2(mod1(2*tmp[i])-res[i]);
}
void work(int n,int K,int* res){
	for(int i=0,t=n+1;i<=K;++i){
		f[i]=(ll)t*ifac[i+1]%MOD;
		g[i]=ifac[i+1];
		t=(ll)t*(n+1)%MOD;
	}
	get_inv(g,K+1,h);
	mul(f,h,K+1);
	for(int i=0;i<=K;++i)res[i]=(ll)f[i]*fac[i]%MOD;
}
}//namespace calcpow

int n,K,f[MAXN+5],s[MAXN+5];

int main() {
	cin>>K>>n>>MOD;
	facinit();
	calcpow::work(K,n,s);
	//for(int i=1;i<=n;++i)cout<<s[i]<<" ";cout<<endl;
	f[0]=1;
	for(int i=1;i<=n;++i){
		for(int j=1;j<=i;++j){
			//s[j]=1^j+2^j+...+K^j
			if(j%2==0)sub(f[i],(ll)f[i-j]*s[j]%MOD);
			else add(f[i],(ll)f[i-j]*s[j]%MOD);
		}
		f[i]=(ll)f[i]*inv[i]%MOD;
	}
	cout<<(ll)f[n]*fac[n]%MOD<<endl;
	return 0;
}