题解 nflsoj99 牛顿的烈焰激光剑（容斥，DP，数学）
题目大意
给定\(n\), \(k\)。对于一个长度为\(n\)的正整数序列\(a\)，它是合法的当且仅当：
- \(\forall i\in[1,n]:a_i\in[1,k]\)。
- \(\forall 1\leq i,j\leq n,i\neq j: a_i\neq a_j\)。
定义一个序列的权值就是所有元素的乘积。求所有长度为\(n\)的合法序列的权值之和。对一个给定的质数取模。
数据范围：\(n\leq 2000\), \(k\leq 10^9\)。
本题题解
因为序列里的数两两不同。我们先求严格上升的这些序列的权值之和，再乘以\(n!\)就是答案了。
DP。设\(f[i]\)表示所有长度为\(i\)的，严格上升的序列，的权值和。
转移时，相当于在长度为\(i-1\)的序列后面新加入一个数。对新加的这个数两个要求：
- 它不能在序列的前\(i-1\)位中出现过。
- 它要严格大于序列的第\(i-1\)个数。
先只考虑第1个要求。转移完成后再把方案数除以\(i\)，就相当于满足了第2个要求。
对第1个要求，可以容斥。先令\(f[i]=f[i-1]\cdot(1+2+\dots+k)\)，也就是说，我们认为第\(i\)个数可以选\(1\dots k\)中任何一个数。但这样会把和前面某个位置上数相同的情况（不符合要求1）也算进去。所以要减去\(f[i-2]\cdot(1^2+2^2+\dots +k^2)\)。然而这样又多减了一些情况（有三个数都相同的情况，这种情况之前没有计算，但也被减掉了）。所以要再加上\(f[i-3]\cdot (1^3+2^3+\dots k^3)\)。以此类推，可知：
其中，设\(s[j]=1^j+2^j+\dots+k^j\) (\(1\leq j\leq n\))，也就是自然数幂求和，这理论上是可以\(O(n\log n)\)求逆预处理出来的，但实际上由于模数不确定又没人想写任意模数NTT，所以直接暴力乘，就变成了\(O(n^2)\)预处理（注意，根据主定理，把普通牛顿迭代求逆里的NTT换成暴力乘，时间复杂度是\(O(n^2)\)而不是\(O(n^2\log 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;
}
